github-actions[bot]
commented
4 years ago WGCNA分析是如何找出基因模块的？ by 果子学生信
一个WGCNA网络的分析过程包括四个方面。
1.确定软阈值，构建邻接矩阵。
WGCNA的分析中为什么要挑选软阈值？
2.把邻接矩阵变成拓扑重叠矩阵
WGCNA有了相关性矩阵为什么还要计算拓扑矩阵？
3.根据拓扑重叠矩阵识别基因模块
这是今天要讲的。
4.探索感兴趣的基因模块：包括性状关联，富集分析，寻找hub基因
WGCNA如何识别模块
靠的是层次聚类。
下面用数据演示一下
首先从表达矩阵计算出拓扑重叠矩阵TOM
library(WGCNA)
### 运行多线程，可能mac用户不适合
enableWGCNAThreads()
### 加载数据
load(file = "FemaleLiver-01-dataInput.RData")
## 确定软阈值
softPower = 6
## 邻接矩阵
adjacency = adjacency(datExpr, power = softPower)
## TOM矩阵
TOM = TOMsimilarity(adjacency)
dissTOM = 1-TOM
获得TOM矩阵(示意图)
现在我把这个数据层次聚类一下
geneTree = hclust(as.dist(dissTOM), method = "average")
层次聚类用的函数是hclust
其中methods有很多种，比如average，single，complete
statquest讲了他们的区别。
在计算两个cluster距离的时候，single选择的离得最近的点，complete选择的是最远的点
然后average是计算每个点然后取平均值。
而WGCNA中这里选择的是average。
现在把数据可视化
sizeGrWindow(12,9)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",labels = FALSE, hang = 0.04);
这一个个像羽毛一样的分支，就是基因模块。
如果给他一条线来切割(红色的横线)
那么这个树就切割成了四个分支。
就对应到四个基因模块。
这种切割分支的方法叫作静态切割。
依靠眼睛看，不方便，更重要的是，他有很大的缺点。
第一，模块太粗糙
比如，第二个模块中，明显还可以继续切割，但是没办法做了。
第二，有些数据会完全失效
像这样的数据，没办法找到一个合适切割值，让主要分支都在其中。
不提供解决方案的提建议都是耍流氓。
你批评一个方法不好时，肯定要有解决方案。
动态切割
所以作者给了他自己的解决方案。Dynamic Tree Cut 动态切割。
在WGCNA中的他是从下往上切割的，先得到小的分支，如果分支相同就合并，不相同就继续提高切割值
找到新的分支。
如果碰到一些模棱两可的基因怎么办？
这时候就用到的K均值聚类的升级版，PAM方法，区别是PAM对异常值不敏感。
把这些模糊的基因分配给最近的模块。
因为这种方法融合了PAM算法，所以又称为动态杂合切割法(Dynamic hybird)
看看这些方案的使用效果吧
静态的切割（Static）,纯PAM法，动态的(Dynamic), 动态联合PAM(Dynamic hybird)
在图中，灰色代表没有被归到模块的基因，我们发现纯PAM方法灰色特别少，
说明，这个方法会把很多没有作用的基因强行归类，所以要结合动态起个来用。
本来我要拆解这个函数的，但是，源码800多行，hold不住。就直接使用了
minModuleSize = 30
dynamicMods = cutreeDynamic(dendro = geneTree,
distM = dissTOM,
deepSplit = 2,
pamStage = TRUE,
pamRespectsDendro = FALSE,
minClusterSize = minModuleSize)
minModuleSize规定最小的模块不小于30
dendro = geneTree, distM = dissTOM 输入的是前面的树和TOM矩阵
deepSplit 是切割强度，越大得到的模块越多
pamStage = TRUE, pamRespectsDendro = FALSE,
这是对PAM方法的限制，pamStage = TRUE是使用PAM方法，pamRespectsDendro = FALSE，是说PAM不会考虑树的感受
翻译就是，如果设置为TRUE，那么PAM在分配模棱两可基因的时候会把他们限定在同一个分支中，这是残血版
如果设置为FALSE，PAM方法就是满血版本了，按照自己的方法去分配。
结果就是，有很多灰色的不属于任何模块的基因，也找到了归属，有时候会导致结果不好解释。
该函数最终返回的结果是各个基因的模块编号
table(dynamicMods)
结果如下
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
88 614 316 311 257 235 225 212 158 153 121 106 102 100 94 91 78 76 65 58 58 48 34
总共找到22个模块，0里面的88个基因不属于任何一个模块，就是垃圾箱。labels2colors这个函数可以把数字转换为颜色。
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
合并模块
好了现在模块就找到了。
假如你觉得模块太多了，还可以把相近的合并
以0.25作为基准(相关性是0.75)，可以看到基准以下的模块可以合并，表示他们很相近
用mergeCloseModules函数就可以搞定了
merge = mergeCloseModules(datExpr,
dynamicColors,
cutHeight = 0.25)
这是合并前后的对比。
找到了模块，我们就要去探索这个模块里面的基因了，
同时跟样本性状的关系也要提上日程。
https://mp.weixin.qq.com/s/s3OhrH7x4Psnlf6ENJm0YQ