WGCNA：构建共表达网络的神器

今天小编给大家分享一个构建共表达网络的神器，WGCNA 。

这个软件从2008年发表至今（截止到2019-05-20）已经被引用3899次。可见这个包多么受欢迎。

在正式介绍这个软件之前，先让我们看看用这个软件做的相关分析实例。

比如预测ncRNA相关功能：

https://www.sciencedirect.com/science/article/pii/S0378111917300446

比如挖掘性状相关基因模块：

当然其他的功能也有很多，但是小编就不一一介绍了。

那么到底怎么使用WGCNA呢？今天小编就以例子实操，一步一步为大家进行介绍。

1.使用WGCNA，首先要准备数据。

数据1:归一化好的基因表达数据 geneExp.txt (基因表达数据首先过滤掉方差为0的基因，如果基因数目还是很多，可以进一步过滤掉方差比较小的基因，不建议直接使用差异表达基因。另外基因表达需要做好归一化。)

数据2:性状数据 trait_use.txt

加载包，并倒入归一化后的数据。

操作命令均在当前文件夹，文件夹中有以下文件:

(1) geneExp.txt (2) trait_use.txt

library(WGCNA)data<-read.table("geneExp.txt ",sep="t",header=T,row.names=1)datExpr0<-t(data)

输入文件横轴是基因名，纵轴是归一化后在各个样本中的表达量，示例如下图：

2.选择软阈值进行模块划分。

powers = c(c(1:10), seq(from = 12, to=20, by=2))sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)

然后最佳阈值：

sft$powerEstimate

3.构建共表达网络

Toplogical Overlap Measure (TOM)

power是软阀值；

maxBlockSize表示在这个数值内的基因将整体被计算，如果调大需要更多的内存；4G内存电脑可处理8000-10000个，16G内存电脑可以处理2万个，32G内存电脑可以处理3万个，计算资源允许的情况下最好放在一个block里面。

TOMType = “unsigned” 一把选用“unsigned”，这个是最常用的。具体解释见下面。

If “unsigned”, the standard TOM will be used (more generally, TOM function will receive the adjacency as input). If “signed”, TOM will keep track of the sign of the adjacency between neighbors.

minModuleSize设定每个modle的最少的基因数目，可以调大；

reassignThreshold：是一个p值阈值，表示这个基因只有pValue低于这个阈值才会被分配到模块中。默认值是0.05.

mergeCutHeight是在计算完所有modules后，将特征量高度相似的modules进行和并，其设定的参数值代表合并模块的阈值，越大模块越少；

numericLabels默认为FALSE返回颜色，设定为TRUE则返回数字；

pamRespectsDendro逻辑值，默认为TRUE，可以设定为FALSE。是和pamStage一起使用的，当pamStage为ture时，这个也是ture时表明模块检测的第二阶段将会执行，一般设置为False.

verbose，如果是0则执行过程中的具体细节就不输出了，数值越大，程序执行细节输出的越详细。

代码：

net = blockwiseModules(datExpr0,power = sft$powerEstimate, maxBlockSize = 6000,TOMType = "unsigned", minModuleSize = 30,reassignThreshold = 0, mergeCutHeight = 0.25,numericLabels = TRUE, pamRespectsDendro = FALSE,verbose = 3)

对网络可视化

mergedColors = labels2colors(net$colors)table(mergedColors)pdf("module_14d.pdf")plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],"Module colors",dendroLabels = FALSE, hang = 0.03,addGuide = TRUE, guideHang = 0.05)dev.off()

图片文件名： module_14d.pdf

图片如下图所示：

4.绘制模块聚类图

MEList = moduleEigengenes(datExpr0, colors = mergedColors) MEs = MEList$eigengenesMET = orderMEs(MEs)pdf("module_cluster_14d.pdf")plotEigengeneNetworks(MET, '', marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab =0.8, xLabelsAngle= 90)dev.off()

图片文件名：module_cluster_14d.pdf

图片如下图所示：

5.绘制gene聚类网络

heatmap中每一行每一列都代表一个gene，注意的是如果绘制全部gene的会占很大内存。和大量的时间。一般会随机选择一些基因进行可视化。如下可以选择400个gene进行可视化。

nSelect = 400set.seed(10);dissTOM = 1-TOMsimilarityFromExpr(datExpr0, power = sft$powerEstimate);#如果基因数目多这个命令大量消耗内存和时间nGenes = ncol(datExpr0)select = sample(nGenes, size = nSelect);selectTOM = dissTOM[select, select];selectTree = hclust(as.dist(selectTOM), method = 'average')selectColors = mergedColors[select];plotDiss = selectTOM^7;diag(plotDiss) = NA;pdf("TOMplot_14d.pdf")TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")dev.off()

图片文件名：TOMplot_14d.pdf

图片如下图所示：

6.将每个模块的gene保存下来

modNames = substring(names(MEs), 3)AllGenes =colnames(datExpr0)for (module in modNames){modGenes = (mergedColors==module)modLLIDs = AllGenes[modGenes]; fileName = paste("14d-", module, ".txt", sep="");write.table(as.data.frame(modLLIDs), file = fileName,row.names = FALSE, col.names = FALSE)}

文件夹gene_module_14d中分别是挖掘出来的模块文件，每个文件中存储的是这个模块的基因名称。如gene_module_14d文件夹下14d-bisque4.txt文件中，存储的是bisque4模块中的gene名称。如下图：

7.计算模块和性状的关系。

moduleColors <- labels2colors(net$colors)

计算模块最具代表性的值的表达量，第一主成分的值.

MEs0 = moduleEigengenes(datExpr0, moduleColors)$eigengenes

读入性状文件

dataTrait<-read.table("trait_use.txt",sep="t",header=T,row.names=1)

性状文件，横轴是样本，纵轴是性状。如下图所示:

计算模块和性状之间的相关性

moduleTraitCor = cor(MEs0, dataTrait, use = "p")

p是一种计算方法，就是数据计算时考虑配对信息。

moduleTraitPvalue = corPvalueStudent(moduleTraitCor, 24)

绘制相关性结果

pdf("module_trait_cor_14d.pdf")textMatrix = paste(signif(moduleTraitCor, 2), "n(",signif(moduleTraitPvalue, 1), ")", sep = "");dim(textMatrix) = dim(moduleTraitCor)par(mar = c(6, 8.5, 3, 3));labeledHeatmap(Matrix = moduleTraitCor,xLabels = names(dataTrait),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"))

图片文件名：module_trait_cor_14d.pdf

图片如下图所示: