献上诚意满满的WGCNA分析

大家好，我是阿琛。今天我们来看看如何进行WGCNA分析。WGCNA分析，又称为加权基因共表达网络分析。该分析方法旨在寻找协同表达的基因模块(module)，并探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。
WGCNA分析适用于复杂的数据模式，根据官方推荐，一般需要5组或者15个样品以上的数据；同时，其适用的研究方向包括：不同器官或组织类型发育调控、同一组织不同发育调控、病原菌侵染后不同时间点应答等等。
下面，我们一起来看一下如何进行WGCNA分析。
R包安装与数据读取

1.1 R包安装

#if (!requireNamespace("BiocManager", quietly = TRUE))# install.packages("BiocManager")#BiocManager::install("preprocessCore")#install.packages("WGCNA")library("WGCNA") #引用WGCNA包
关于WGCNA的分析方法，官方存在专门的R包，即WGCNA包；该包储存在CRAN上，直接使用install.packages()命令即可完成安装。此外，由于在分析过程中涉及多线程分析，需要安装R的preprocessCore包。

1.2 数据读取

接下来，是输入数据的准备工作。对于WGCNA分析，测序数据和芯片数据均适用。如果是转录组数据，最好是RPKM，TPM或者其他归一化后的表达量；而如果是芯片数据的话，那么使用常规的归一化矩阵即可。在此，我们使用转录组测序的TPM数据来进行展示。
rt=read.table("TCGA.exp.txt",sep="\t",header=T,check.names=F)rt=as.matrix(rt)rownames(rt)=rt[,1]exp=rt[,2:ncol(rt)]dimnames=list(rownames(exp),colnames(exp))exp=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
接下来，我们可以查看一下具体的数据结构。
dim(exp)length(colnames(exp)
结果显示，该数据集由407例样本，8840个基因组成。
检查样本和基因
准备好输入数据以后，我们需要进一步对输入数据进行检查。

2.1 数据格式整理

datExpr0 = as.data.frame(t(exp))####截取部分数据查看datExpr0[1:3,1:3]
由于WGCNA分析是基于基因表达的聚类分析，所以在分析之前，我们需要对其进行转置。

2.2 检查缺失值

此外，我们需要判断一下基因和样品质量的好坏。使用goodSamplesGenes()函数，可以检查数据中缺失的条目，权重低于阈值的条目以及零方差基因，并返回具有最大缺失值或低权重值标准的样本和基因列表。
gsg = goodSamplesGenes(datExpr0, verbose = 3)gsg$allOK
可以看到，经过判断，结果显示TRUE，表明其中不包含有低质量的样品或基因名。
接下来，我们需要进一步判断和筛选质量不行的样品和基因。
if (!gsg$allOK){if (sum(!gsg$goodGenes)>0) printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")));if (sum(!gsg$goodSamples)>0) printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]}
当然，如果gsg$allOK的结果为FALSE，经过取反以后，!gsg$allOK的结果为TRUE，因此需要进一步分别筛选基因和样品质量。
最后，得到407个样品，8840个基因的表达矩阵。
样品聚类

3.1 聚类

随后，我们使用dist()函数计算样品之间的欧式距离，并通过聚类查看是否存在异常样品。
sampleTree = hclust(dist(datExpr0), method = "average")

3.2 画图

运算完毕后，将得到的聚类结果进行可视化。
par(cex = 0.6)par(mar = c(0,4,2,0))plot(sampleTree)plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2)
结果显示：可以看到，最左边明显有几个样品和别人不在一条枝叶上；因此，我们需要将那几个明显偏离的分支给剪去。
abline(h = 155, col = "red")
根据聚类树的结果，我们需要将异常样品进行去除；在此，将剪切线的地方定于高度为155的位置。
clust = cutreeStatic(sampleTree, cutHeight = 155, minSize = 10)table(clust)keepSamples = (clust==1)datExpr0 = datExpr0[keepSamples, ]
结果显示，以剪切线所在位置作为分界，异常样品为1例；提取适合样品的表达信息，最终得到了406例样品纳入随后的分析。
临床数据的准备
接下来，我们需要根据样品的分组信息，制作临床信息表格。当然，除了自己制作外，也可以从工作目录中提取准备好的表格。
traitData=data.frame(Normal=c(rep(1,32),rep(0,375)), Tumor=c(rep(0,32),rep(1,375)))dim(traitData)row.names(traitData)=colnames(tpms)
根据正常和肿瘤样品信息，得到了新的数据框储存相应的临床信息。
sameSample=intersect(rownames(datExpr0), rownames(traitData))datExpr0=datExpr0[sameSample,]datTraits=traitData[sameSample,]
由于前面以及删除了1例异常样品，我们通过intersect()函数获取共同的样品名，接着通过取子集的方式，最终得到了406例样品的表达和临床矩阵。
重新聚类
随后，我们对新得到的样品重新聚类。
sampleTree2 = hclust(dist(datExpr0), method = "average")plot(sampleTree2)
结果显示：
traitColors = numbers2colors(datTraits, signed = FALSE)##画图plotDendroAndColors(sampleTree2, traitColors, groupLabels = names(datTraits), main = "Sample dendrogram and trait heatmap")save(datExpr0, datTraits, file = "01-dataInput.RData")
根据临床分组信息，使用number2colors()函数为临床信息创建颜色集，其中红色表示阳性，白色表示阴性；接着，使用plotDendroAndColors()函数，在树状图中绘制了层次聚类树状图以及基于颜色注释的临床分组特征。
构建网络，识别模块
1、选择合适的软阈值（soft thresholding power）即beta值
rm(list = ls())load(file = "01-dataInput.RData")enableWGCNAThreads() #多线程工作powers = c(1:20) #幂指数范围1:20sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)
选取1到20为候选的power值，使用pickSoftThreshold()函数对相应的power值进行计算；
par(mfrow = c(1,2))cex1 = 0.9###拟合指数与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");abline(h=0.90,col="red") #可以修改###平均连通性与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")
根据结果，绘制power值对应的散点图；结果显示，1到2，2到3，整个变化趋势都还很明显，而3往后，随着power值的增大，其变化趋势已经不明显了；因此，选取3或者4为最佳的power值。
2、选取最佳的power值，建立邻近矩阵，根据连接度使我们的基因分布符合无尺度网络
sft #查看最佳power值softPower =sft$powerEstimate #最佳power值adjacency = adjacency(datExpr0, power = softPower)
3、计算拓扑矩阵（TOM矩阵）
TOM = TOMsimilarity(adjacency);dissTOM = 1-TOM
4、基因聚类
geneTree = hclust(as.dist(dissTOM), method = "average")plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity", labels = FALSE, hang = 0.04)
使用基因表达得到的TOM矩阵，再次对基因继续聚类；
5、动态剪切模块识别
设定最小模块基因数量，对基因聚类结果进行剪切，得到不同的基因模块。
minModuleSize = 50 #模块基因数目dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize);table(dynamicMods)
使用cutreeDynamic()函数，最终得到了13个不同的基因模块；
dynamicColors = labels2colors(dynamicMods)table(dynamicColors)
接着，我们通过labels2colors()函数，对每一个基因模块进行颜色赋值；
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors")
最后，基于plotDendroAndColors()函数，将剪切结果添加到聚类树上，并使用前面赋予的模块颜色来进行区分。
相似模块聚类
得到不同的基因模块后，我们需要对表达模式相似的模块进行聚类。
MEList = moduleEigengenes(datExpr0, colors = dynamicColors)MEs = MEList$eigengenesMEDiss = 1-cor(MEs);METree = hclust(as.dist(MEDiss), method = "average")plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")
结果显示：
当然，如果对聚类结果还需要手工调整模块的话，我们可以根据Height，设定MEDissTheres的阈值，对模块进行人工的合并。根据这个结果，不同模块之间的距离相对可以，就不进行合并了。
MEDissThres = 0.25 #剪切高度可修改abline(h=MEDissThres, col = "red")
结果显示：
最后，根据人工设定的剪切高度，对相似的基因模块进行合并。
merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)mergedColors = merge$colorsmergedMEs = merge$newMEsplotDendroAndColors(geneTree, mergedColors,"Dynamic Tree Cut", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors")
结果显示：
moduleColors = mergedColorstable(moduleColors)
可以看到，此时的模块和前面是保持一致的。
colorOrder = c("grey", standardColors(50))moduleLabels = match(moduleColors, colorOrder)-1MEs = mergedMEs
模块与性状数据热图
随后，我们需要计算不同的基因模块与临床性状之间的相关性系数（cor）和对应的P值（Pvalue）。
nGenes = ncol(datExpr0)nSamples = nrow(datExpr0)moduleTraitCor = cor(MEs, datTraits, use = "p")moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
然后，根据相关系数和p值，将其使用热图进行展示。
##绘图展示textMatrix = paste(signif(moduleTraitCor, 2), "\n(",signif(moduleTraitPvalue,1), ")", sep = "")dim(textMatrix) = dim(moduleTraitCor)par(mar = c(5, 10, 3, 3))labeledHeatmap(Matrix = moduleTraitCor,xLabels = names(datTraits),yLabels = names(MEs),ySymbols = names(MEs),colorLabels = FALSE,colors = blueWhiteRed(50),textMatrix = textMatrix,setStdMargins = FALSE,cex.text = 0.5,zlim = c(-1,1),main = paste("Module-trait relationships"))
结果显示了每一个基因模块与正常或肿瘤样品之间的相关性结果。我们可以根据相关系数和p值，挑选与临床性状之间最显著相关的模块进行后续分析。
不同模块与临床形状的具体分析
尽管计算得到了临床性状与基因模块之间的相关性，可以挑选最相关的那些模块来分析，但是模块本身仍然包含非常多的基因，还需进一步的寻找最重要的基因。所有的模块都可以跟基因算出相关系数，所有的连续型性状也可以跟基因的表达值算出相关系数。
1、计算模块与基因之间的相关性矩阵
modNames = substring(names(MEs), 3)geneModuleMembership = as.data.frame(cor(datExpr0, MEs, use = "p"))MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))names(geneModuleMembership) = paste("MM", modNames, sep="")names(MMPvalue) = paste("p.MM", modNames, sep="")
2、计算临床性状与基因之间的相关性矩阵
traitNames=names(datTraits)geneTraitSignificance = as.data.frame(cor(datExpr0, datTraits, use = "p"))GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))names(geneTraitSignificance) = paste("GS.", traitNames, sep="")names(GSPvalue) = paste("p.GS.", traitNames, sep="")
3、批量输出性状与模块散点图：把两个相关性矩阵联合起来,指定感兴趣模块进行分析
for (traitintraitNames){ traitColumn=match(trait,traitNames) for (module in modNames){ column = match(module, modNames) moduleGenes = moduleColors==moduleif (nrow(geneModuleMembership[moduleGenes,]) > 1){ outPdf=paste(trait, "_", module,".pdf",sep="") pdf(file=outPdf,width=7,height=7) par(mfrow = c(1,1)) verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]), abs(geneTraitSignificance[moduleGenes, traitColumn]), xlab = paste("Module Membership in", module, "module"), ylab = paste("Gene significance for ",trait),main = paste("Module membership vs. gene significance\n"), cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module) abline(v=0.8,h=0.5,col="red") dev.off() } }}
输出每个模块的基因
最后，我们需要把每个模块里面的基因进行输出，以用于后续的分析。
for(mod in 1:nrow(table(moduleColors))){modules = names(table(moduleColors))[mod]probes = colnames(datExpr0)inModule = (moduleColors == modules)modGenes = probes[inModule]write.table(modGenes,file =paste0(modules,".txt"),sep="\t",row.names=F,col.names=F,quote=F)}
逐个输出每个颜色模块的基因，并且根据与临床性状之间的相关性，选取最相关的基因模块，对这个基因模块进行后续的一系列分析。到此，WGCNA分析的主线部分基本就完成了。当然，WGCNA分析不只是如此，还可以指定基因模块，绘制表达与性状之间的热图等等。
—END—
继续阅读
阅读原文