初探mRNA、lncRNA联合分析之下游

生信菜鸟团

发布于 2023-08-23 08:55:54

5130

发布于 2023-08-23 08:55:54

文章被收录于专栏：生信菜鸟团

基因水平注释分类

虽然这个项目是在转录本水平上开展的研究，但既然我们拿到了基因表达矩阵，也干脆看一看一些基本情况，这个部分代码此处省略，基本上和后面的转录本水平对应代码，包括使用的封装函数，是一致的

DEGs结果

mRNA:

> deg_mrna <- deg(filter_count[mRNA_index,])
Warning: some variables in design formula are characters, converting to factorsestimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

  down normal     up 
   110  15061    380

lncRNA:

> deg_lncrna <- deg(filter_count[lncRNA_index,])
Warning: some variables in design formula are characters, converting to factorsestimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

  down normal     up 
    22   1913     43

可以发现基因水平和作者在转录水平得到的结果相差很多

联系后面我们做的转录本差异表达结果，可以验证我们前面一篇推文中的结论

明明PCA区分非常好，但是差异基因数量很少？ ... 那么转录本表达量差异分析和我们常见的基因有什么区别呢？ ... 这篇文章虽然作者用的不是stringtie但其从同一数据集获得的差异表达结果中转录本确实更少 ... 与RNA测序转录本差异分析相比，基因水平的差异表达分析更稳健，并且在实验上更可行。然而，使用基因计数进行统计分析可以掩盖转录物水平的动态

转录本水平注释分类

读取数据

rm(list = ls())
library(edgeR)
library(DESeq2)
library(FactoMineR)
library(factoextra)
library(org.Hs.eg.db)
library(stringr)
library(stringi)
library(tidyverse)
library(ggplot2)
library(patchwork)
library(pheatmap)
library(VennDiagram)
library(RColorBrewer)

exp <- read.table("6.exprs_mat/transcript_count_matrix.csv",header = T,sep = ",")
exp[1:4,1:4]

去重，低表达

exp <- na.omit(exp)
range(exp[,2:7])
exp$median <- apply(exp[2:7],1,median)
exp <- exp[order(exp$transcript_id,exp$median,decreasing = T),]

keep1 <- !duplicated(exp$transcript_id)
exp <- exp[keep1,]

keep2 <- rowSums(exp[,2:7]>0) >= floor(0.5*ncol(exp[,2:7]))
exp <- exp[keep2,]

table(duplicated(exp$transcript_id))
table(rowSums(exp[,2:7]>0))

rownames(exp) <- exp$transcript_id
exp <- exp[,c(-1,-8)]
exp <- exp[grepl("ENST",rownames(exp)),]
split_list <- rownames(exp) %>% str_split("\\.",simplify = T)
rownames(exp) <- split_list[,1]

group_list <-  c(rep("Ctrl",3),rep("ADS",3))
group_list
# table(group_list)

filter_count

filter_count <- exp
express_cpm <- log2(cpm(filter_count)+1)
range(express_cpm)

exp

boxplot(express_cpm, las=2)
exp <- normalizeBetweenArrays(express_cpm)
boxplot(exp, las=2)

注释mRNA和lncRNA

本文所有id的转换都是用biomaRt包进行

library(biomaRt)
mart<-useMart("ensembl")
data=listDatasets(mart)
data2=useDataset("hsapiens_gene_ensembl",mart=mart)
listFilters(data2) %>% View()

transcript_id <- rownames(exp)
# write.table(transcript_id,file = "transcript_id.txt",quote = F,sep = ",",col.names = F,row.names = F)
# magic needed
# annos <- read.table("annos.txt",sep = ",",header = F)

annos <- getBM(attributes=c('ensembl_transcript_id','biotype'),
       filters = 'ensembl_transcript_id', values = transcript_id,
       mart = data2)

table(annos$V2) %>% sort() %>% tail()

mRNA_index <- rownames(exp) %in% annos$V1[annos$V2=="protein_coding"]
table(mRNA_index)
lncRNA_index <- rownames(exp) %in% annos$V1[annos$V2=="lncRNA"]
table(lncRNA_index)

PCA和聚类热图

draw_pca <- function(exp){
  exp=t(exp)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
  exp=as.data.frame(exp)#将matrix转换为data.frame 
  library("FactoMineR")#画主成分分析图需要加载这两个包
  library("factoextra")  
  exp.pca <- PCA(exp , graph = FALSE)#现在exp最后一列是group_list，需要重新赋值给一个exp.pca,这个矩阵是不含有分组信息的
  pca_p<-fviz_pca_ind(exp.pca,
                   geom.ind = "point", # show points only (nbut not "text")
                   col.ind =  group_list, # color by groups 
                   addEllipses = T,
                   legend.title = "Groups"
  )
  return(pca_p)
}

p1 <- draw_pca(exp[mRNA_index,])
p2 <- draw_pca(exp[lncRNA_index,])
p1+p2

draw_heatp <- function(exp){
  cg=names(tail(sort(apply(exp,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
  library(pheatmap)
  pheatmap(exp[cg,],show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
  n=t(scale(t(exp[cg,]))) # 'scale'可以对log-ratio标准化数值进行归一化
  n[n>2]=2 
  n[n< -2]= -2
  n[1:4,1:4]
  pheatmap(n,show_colnames =F,show_rownames = F)
  
  # 绘制这些基因的热图，并将分组信息加入到热图中
  ac=data.frame(group=group_list)
  rownames(ac)=colnames(n)
  pheatmap(n,show_colnames =F,show_rownames = F,
           annotation_col=ac)
  pheat_p<-pheatmap(n,show_colnames =F,show_rownames = F,
               annotation_col=ac,
               filename = 'heatmap_top1000_sd.png')
  return(pheat_p)
}

library(ggplotify)
p1 <- draw_heatp(exp[mRNA_index,]) %>% as.ggplot()
p2 <- draw_heatp(exp[lncRNA_index,]) %>% as.ggplot()
p1+p2

mRNA、lncRNA

sd前1000转录本：

其实可以发现，我们的热图和作者一样，都有些奇怪，没有较好地聚在一起

但毕竟每组样本只有三个，PCA分组勉勉强强，还是接着往下做了

差异基因分析

deg <- function(filter_count){
  # 默认counts所以要求整数 内部有log2
  #1 查看分组信息和表达矩阵数据
  exprSet <- filter_count
  dim(exprSet)
  table(group_list)
  # 加载包
  library(DESeq2)
  #2 第一步，构建DESeq2的DESeq对象dds(建立DEseq数据矩阵)
  colData <- data.frame(row.names=colnames(exprSet),group_list=group_list)
  colData
  dds <- DESeqDataSetFromMatrix(countData = exprSet,
                                colData = colData,
                                design = ~ group_list)
  dds
  #3 第二步，进行差异表达分析
  dds2 <- DESeq(dds)
  #4 提取差异分析结果，trt组对untrt组的差异分析结果
  tmp <- results(dds2,contrast=c("group_list","ADS","Ctrl"))
  DEG_DESeq2 <- as.data.frame(tmp[order(tmp$padj),])
  head(DEG_DESeq2)  # 联系之前的内容DEseq2和edge会存在p值为0的情况
  # 去除差异分析结果中包含NA值的行
  DEG_DESeq2 = na.omit(DEG_DESeq2)
  # 筛选上下调差异基因，设定阈值
  fc_cutoff <-2.0
  fdr <- 0.05#p值
  DEG_DESeq2$regulated <- "normal"
  loc_up <- intersect(which(DEG_DESeq2$log2FoldChange>log2(fc_cutoff)),
                      which(DEG_DESeq2$pvalue<fdr))
  loc_down <- intersect(which(DEG_DESeq2$log2FoldChange< (-log2(fc_cutoff))),
                        which(DEG_DESeq2$pvalue<fdr))
  DEG_DESeq2$regulated[loc_up] <- "up"
  DEG_DESeq2$regulated[loc_down] <- "down"
  print(table(DEG_DESeq2$regulated))
  return(DEG_DESeq2)
}

deg_mrna <- deg(filter_count[mRNA_index,])
deg_lncrna <- deg(filter_count[lncRNA_index,])

火山图

library(EnhancedVolcano)
draw_volc <- function(deg_res){
EnhancedVolcano(deg_res,x='log2FoldChange',y='pvalue',
                lab=rownames(deg_res),
                pCutoff = 0.05,
                FCcutoff = log2(2))
}

p1 <- draw_volc(deg_mrna)
p2 <- draw_volc(deg_lncrna)
p1+p2

check

检查一下作者说的"Top differentially expressed"的转录本:

差别很大

但可能是由于这是在转录水平上

手动查看我们的转录本差异表达结果中的top，发现排第一的ENST00000343067其实也是EPB41

于是回到基因水平查看DEGs：

看看和原文差异表达结果logfc的散点图：

library(data.table)
deg_mrna_au <- fread("./author_supply/GSE209825_mRNA_ADS_vs_Normal.DE_ALL.annotaion.xls.gz")
deg_lncrna_au <- fread("./author_supply/GSE209825_lncRNA_ADS_vs_Normal.DE_ALL.annotaion.xls.gz")

range(deg_mrna_au$`log2(foldchange)`)
# 转换-Inf和Inf
# mRNA
topfc <- sort(deg_mrna_au$`log2(foldchange)`) %>% unique() %>% head()
tailfc <- sort(deg_mrna_au$`log2(foldchange)`) %>% unique() %>% tail()
deg_mrna_au$`log2(foldchange)`[deg_mrna_au$`log2(foldchange)`==Inf] <- tailfc[5]+1
deg_mrna_au$`log2(foldchange)`[deg_mrna_au$`log2(foldchange)`==-Inf] <- topfc[2]-1
# lncRNA
topfc <- sort(deg_lncrna_au$`log2(foldchange)`) %>% unique() %>% head()
tailfc <- sort(deg_lncrna_au$`log2(foldchange)`) %>% unique() %>% tail()
deg_lncrna_au$`log2(foldchange)`[deg_lncrna_au$`log2(foldchange)`==Inf] <- tailfc[5]+1
deg_lncrna_au$`log2(foldchange)`[deg_lncrna_au$`log2(foldchange)`==-Inf] <- topfc[2]-1

range(deg_mrna_au$`log2(foldchange)`)
range(deg_lncrna_au$`log2(foldchange)`)

overlap_mrna <- intersect(rownames(deg_mrna), deg_mrna_au$transcript_id)
deg_lncrna_au$transcript_id <- str_split(deg_lncrna_au$transcript_id,"\\.",simplify = T)[,1]
overlap_lncrna <- intersect(rownames(deg_lncrna), deg_lncrna_au$transcript_id)

draw_scatter <- function(overlap, deg_my, deg_au){
  logfc_my <-  deg_my$log2FoldChange[match(overlap, rownames(deg_my))]
  logfc_author <-  deg_au$`log2(foldchange)`[match(overlap, deg_au$transcript_id)]
  
  overlap_logfc <- data.frame(
    genes = overlap,
    logfc_my,
    logfc_author
  )
  
  # 计算相关系数
  cor_1 <- cor(logfc_author[logfc_author>0],logfc_my[logfc_author>0])
  cor_2 <- cor(logfc_author[logfc_author<0],logfc_my[logfc_author<0])
  
  # basic scatterplot
  scatter_plot <- ggplot(overlap_logfc, aes(x=logfc_my, y=logfc_author)) + 
      geom_point()
  scatter_plot+
    annotate("text",2,2,size=6,label=round(cor_1,3),color="red")+
    annotate("text",-2,-2,size=6,label=round(cor_2,3),color="red")
}

p1 <- draw_scatter(overlap_mrna, deg_mrna, deg_mrna_au)
p2 <- draw_scatter(overlap_lncrna, deg_lncrna, deg_lncrna_au)
p1+p2

看看我们常规featureCounts定量

NCBI pipeline（check）

NCBI最近使用统一的pipeline：hisat2+featureCounts，对数据库中的bulk转录组数据进行了处理并公开了所有项目的原始计数矩阵+归一化计数矩阵+基因注释，非常方便

详情参见：

https://www.ncbi.nlm.nih.gov/geo/info/rnaseqcounts.html 足不出户，GEO能进行RNA-seq差异表达分析啦（测试版）

我们这里获取本数据集的NCBI-generated RNA-seq count data

https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE209825

也进行我们上述流程看看

NCBI给的其实是基因水平的定量结果

可以发现数量关系和我们前面的结果是颠倒的：

于是手动检查EPB41这个基因：

发现确实是反的

我们前面主要是对热图存疑，这里就主要看一下热图：

也不咋地，所以这应该是实验样本或者测序数据本身的问题

寻找lncRNA共表达mRNA

This study analyzed the correlation between lncRNA and mRNA in the samples using the Pearson correlation coefficient method. The absolute values of correlation coefficients >0.95 and p < 0.001 were specified as screening criteria. Then, the biological functions of lncRNAs were predicted by performing a functional enrichment analysis on mRNAs.
In brief, the coexpressed mRNAs of each DE lncRNA were identified by the correlation with protein-coding gene expression, which was used for functional enrichment analysis. ...

mat_mrna <- filter_count[mRNA_index,]
mat_lncrna <- filter_count[lncRNA_index,]

up_lncrnas <- rownames(deg_lncrna)[deg_lncrna$regulated=="up"]
down_lncrnas <- rownames(deg_lncrna)[deg_lncrna$regulated=="down"]

get_co_mrna <- function(lncrnas){
  
  co_mrna <- c()
  
  for (i in seq(1,length(lncrnas))){
    
    for (j in seq(1, length(rownames(mat_mrna)))){
          res <- cor.test(unlist(mat_lncrna[lncrnas[i],]),
                    unlist(mat_mrna[j,]),method = "p")
          if (!is.na(res$estimate) & res$estimate>0.95 & res$p.value<0.001){
            co_mrna <- c(co_mrna, rownames(mat_mrna)[j])
          }
    }
    
  }
  
  return(co_mrna)
  
}

up_lncrnas_co_mrnas <- unique(get_co_mrna(up_lncrnas))
down_lncrnas_co_mrnas <- unique(get_co_mrna(down_lncrnas))

获取对应gene symbol

save(up_lncrnas_co_mrnas,down_lncrnas_co_mrnas,file = "co_mrna.Rdata")
# biomart magic
# load("co_genes_ENTREZID.Rdata")
up_co_genes <- getBM(attributes=c('ensembl_transcript_id','hgnc_symbol'),
                     filters = 'ensembl_transcript_id', up_lncrnas_co_mrnas,
                     mart = data2)
down_co_genes <- getBM(attributes=c('ensembl_transcript_id','hgnc_symbol'),
                       filters = 'ensembl_transcript_id', down_lncrnas_co_mrnas,
                       mart = data2)

up_genes <- unique(up_co_genes$entrezgene_id)
down_genes <- unique(down_co_genes$entrezgene_id)

富集分析

我们只看看lncRNA共表达mRNAmap到的基因富集分析情况，作者还对所有差异表达的mRNA进行了分析

library(clusterProfiler)
library(ggthemes)
library(org.Hs.eg.db)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)

1.GO 富集分析

up_go <- enrichGO(gene = up_genes,
                  OrgDb = org.Hs.eg.db,
                  ont = "ALL",
                  readable = TRUE)
down_go <- enrichGO(gene = down_genes,
                  OrgDb = org.Hs.eg.db,
                  ont = "ALL",
                  readable = TRUE)
save(up_go,down_go,file="GO.Rdata")

# 可视化
#条带图
p1<-barplot(up_go)
p2 <- barplot(down_go)
p1+p2

# 气泡图
p1 <- dotplot(up_go)
p2 <- dotplot(down_go)
p1+p2

p1<-dotplot(up_go, split = "ONTOLOGY",
            font.size = 10, showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") + 
  scale_y_discrete(labels = function(x) str_wrap(x, width = 45))
p2 <- dotplot(down_go, split = "ONTOLOGY",
            font.size = 10, showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") + 
  scale_y_discrete(labels = function(x) str_wrap(x, width = 45))
p1+p2

一个通路居然在上下调基因集里面都富集到了

GO:up&down

作者：

2.KEGG 通路富集分析

up_kegg <- enrichKEGG(up_genes, organism = "hsa")
down_kegg <- enrichKEGG(down_genes, organism = "hsa")
save(up_kegg,down_kegg,file = "KEGG.Rdata")

看看富集到了吗

table(up_kegg@result$p.adjust<0.05)
table(down_kegg@result$p.adjust<0.05)

p1 <- barplot(up_kegg,showCategory=10,title="KEGG Pathway",
              colorBy="p.adjust",font.size = 8,label_format = 15)+
  scale_fill_gradient(low="#DC0000B2",high ="#00A087B2" )
p2 <- barplot(down_kegg,showCategory=10,title="KEGG Pathway",
              colorBy="p.adjust",font.size = 8,label_format = 15)+
  scale_fill_gradient(low="#DC0000B2",high ="#00A087B2" )
p1+p2

KEGG:up&down

作者：

lncRNA–mRNA共表达网络

LncRNAs play a biological role by regulating mRNAs. A coexpression network of lncRNA/mRNA was constructed to investigate the potential interactions between lncRNA and mRNA, which could identify the key lncRNAs involved in ADS and their potential functional. This study analyzed the correlation between lncRNA and mRNA in the samples using the Pearson correlation coefficient method. The absolute values of correlation coefficients >0.95 and p < 0.001 were specified as screening criteria
Ultimately, eight interested lncRNAs were generated in the analysis, and a lncRNA-mRNA coexpression network was constructed for visualization

这里作者只根据前面那样计算的相关性确定共表达的mRNA和lncRNA构建网络（所以可以看到mRNA之间没有连线）

我们这里尝试用WGCNA寻找共表达模块（mRNA+lncRNA），看看结果能不能有overlap

参考：

4个发育时间点的总共12个鸡转录组测序样本的长非编码RNA的鉴定 WGCNA仅仅是划分基因模块，其它都是附加分析

WGCNA

library(tidyverse)
library(edgeR)
library(WGCNA)
library(FactoMineR)
library(factoextra)  
library(tidyverse) # ggplot2 stringer dplyr tidyr readr purrr  tibble forcats
library(data.table) #多核读取文件

这时需要把mRNA和lncRNA的标准化后的cpm矩阵放在一起：

express_cpm[1:4,1:4]
range(express_cpm)
group_list

PCA和聚类树参考前面，一组就三个样，懒得看了，还是不拿DEGs做WGCNA了

将近12W个转录本太大了，取子集:MAD前2W的转录本，进行

tmp <- as.data.frame(express_cpm)
tmp$mad <- apply(tmp,1,mad)
keep <- tmp %>% arrange(-mad) %>% head(20000) %>% rownames()

# 挑选最佳阈值power
R.sq_cutoff = 0.8  #设置R^2 cut-off
datExpr <- t(express_cpm[keep,]) %>% as.data.frame()
nGenes <- ncol(datExpr)
nSamples <- nrow(datExpr)
if(T){
  # Call the network topology analysis function
  #设置power参数选择范围 1到20，每次增加1，以及22到30，每次增加2
  # 为了在选择最佳power时，尽可能涵盖可能的取值范围
  powers <- c(seq(1,20,by = 1), seq(22,30,by = 2)) 
  # 用来选择最佳power参数的，
  # 输入参数包括表达数据datExpr、网络类型（这里是无向网络）、power取值的向量、R平方的截断值、以及是否输出过程信息
  # 输出结果为一个包含有关于不同power下拟合优度、均值连接度、平均k值等的信息表
  sft <- pickSoftThreshold(datExpr, 
                           networkType = "unsigned",
                           powerVector = powers, 
                           RsquaredCut = R.sq_cutoff,  
                           verbose = 5)
  #SFT.R.sq > 0.8 , slope ≈ -1
  # 将拟合优度图和均值连接度图保存至PDF文件中
  pdf("step2_power-value.pdf",width = 16,height = 12)
  # Plot the results: 寻找拐点，确认最终power取值
  par(mfrow = c(1,2));
  cex1 = 0.9;
  
  # Scale-free topology fit index as a function of the soft-thresholding power
  # 画出在不同power下拟合优度的趋势图，x轴为power，y轴为拟合优度（signed R2值）
  # 拟合优度要乘以正负1，这是因为负的拟合优度意味着模型拟合得更差
  plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
       xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n")
  # 在每个数据点上添加power的数值标注，颜色为红色
  text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
       labels=powers,cex=cex1,col="red")
  # this line corresponds to using an R^2 cut-off of h
  # 画一条水平的线，标注R平方的截断点，颜色为红色
  abline(h=R.sq_cutoff ,col="red")
  
  # Mean connectivity as a function of the soft-thresholding power
  # 画出在不同power下均值连接度的趋势图，x轴为power，y轴为均值连接度
  plot(sft$fitIndices[,1], sft$fitIndices[,5],
       xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n")
  # 在每个数据点上添加power的数值标注，颜色为红色
  text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
  # 画一条水平的线，标注均值连接度的截断点，均值连接度大于100的边自动被过滤掉，颜色为红色
  abline(h=100,col="red")
  dev.off()
}  # 可看出拐点大致在power为16时出现，且各项数值基本满足挑选标准，因此设定power=16

# 打印出在不同power下估计的最佳power取值
sft$powerEstimate  #查看估计的最佳power [1] NA
power = sft$powerEstimate

# 若无向网络在power小于15或有向网络power小于30内，没有一个power值使
# 无标度网络图谱结构R^2达到0.8且平均连接度在100以下，可能是由于
# 部分样品与其他样品差别太大。这可能由批次效应、样品异质性或实验条件对
# 表达影响太大等造成。可以通过绘制样品聚类查看分组信息和有无异常样品。
# 如果这确实是由有意义的生物变化引起的，也可以使用下面的经验power值。

# 如果最佳power取值为空
# 根据数据样本数确定经验power取值：
# 如果样本数小于20，则unsigned类型网格选择power=9，否则选择power=等一系列取值；
# 如果样本数大于等于20，则unsigned类型网格选择power=6，否则选择power=12。这个默认给出较为通用的power取值方案，但是具体的数据情况下是否适用需再加以考虑。
if(is.na(power)){
  # 官方推荐 "signed" 或 "signed hybrid"
  # 为与原文档一致，故未修改
  type = "unsigned"
  nSamples=nrow(datExpr)
  power = ifelse(nSamples<20, ifelse(type == "unsigned", 9, 18),
                 ifelse(nSamples<30, ifelse(type == "unsigned", 8, 16),
                        ifelse(nSamples<40, ifelse(type == "unsigned", 7, 14),
                               ifelse(type == "unsigned", 6, 12))  
                 )
  )
}
# 将选择好的最佳power值和soft-thresholding的计算结果保存至RData文件中
save(sft, power, file = "step2_power_value.Rdata")

######一步法构建加权共表达网络，识别基因模块######
if(T){
  cor <- WGCNA::cor
  # blockwiseModules函数构建加权共表达网络，并识别基因模块
  net <- blockwiseModules(
    datExpr,  # datExpr是包含基因表达数据的矩阵
    power = power,  # power是上一步计算得到的软阈值幂次
    maxBlockSize = ncol(datExpr),  # 计算机能处理的最大模块的基因数量
    corType = "pearson", #计算相关性的方法默认为"pearson","bicor"则更能考虑离群点的影响
    networkType = "unsigned",  # 无向网络 
    TOMType = "unsigned",   # 建立邻接矩阵和TOM矩阵时，是否考虑正负相关性，unsigned表示只考虑正相关性
    # minModuleSize指每个模块中基因的最少数目
    minModuleSize = 10,    ##越大模块越少
    # mergeCutHeight指合并模块的阈值
    mergeCutHeight = 0.25, ##越大模块越少
    numericLabels = TRUE,  # 返回数字作为模块的名称，后面可以再转换为颜色
    saveTOMs = F,  # 是否存储TOM矩阵，TOM矩阵计算最耗费时间的步骤之一，如果需要多次分析，可以将其存储起来供后续使用
    verbose = 3  # 控制输出信息的详细程度，数值越大输出的信息越多
  )
  # 统计每个模块中基因的数量
  table(net$colors) 
  # 0  1  2  3  4 
  # 18 56 22 16 10
  # power: 上一步计算的软阈值
  # maxBlockSize:计算机能处理的最大模块的基因数量(默认5000),16G内存可以处理2万个，
  # 计算资源允许的情况下最好放在一个block里面。
  # corType：计算相关性的方法；可选pearson(默认)，bicor。后者更能考虑离群点的影响。
  # networkType:计算邻接矩阵时，是否考虑正负相关性；默认为"unsigned",可选"signed", "signed hybrid"
  # TOMType：计算TOM矩阵时，是否考虑正负相关性；默认为"signed",可选"unsigned"。但是根据幂律转换的邻接矩阵(权重)的非负性，所以认为这里选择"signed"也没有太多的意义。
  # numericLabels: 返回数字而不是颜色作为模块的名字，后面可以再转换为颜色
  # saveTOMs：最耗费时间的计算，可存储起来供后续使用，
  # mergeCutHeight: 合并模块的阈值，越大模块越少,一般为0.25
  # minModuleSize: 每个模块里最少放多少个基因，设定越大模块越少
  # 输出结果根据模块中基因数目的多少，降序排列，依次编号为 `1-最大模块数`。
  # **0 (grey)**表示**未**分入任何模块的基因。
}

## 模块可视化，层级聚类树展示各个模块
# 模块构建完成后，将模块颜色labels转换为可视化所需的颜色，最后用层级聚类树展示各个模块
if(T){
  # Convert labels to colors for plotting
  moduleColors <- labels2colors(net$colors)
  table(moduleColors)
  # Plot the dendrogram and the module colors underneath
  pdf("step3_genes-modules_ClusterDendrogram.pdf",width = 16,height = 12)
  plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
                      "Module colors",
                      dendroLabels = FALSE, hang = 0.03,
                      addGuide = TRUE, guideHang = 0.05)
  dev.off()
}
# 最后将得到的模块信息和模块颜色保存起来，方便后续的分析和可视化使用
save(net, moduleColors, file = "step3_genes_modules.Rdata")

关联基因模块与表型

# 分组
datTraits <- data.frame(
  row.names = rownames(datExpr),
  exp_group = group_list,
  exp_group_No = rep(c(1,2),each=3)
)
datTraits

if(T){
  datTraits$exp_group <- as.factor(datTraits$exp_group)
  # 表达矩阵对应的表型变量为设计矩阵，并设置组别
  design <- model.matrix(~0+datTraits$exp_group)
  colnames(design) <- levels(datTraits$exp_group) #获取组别信息
  # 计算所有模块的特征基因eigengenes.
  MES0 <- moduleEigengenes(datExpr,moduleColors)$eigengenes
  # MEblue      MEbrown       MEgrey   MEturquoise     MEyellow
  # LH_41_hg19 -0.054224196 -0.113190330  0.017756936 -0.0596382729  0.202767908
  # LH_42_hg19 -0.088915385 -0.133718403 -0.078586736 -0.0917636118 -0.026137421
  MEs <- orderMEs(MES0)  #将相关的特征基因放在一起
  # > dim(MEs)
  # [1] 40  5
  # > MEs    # 其实就换了下module的顺序
  # MEblue     MEyellow      MEbrown   MEturquoise       MEgrey
  # LH_41_hg19 -0.054224196  0.202767908 -0.113190330 -0.0596382729  0.017756936
  # LH_42_hg19 -0.088915385 -0.026137421 -0.133718403 -0.0917636118 -0.078586736
  
  # 可操作地方！！！！！
  # 计算模块与表型之间的相关性矩阵
  moduleTraitCor <- cor(MEs,design,use = "p")
  # > moduleTraitCor
  # [,1]       [,2]
  # MEblue       0.6998699 -0.6998699
  # MEyellow    -0.4405819  0.4405819
  # MEbrown      0.4076938 -0.4076938
  # MEturquoise  0.5932559 -0.5932559
  # MEgrey       0.1919625 -0.1919625
  moduleTraitPvalue <- corPvalueStudent(moduleTraitCor,nSamples)
  # [,1]         [,2]
  # MEblue      5.014994e-07 5.014994e-07
  # MEyellow    4.437355e-03 4.437355e-03
  # MEbrown     9.020902e-03 9.020902e-03
  # MEturquoise 5.474019e-05 5.474019e-05
  # MEgrey      2.353641e-01 2.353641e-01
  
  # 设定热图的文本标签为矩阵中每个元素的相关性值及其对应的显著性水平
  textMatrix <- paste0(signif(moduleTraitCor,2),"\n(",
                       signif(moduleTraitPvalue,1),")")
  # > textMatrix
  # [1] "0.7\n(5e-07)"   "-0.44\n(0.004)" "0.41\n(0.009)"  "0.59\n(5e-05)" 
  # [5] "0.19\n(0.2)"    "-0.7\n(5e-07)"  "0.44\n(0.004)"  "-0.41\n(0.009)"
  # [9] "-0.59\n(5e-05)" "-0.19\n(0.2)" 
  # signif 保留位数
  # dim函数用于更改或显示矩阵的维数，包括行数和列数。
  # 在这里，使用dim(textMatrix)函数来设置textMatrix的维度，
  # 其行数和列数与模块与表型相关性的矩阵相同，以确保文本正确地添加到每个网格中。
  dim(textMatrix) <- dim(moduleTraitCor)
  # [,1]             [,2]            
  # [1,] "0.7\n(5e-07)"   "-0.7\n(5e-07)" 
  # [2,] "-0.44\n(0.004)" "0.44\n(0.004)" 
  # [3,] "0.41\n(0.009)"  "-0.41\n(0.009)"
  # [4,] "0.59\n(5e-05)"  "-0.59\n(5e-05)"
  # [5,] "0.19\n(0.2)"    "-0.19\n(0.2)" 
  
  pdf("step4_Module-trait-relationship_heatmap.pdf",
      width = 2*length(colnames(design)), 
      height = 0.6*length(names(MEs)) )
  par(mar=c(5, 9, 3, 3)) #留白：下、左、上、右
  # 绘制
  labeledHeatmap(Matrix = moduleTraitCor,  # 要绘制热图的矩阵
                 xLabels = colnames(design),  # X轴上的标签
                 yLabels = names(MEs),  # Y轴上的标签
                 ySymbols = names(MEs),  # 线条的标记符号
                 colorLabels = F,  # 颜色标签 如果为TRUE，则在热图下方添加一个注释区，表示颜色对应的值的范围
                 colors = blueWhiteRed(50),  # 颜色网格的颜色图谱
                 textMatrix = textMatrix,  # 添加到每个网格中的文本（可以是数字或字符）
                 setStdMargins = F,  # 是否将margins调整为标准大小
                 cex.text = 0.5,  # 文本大小
                 zlim = c(-1,1),   # 颜色图谱的值范围
                 main = "Module-trait relationships")  # 热图的标题
  dev.off()
  save(design, file = "step4_design.Rdata")
}

取MAD前2W转录本进行：

可以发现MRgreen模块与表型关联很强

可以进一步找hub，这里就直接看看作者的网络中的34个节点在不在这MEgreen 958个转录本对应的gene symbol中

id转换:

lncRNA的一些基础知识

LNCipedia进不去

询问曾老师：这是软件给的编号

文章作者也没给出对应的ENSEMBL id

看看mRNA算了

> aut_hub_mrna <- c("SWT1","AXIN2","HARS2","MMS19","HYLS1","PTRH1","MXRA7","PUF60","ASAH1","CLEC16A","NFIX","ZNF638","FIGNL1","BID","CNOT2","PTPN12","CD44","NCAM1","CANT1","DUSP14","REPS1","MFAP3L","EPB41","TCF4","IFNGR1","RNASEH2B")
> table(aut_hub_mrna %in% annos$hgnc_symbol)

FALSE  TRUE 
   21     5 
> aut_hub_mrna[aut_hub_mrna %in% annos$hgnc_symbol]
[1] "NFIX"     "EPB41"    "TCF4"     "IFNGR1"   "RNASEH2B"

虽然我们只取MAD前2W转录本进行WGCNA，发现作者的mRNA节点中只有5个在我们WGCNA结果与表型最显著的模块中，但是可以清晰地发现这五个基因在作者共表达网络中的关系

顺便看看PPI 的hub gene：

> PPI_hub <- c("MYC","EGF","AKT1","SUPT20H","UBA52")
> PPI_hub[PPI_hub %in% annos$hgnc_symbol]
[1] "SUPT20H" "UBA52"

所以无论是共表达网络、PPI还是WGCNA都是我们寻找hub gene进行narrow down的手段

以上就是初探mRNA、lncRNA联合分析下游的全部内容，我们部分复现了作者的流程并联系一些参考资料和其它手段进行了学习、解读，但这篇文章只对目前已知lncRNA定量，不涉及鉴定的新lncRNA的流程，毕竟这是人的数据，下周我们将一起看看lncRNA鉴定的全套流程

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2023-08-08，如有侵权请联系 cloudcommunity@tencent.com 删除

函数

本文分享自生信菜鸟团微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！