WGCNA分析是如何找出基因模块的?

一个WGCNA网络的分析过程包括四个方面。

1.确定软阈值,构建邻接矩阵。
2.把邻接矩阵变成拓扑重叠矩阵
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 = FALSEhang = 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)

这是合并前后的对比。

找到了模块,我们就要去探索这个模块里面的基因了,
同时跟样本性状的关系也要提上日程。

生物信息学

第一次见3D的火山图

2020-6-25 20:24:40

生物信息学

如何找到基因组中的模糊碱基呢

2020-6-25 20:27:51

声明 本网部分文章源于互联网,转载出于传递更多信息和学习之目的,并不意味着赞同其观点。
如转载稿涉及版权等问题,请立即联系管理员,我们会予以修改、删除相关文章,请留言反馈
When your legal rights are being violated, please send an email to: sci666net@qq.com
0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
购物车
优惠劵
今日签到
有新私信 私信列表
搜索