一文看懂 WGCNA 分析

生信技能树2019-10-12

16119

基本概念

WGCNA 其译为加权基因共表达网络分析。该分析方法旨在寻找协同表达的基因模块 (module)，并探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。

适用于复杂的数据模式，推荐 5 组 (或者 15 个样品) 以上的数据。

一般可应用的研究方向有：

不同器官或组织类型发育调控；
同一组织不同发育调控；
非生物胁迫不同时间点应答；
病原菌侵染后不同时间点应答。

基本原理

从方法上来讲，WGCNA 分为表达量聚类分析和表型关联两部分，主要包括基因之间相关系数计算、基因模块的确定、共表达网络、模块与性状关联四个步骤。

第一步计算任意两个基因之间的相关系数（Person Coefficient）。为了衡量两个基因是否具有相似表达模式，一般需要设置阈值来筛选，高于阈值的则认为是相似的。但是这样如果将阈值设为 0.8，那么很难说明 0.8 和 0.79 两个是有显著差别的。因此，WGCNA 分析时采用相关系数加权值，即对基因相关系数取 N 次幂，使得网络中的基因之间的连接服从无尺度网络分布 (scale-freenetworks)，这种算法更具生物学意义。

第二步通过基因之间的相关系数构建分层聚类树，聚类树的不同分支代表不同的基因模块，不同颜色代表不同的模块。基于基因的加权相关系数，将基因按照表达模式进行分类，将模式相似的基因归为一个模块。这样就可以将几万个基因通过基因表达模式被分成了几十个模块，是一个提取归纳信息的过程。

WGCNA 术语

权重 (weghted)

基因之间不仅仅是相关与否，还记录着它们的相关性数值，数值就是基因之间的联系的权重 (相关性)。

Module

模块 (module)：表达模式相似的基因分为一类，这样的一类基因成为模块；

Eigengene

Eigengene（eigen +‎ gene）：基因和样本构成的矩阵，https://en.wiktionary.org/wiki/eigengene

Adjacency Matrix

邻近矩阵：是图的一种存储形式，用一个一维数组存放图中所有顶点数据；用一个二维数组存放顶点间关系（边或弧）的数据，这个二维数组称为邻接矩阵；在 WGCNA 分析里面指的是基因与基因之间的相关性系数矩阵。如果用了阈值来判断基因相关与否，那么这个邻近矩阵就是 0 / 1 矩阵，只记录基因相关与否。但是 WGCNA 没有用阈值来卡基因的相关性，而是记录了所有基因之间的相关性。

Topological Overlap Matrix (TOM)

WGNA 认为基因之间的简单的相关性不足以计算共表达，所以它利用上面的邻近矩阵，又计算了一个新的邻近矩阵。一般来说，TOM 就是 WGCNA 分析的最终结果，后续的只是对 TOM 的下游注释。

下游分析

得到模块之后的分析有：

模块的功能富集
模块与性状之间的相关性
模块与样本间的相关系数

挖掘模块的关键信息：

找到模块的核心基因
利用关系预测基因功能

代码示例

其中第一步数据准备反而是最复杂的，取决于大家的 R 语言水平，这个数据 GSE48213 -wgcna-input.RData 我已经保存下来咯，如果大家不会做，又想体验一下这个 WGCNA 流程，就可以直接 load 我保存好的数据文件即可。

step1: 输入数据的准备

这里主要是表达矩阵，如果是芯片数据，那么常规的归一化矩阵即可，如果是转录组数据，最好是 RPKM/TPM 值或者其它归一化好的表达量。然后就是临床信息或者其它表型，总之就是样本的属性。

为了保证后续脚本的统一性，表达矩阵统一用 datExpr 标识，临床信息统一用 datTraits 标识。(PS: 如果你 R 语言很差，变量名不要轻易修改)

上面代码里面的 rpkm 就是我们的转录组数据的表达矩阵，以 rpkm 为单位。而 datTraits 就是所有样本对应的表型信息。需要自己制作，这个是学习 WGCNA 的基础，本次实例代码都是基于这两个数据。至于如何做出上面代码的两个例子，取决于大家自己的项目，我这里给出自己的代码，仅供参考哈！

很明显，这个数据 GSE48213 -wgcna-input.RData 我已经保存下来咯，如果大家不会做，又想体验一下这个 WGCNA 流程，那么可以找我 email 求取这个数据哦。我的邮箱是 jmzeng1314@163.com

我给大家演示的示例数据大概是下面这个样子：

这个数据集里面的 56 种细胞系被分成了 5 组，如果要分开两两做差异分析，有 10 种组合，也就是说需要做 10 次差异分析，每个差异分析结果都需要去注释，会比较麻烦，这个时候 WGCNA 就派上用场啦。当然，如果你一定要去做差异分析，我也给你代码：https://github.com/jmzeng1314 /my-R/blob/master/ 10 -RNA-seq- 3 -groups/hisat2_mm10_htseq.R

实际上多个分组，差异分析策略是非常个性化的，比如：

https://mp.weixin.qq.com/s/hc6JkKxyelc7b1M1MRiHRQ

step2: 确定最佳 beta 值

选择合适「软阀值（soft thresholding power）」beta，同样的，也是使用教程标准代码即可：

关键就是理解 pickSoftThreshold 函数及其返回的对象，最佳的 beta 值就是 sft$powerEstimate

最佳 beta 值

参数 beta 取值默认是 1 到 30，上述图形的横轴均代表权重参数β，左图纵轴代表对应的网络中 log(k) 与 log(p(k)) 相关系数的平方。相关系数的平方越高，说明该网络越逼近无网路尺度的分布。右图的纵轴代表对应的基因模块中所有基因邻接函数的均值。最佳的 beta 值就是 sft$powerEstimate，已经被保存到变量了，不需要知道具体是什么，后面的代码都用这个即可，在本例子里面是 6。

即使你不理解它，也可以使用代码拿到合适「软阀值（soft thresholding power）」beta 进行后续分析。

step3：一步法构建共表达矩阵

有了表达矩阵和估计好的最佳 beta 值，就可以直接构建共表达矩阵了。

所有的核心就在这一步，把输入的表达矩阵的几千个基因组归类成了几十个模块。大体思路：计算基因间的邻接性，根据邻接性计算基因间的相似性，然后推出基因间的相异性系数，并据此得到基因间的系统聚类树。然后按照混合动态剪切树的标准，设置每个基因模块最少的基因数目为 30。

根据动态剪切法确定基因模块后，再次分析，依次计算每个模块的特征向量值，然后对模块进行聚类分析，将距离较近的模块合并为新的模块。

step4: 模块可视化

这里用不同的颜色来代表那些所有的模块，其中灰色默认是无法归类于任何模块的那些基因，如果灰色模块里面的基因太多，那么前期对表达矩阵挑选基因的步骤可能就不太合适。

基因的模块可视化

这里的重点就是 plotDendroAndColors 函数，它接受一个聚类的对象，以及该对象里面包含的所有个体所对应的颜色。比如对表达矩阵进行 hclust 之后，加上表达矩阵里面所有样本的分组信息对应的颜色，也是可以用 plotDendroAndColors 函数可视化的，比如下面的样品图：

上面给样本进行聚类的代码可以不运行，其实跟 WGCNA 本身关系不大。

样本的聚类可视化

可以看到这些乳腺癌的细胞系的表达谱聚类情况并不是完全与其分类匹配，所以仅仅是根据样本的分组信息做差异分析并不完全准确。

step5: 模块和性状的关系

通过模块与各种表型的相关系数，可以很清楚的挑选自己感兴趣的模块进行下游分析了。这个图就是把 moduleTraitCor 这个矩阵给用热图可视化一下。

模块和性状的关系

因为一些历史遗留问题，这个图片缺乏 X 轴的标记。

从上图已经可以看到跟乳腺癌分类相关的基因模块了，包括"Basal" "Claudin-low" "Luminal" "Non-malignant" "unknown" 这 5 类所对应的不同模块的基因列表。可以看到每一种乳腺癌都有跟它强烈相关的模块，可以作为它的表达 signature，模块里面的基因可以拿去做下游分析。我们看到 Luminal 表型跟棕色的模块相关性高达 0.86，而且极其显著的相关，所以值得我们挖掘，这个模块里面的基因是什么，为什么如此的相关呢？

step6: 感兴趣性状的模块的具体基因分析

性状跟模块虽然求出了相关性，可以挑选最相关的那些模块来分析，但是模块本身仍然包含非常多的基因，还需进一步的寻找最重要的基因。所有的模块都可以跟基因算出相关系数，所有的连续型性状也可以跟基因的表达值算出相关系数。主要参考资料：PDF document, R script 如果跟性状显著相关基因也跟某个模块显著相关，那么这些基因可能就非常重要。

首先计算模块与基因的相关性矩阵