lncRNA实战项目-第六步-WGCNA相关性分析 / 开普饭

响应生信技能树的号召：lncRNA数据分析传送门 , 一起来一个lncRNA数据分析实战！

往期目录

lncRNA实战项目-第二步-了解文章及数据

lncRNA实战项目-第三步-了解参考基因组及注释文件

lncRNA实战项目-第四步-得到表达矩阵的流程

lncRNA实战项目-第五步-差异表达的mRNA和lncRNA

因为样本数量比较可观，所以可以进行WGCNA分析。这里是并不需要选取所有的基因来做WGCNA分析，挑选的标准可以是top变异程度大的基因集合，或者显著差异表达的基因集合等等。

这里可以参考：https://github.com/jmzeng1314/my_WGCN

WGCNA将lncRNA分成18个模块（3635个lncRNA），空间模块中lncRNA表达呈现明显的组织区域特异性，如：CB (M1, 794个lncRNAs),DG/CA1 (M2, 443个lncRNAs), CA1 (M4, 369个lncRNAs),neocortex (M7, 123个lncRNAs)和OC (M10,57个lncRNAs)。时间模块中lncRNA表达与年龄有关，而与组织区域不明显；性别模块中lncRNA表达与性别和年龄都相关。每个模块就必须做pathway/go等数据库的注释分析咯！

资料收集：

google搜索或在生信技能树和生信菜鸟团搜索WGCNA ，能找到很多教程，下面列出几个中文教程和英文教程，强烈推荐中文教程1和英文教程3。

一文学会WGCNA分析
学习WGCNA总结
Tutorials for the WGCNA package

WGCNA Background and glossary
Data input and cleaning
Network construction and module detection
Relating modules to external information and identifying important genes
Interfacing network analysis with other data such as functional annotation and gene ontology
Network visualization using WGCNA functions
Exporting a gene network to external visualization software

背景知识

基本概念:

WGCNA(weighted correlation network analysis)加权基因共表达网络分析, 用于提取与性状或临床特征相关的基因模块，解析基础代谢途径，转录调控途径、翻译水平调控等生物学过程。WGCNA适合于复杂的数据模式，推荐5组以上的数据，如：
不同器官、组织类型发育调控；
同一组织不同时期发育调控；
非生物胁迫不同时间点应答；
病原物侵染后不同时间点应答。

基本步骤：

WGCNA分为表达量聚类分析和表型关联两部分，具体步骤包括基因之间相关系数的计算，共表达网络的构建，筛选特定模块，模块与性状关联，核心基因的筛选。

术语：

Co-expression weighted network: 是一个无向有权重（undirected, weighted）的网络。“无权重（unweighted network）”,基因与基因之间的相关度只能是0或者1,0表示两个基因没有联系，而1表示有。“有权重（weighted network）”基因之间不仅仅是相关与否，还记录着它们的相关性数值，数值就是基因之间的联系的权重(相关性)。
Module:（模块）指表达模式相似的基因聚为一类，这样的一类基因称为模块。
Connectivity（连接度）:连接度是指某个基因与该网络中其他基因的相关程度，常为“相关度”之和。
Eigengene（eigen- +‎ gene）：基因和样本构成的矩阵。
Hub gene: 模块中的高度相关的核心基因。
Gene signicance ，GS:基因显著性参数，为非负数字信息，比如基于样本的临床数据（clinical trails）和基于每个基因的-log(p-value)等。

分析流程

WGCNA输入文件需要一个表达矩阵，最好是RPKM或其他归一化好的表达量，还需要一个矩阵关于临床信息或者其它关于样本属性的信息。

STEP1: 输入数据的准备

表达矩阵可以从作者GitHub上下载（https://github.com/DChenABLife/RhesusLncRNA），我这里只尝试对lncRNA的表达矩阵(https://github.com/DChenABLife/RhesusLncRNA/blob/master/All_sample_LncRNA_exp_RPKM.xlsx)做后续分析。下载的表达矩阵文件是Excel格式的，需要转为csv格式方便后续用R处理，可以直接打开这个excel文件，然后另存为csv格式即可。

读入原始表达数据

原始数据包含64个样本，9904个lncRNA表达量，这时的矩阵行为基因，列为样本信息，其中第一列是lncRNA ID，(这里的lncRNA ID是cufflinks 组装时给的自由的ID, 是需要和已有的ID对应, 对于新的转录本再通过nr/nt等数据库注释），第66列是作者给出的注释（我查了几个注释有的也查不出来是什么意思）。

setwd("G:/My_exercise/WGCNA/") lncRNAexpr <- read.csv("All_sample_LncRNA_exp_RPKM.csv",sep=",",header = T) head(lncRNAexpr) dim(lncRNAexpr) #[1] 9904 66

重命名数据列表，行名和列名

##去掉Annotation这列 fpkm <- lncRNAexpr[,-66] head(fpkm) ##重命名行名和列名 rownames(fpkm)=fpkm[,1] fpkm=fpkm[,-1] fpkm[1:4,1:4]

准备表型信息

这里有64个样本，包含猕猴脑不同空间区域，不同发育时期，以及性别，因为每个样本都交叉包含着三种不同的信息，如果选择全部表型信息，我试了试，后续的模块和性状完全看不清关系，所以我这里仅选择脑不同区域的表型信息,包括CB、DG、PFC、PCC、CA1、OC、PC、TC。

##Sample Info subname=sapply(colnames(fpkm),function(x) strsplit(x,"_")[[1]][1]) datTraits = data.frame(gsm=names(fpkm), subtype=subname) rownames(datTraits)=datTraits[,1] head(datTraits)

下载并载入WGCNA包

source("http://bioconductor.org/biocLite.R") #biocLite(c("AnnotationDbi", "impute", "GO.db", "preprocessCore")) ##如果上面的包已经下载过了，就不用下载 biocLite("WGCNA") library(WGCNA)

行列转置

WGCNA针对的是基因进行聚类，而一般我们的聚类是针对样本用hclust即可，也就是说要转置为行表示样本, 列表示基因。

RNAseq_voom <- fpkm WGCNA_matrix = t(RNAseq_voom[order(apply(RNAseq_voom,1,mad), decreasing = T)[1:5000],]) datExpr <- WGCNA_matrix ## top 5000 mad genes datExpr[1:4,1:4]

确定临床表型与样本名字

sampleNames = rownames(datExpr); traitRows = match(sampleNames, datTraits$gsm) rownames(datTraits) = datTraits[traitRows, 1]

datExpr和datTraits准备好后，接下来就是构建基因网络，鉴定模块。网络构建有三种方法：1）一步法构建网络；2）多步法构建网络；3）block-wise构建网络（主要针对大数据集）。下面的介绍的步骤是一步法构建网络。

STEP2：确定最佳soft-thresholding power

选择合适的“软阀值（soft thresholding power）”beta, 用到的函数是pickSoftThreshold，pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)，powerVector可以是一系列数值，从而选择最优值。这个函数返回一个列表，第一项是powerEstimate是估计的最优power；第二项是fitIndices是详细的矩阵数据，其中第五列是mean.k表示平均“连接度（connectivity）”。

Constructing a weighted gene network entails the choice of the soft thresholding power to which co-expression similarity is raised to calculate adjacency.

# Choose a set of soft-thresholding powers powers = c(c(1:10), seq(from = 12, to=20, by=2)) # Call the network topology analysis function sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5) # Plot the results: par(mfrow = c(1,2)); cex1 = 0.9; # Scale-free topology fit index as a function of the soft-thresholding power 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", main = paste("Scale independence")); 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 abline(h=0.90,col="green") # Mean connectivity as a function of the soft-thresholding power plot(sft$fitIndices[,1], sft$fitIndices[,5], xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n", main = paste("Mean connectivity")) text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red") # 选择beta值 best_beta=sft$powerEstimate #> best_beta [1] 3

最佳beta值是3。

STEP3: 一步法构建共表达矩阵

一步法构建网络，使用函数blockwiseModules(), 这个函数包含很多参数，其中power=sft$powerEstimate=3 即上一步得到的最佳软阈值；maxBlockSize 默认为5000，表示在这个数值内的基因将整体被计算，如果调大需要更多的内存；numerricLabels 默认为返回颜色，设置为TRUE则返回数字；mergeCutHeight是合并模块阈值的一个参数。

net = blockwiseModules(datExpr, power = sft$powerEstimate, maxBlockSize = 6000,TOMType = "unsigned", minModuleSize = 30,reassignThreshold = 0, mergeCutHeight = 0.25, numericLabels = TRUE, pamRespectsDendro = FALSE, saveTOMs = TRUE, saveTOMFileBase = "AS-green-FPKM-TOM", verbose = 3)

STEP4:模块鉴定及可视化

模块鉴定

上一步的返回结果是一个列表，可以用table()函数查看，0表示没有任何module接受。table(net$colors) 可以看总共有多少模块，每个模块的大小，这里共有9个模块，从1-9每个模块的大小是递减的，从2254-115，0表示这些基因不在所有模块内。

table(net$colors)

可视化

dendrograms表示在一个block中所有基因的进化树图, 使用函数plotDendroAndColors()查看系统发生树；blockGenes是一个block中所有的基因。

# Convert labels to colors for plotting mergedColors = labels2colors(net$colors) table(mergedColors) # Plot the dendrogram and the module colors underneath plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]], "Module colors", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05) ## assign all of the gene to their corresponding module ## hclust for the genes.

#明确样本数和基因数 nGenes = ncol(datExpr) nSamples = nrow(datExpr) #首先针对样本做个系统聚类树 datExpr_tree<-hclust(dist(datExpr), method = "average") par(mar = c(0,5,2,0)) plot(datExpr_tree, main = "Sample clustering", sub="", xlab="", cex.lab = 2, cex.axis = 1, cex.main = 1,cex.lab=1) ## 如果这个时候样本是有性状，或者临床表型的，可以加进去看看是否聚类合理 #针对前面构造的样品矩阵添加对应颜色 sample_colors <- numbers2colors(as.numeric(factor(datTraits$subtype)), colors = c("grey","blue","red","green"),signed = FALSE) ## 这个给样品添加对应颜色的代码需要自行修改以适应自己的数据分析项目。 # sample_colors <- numbers2colors( datTraits ,signed = FALSE) ## 如果样品有多种分类情况，而且 datTraits 里面都是分类信息，那么可以直接用上面代码， ##当然，这样给的颜色不明显，意义不大。 #构造10个样品的系统聚类树及性状热图 par(mar = c(1,4,3,1),cex=0.8) plotDendroAndColors(datExpr_tree, sample_colors, groupLabels = colnames(sample), cex.dendroLabels = 0.8, marAll = c(1, 4, 3, 1), cex.rowText = 0.01, main = "Sample dendrogram and trait heatmap")

STEP5:模块和性状的关系

MEs是一个关于modules的特征量矩阵，行数等于筛选的modules数，列数等于样本数；

design=model.matrix(~0+ datTraits$subtype) colnames(design)=levels(datTraits$subtype) moduleColors <- labels2colors(net$colors) # Recalculate MEs with color labels MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes MEs = orderMEs(MEs0); ##不同颜色的模块的ME值矩阵(样本vs模块) moduleTraitCor = cor(MEs, design , use = "p"); moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

sizeGrWindow(10,6) # Will display correlations and their p-values textMatrix = paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")", sep = ""); dim(textMatrix) = dim(moduleTraitCor) par(mar = c(6, 8.5, 3, 3)); # Display the correlation values within a heatmap plot labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(design), yLabels = names(MEs), ySymbols = names(MEs), colorLabels = FALSE, colors = greenWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.5, zlim = c(-1,1), main = paste("Module-trait relationships"))

图中第二列第五行，即CB/turquoise相关性最强（0.97），pvalue=1e-41,后续分析可以挑选这个模块。
每一列对应的样本特征可以从design这里查看。

STEP6:感兴趣性状的模块的具体基因分析

选定CB/turquoise这个模块后，接下来看看模块与样本特征和基因的相关性，1. 首先计算全模块与基因的相关性矩阵，2. 再计算性状与基因的相关性矩阵，3. 使用verboseScatterplot()绘制相关性图：

首先计算模块与基因的相关性矩阵

# names (colors) of the modules modNames = substring(names(MEs), 3) geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p")); ## 算出每个模块跟基因的皮尔森相关系数矩阵 ## MEs是每个模块在每个样本里面的值 ## datExpr是每个基因在每个样本的表达量 MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples)); names(geneModuleMembership) = paste("MM", modNames, sep=""); names(MMPvalue) = paste("p.MM", modNames, sep="");

再计算性状与基因的相关性矩阵

## 只有连续型性状才能只有计算 ## 这里把是否属于 CB 表型这个变量用0,1进行数值化。 CB = as.data.frame(design[,2]); names(CB) = "CB" geneTraitSignificance = as.data.frame(cor(datExpr, CB, use = "p")); GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples)); names(geneTraitSignificance) = paste("GS.", names(CB), sep=""); names(GSPvalue) = paste("p.GS.", names(CB), sep="")

最后把两个相关性矩阵联合起来,指定感兴趣模块进行分析

module = "turquoise" column = match(module, modNames); moduleGenes = moduleColors==module; sizeGrWindow(7, 7); par(mfrow = c(1,1)); verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]), abs(geneTraitSignificance[moduleGenes, 1]), xlab = paste("Module Membership in", module, "module"), ylab = "Gene significance for CB", main = paste("Module membership vs. gene significance\n"), cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

上图可以看出基因跟其对应的性状高度相关，可以导出做个GO/KEGG注释，看看这些基因的具体功能。

STEP7:网络的可视化

Topological Overlap Matrix (TOM)图是指基因和筛选模型的可视化表示。WGCNA认为基因之间的简单的相关性不足以计算共表达，所以它利用邻近矩阵，又计算了一个新的邻近矩阵。一般来说，TOM就是WGCNA分析的最终结果，后续的只是对TOM的下游注释。

生成TOM图的步骤：1. 生成全基因不相似TOM矩阵，比如dissTOM=1-TOMsimilarityFromExpr(datExpr, power = 6)，可以把得到的不相关矩阵加幂，这样绘制的TOM图色彩差异会比较明显。2. 使用TOMplot()绘制TOM图。

#首先针对所有基因画热图 nGenes = ncol(datExpr) nSamples = nrow(datExpr) geneTree = net$dendrograms[[1]]; #生成全基因不相似TOM矩阵 dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6); plotTOM = dissTOM^7; diag(plotTOM) = NA; #TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")

#然后随机选取部分基因作图 nSelect = 400 # For reproducibility, we set the random seed set.seed(10); select = sample(nGenes, size = nSelect); selectTOM = dissTOM[select, select]; # There’s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster. selectTree = hclust(as.dist(selectTOM), method = "average") selectColors = moduleColors[select]; # Open a graphical window sizeGrWindow(9,9) # Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing # the color palette; setting the diagonal to NA also improves the clarity of the plot plotDiss = selectTOM^7; diag(plotDiss) = NA; TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")

图上和图左是全基因系统发育树，不同颜色亮带表示不同的module，每一个亮点表示每个基因与其他基因的相关性强弱（越亮表示相关性越强）。

#最后画模块和性状的关系 # Recalculate module eigengenes MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes ## 只有连续型性状才能只有计算 ## 这里把是否属于 Luminal 表型这个变量用0,1进行数值化。 CB = as.data.frame(design[,2]); names(CB) = "CB" # Add the weight to existing module eigengenes MET = orderMEs(cbind(MEs, CB)) # Plot the relationships among the eigengenes and the trait sizeGrWindow(5,7.5); par(cex = 0.9) plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle= 90) # Plot the dendrogram sizeGrWindow(6,6); par(cex = 1.0) ## 模块的聚类图 plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0), plotHeatmaps = FALSE) # Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot) par(cex = 1.0) ## 性状与模块热图 plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2), plotDendrograms = FALSE, xLabelsAngle = 90)

STEP8:提取指定模块的基因名

提取基因信息,可以做GO/KEGG等分析，进而解释这些module的生物学意义。这里的lncRNA ID转换着有点麻烦，这一步先略过，之后再看看。

# Select module module = "turquoise"; # Select module probes probes = colnames(datExpr) ## 我们例子里面的probe就是基因名 inModule = (moduleColors==module); modProbes = probes[inModule];

STEP9: 模块的导出

主要模块里面的基因直接的相互作用关系信息可以导出到cytoscape,VisANT等网络可视化软件。

# Recalculate topological overlap TOM = TOMsimilarityFromExpr(datExpr, power = 6); # Select module module = "turquoise"; # Select module probes probes = colnames(datExpr) ## 我们例子里面的probe就是基因名 inModule = (moduleColors==module); modProbes = probes[inModule]; ## 也是提取指定模块的基因名 # Select the corresponding Topological Overlap modTOM = TOM[inModule, inModule]; dimnames(modTOM) = list(modProbes, modProbes)

首先是导出到VisANT

vis = exportNetworkToVisANT(modTOM, file = paste("VisANTInput-", module, ".txt", sep=""), weighted = TRUE, threshold = 0)

也可以输出模块里的部分基因：

nTop = 30; IMConn = softConnectivity(datExpr[, modProbes]); top = (rank(-IMConn) <= nTop) vis = exportNetworkToVisANT(modTOM[top, top], file = paste("VisANTInput-", module, "-top30.txt", sep=""), weighted = TRUE, threshold = 0, probeToGene = data.frame(annot$substanceBXH, annot$gene_symbol) )

然后是导出到cytoscape

cyt = exportNetworkToCytoscape( modTOM, edgeFile = paste("CytoscapeInput-edges-", paste(module, collapse="-"), ".txt", sep=""), nodeFile = paste("CytoscapeInput-nodes-", paste(module, collapse="-"), ".txt", sep=""), weighted = TRUE, threshold = 0.02, nodeNames = modProbes, nodeAttr = moduleColors[inModule] );

编辑：jimmy

后续分析，请大家持续关注~

lncRNA实战项目-第六步-WGCNA相关性分析