WGCNA如何挖掘潜在的共表达基因

WGCNA如何挖掘潜在的共表达基因

欢迎关注”生信修炼手册”!

共表达基因指的是表达量具有协同变化趋势的基因集合，通常认为这些基因参与相同的生物学过程，比如参与同一个代谢通路，正是由于功能上的协同作用，导致表达量呈现出高度相关性。

在WGCNA中，对传统的相关系数进行乘方运算，用最终得到的值来表征基因间的相关性。在计算出这样的相关性统计量值之后，如何确定哪些基因是共表达的呢？

WGCNA的做法是聚类分析，聚类分析属于一种非监督的机器学习算法，通过聚类树，可以观察到哪些基因在聚类树中属于同一分支，属于同一分支的基因可以归为一类。实际操作中，考虑到基因数目较多等情况，肯定需要算法来自动化的进行分类，WGCNA采用的是​​dynamicTreeCut​​这个R包。

对于聚类算法而言，需要输入基因间的距离矩阵，首先就需要将基因间的邻接矩阵转换为距离矩阵，对相关系数进行乘方运算，可以计算出邻接矩阵，但是这个值本质上反映的是基因间的相似度，并不是距离。在计算距离矩阵时，WGCNA采用了​​TOM​​这种统计量，该统计量可以表征网络中节点的相似性，计算公式如下

借助​​TOM​​值，将基因间的相关系数转换为了距离，然后就可以用该距离矩阵进行聚类。上述的计算方法在WGCNA中都有对应的公式，代码如下

# 确定乘方运算中power的最佳取值powers <- c(c(1:10), seq(from = 12, to=20, by=2))sft <- pickSoftThreshold(datExpr,powerVector = powers,verbose = 5)softPower <- sft$powerEstimate# 计算邻接矩阵adjacency <- adjacency(datExpr, power = softPower)# 计算TOM相似度矩阵TOM <- TOMsimilarity(adjacency)# 计算距离矩阵dissTOM <- 1-TOM# 聚类geneTree <- hclust(as.dist(dissTOM), method = "average")

根据聚类结果和距离矩阵，就可以调用​​dynamicTreeCut​​的算法来识别modules, 代码如下

# 指定每个module中基因数目的最小值minModuleSize <- 30# 识别modulesdynamicMods <- cutreeDynamic(dendro = geneTree,distM = dissTOM,deepSplit = 2,pamRespectsDendro = FALSE,minClusterSize = minModuleSize)

通过​​table​​函数可以查看modules的结果，用法如下

> table(dynamicMods)dynamicMods 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 1988 614 316 311 257 235 225 212 158 153 121 106 102 100 94 91 78 76 65 5820 21 2258 48 34

可以看到，识别出22个modules, ​​0​​​代表那些没有归入任何modules的基因。通过​​plotDendroAndColors​​函数可视化聚类树对应的modules, 代码如下

dynamicColors = labels2colors(dynamicMods)plotDendroAndColors(geneTree,dynamicColors,"Dynamic Tree Cut",dendroLabels = FALSE, hang = 0.03,addGuide = TRUE,guideHang = 0.05,main = "Gene dendrogram and module colors")

生成的图片如下

整个图片分为两个部分，上方为基因的聚类树，下方为识别到的modules， 不同的modules对应不同的颜色，其中灰色对应那些没有归入任何modules的基因。

通过​​dynamicTreeCut​​​识别到modules之后，还会结合每个modules的基因表达量数据，来识别相关性很高的modules, 从而进行合并，其原理是对modules进行聚类，每个module下的基因表达量是一个二维矩阵，做相关性分析我们只需要一个一维向量就可以了，可以利用PCA分析提取第一主成分来表征原始的矩阵，在WGCNA中，把每个module的表达谱数据对应的一维向量称之为​​Module eigengene E​​。获取一维向量之后，就可以计算相关性，直接用1减去相关性作为距离，进行聚类，代码如下

MEList <- moduleEigengenes( datExpr, colors = dynamicColors)MEs <- MEList$eigengenesMEDiss <- 1-cor(MEs)METree <- hclust(as.dist(MEDiss), method = "average")plot(METree,main = "Clustering of module eigengenes",xlab = "", sub = "")

modules的聚类树示意如下

每个modules的名字用对应的颜色表示，在该聚类数中，分支长度为1减去两个module间的相关系数，在合并modules时，将高相关性的合并为一类，可以指定一个阈值，比如将相关系数大于0.8的合并为一类，在该聚类树中，对应的就是height小于0.2的modules,  对应下图红色的线

可以看到有8个modules都满足条件，在合并时，会将原本属于同一分支的modules直接合并为一个，从图上可以看出，合并后会减少4个modules。合并的代码如下

MEDissThres <- 0.2merge <- mergeCloseModules( datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)mergedColors <- merge$colorsmergedMEs <- merge$newMEsplotDendroAndColors(geneTree,cbind(dynamicColors, mergedColors),c("Dynamic Tree Cut", "Merged dynamic"),dendroLabels = FALSE,hang = 0.03,addGuide = TRUE,guideHang = 0.05)

合并之后的modules 对应的图片如下

最后总结一下，WGCNA在挖掘共表达基因时，首先通过​​TOM​​​统计量将邻接矩阵转换为距离矩阵，然后聚类，利用​​dynamicTreeCut​​的算法识别modules, 最后根据modules之间的相关性，合并modules。

·end·

—如果喜欢，快分享给你的朋友们吧—

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。