WGCNA精讲

最新推荐文章于 2024-07-29 20:44:07 发布

Xiaofei@IDO

最新推荐文章于 2024-07-29 20:44:07 发布

阅读量8.2k

点赞数 24

分类专栏：生信

本文链接：https://blog.csdn.net/nixiang_888/article/details/114867613

版权

生信专栏收录该内容

8 篇文章 10 订阅

订阅专栏

背景

加权基因共表达网络（wighted gene co-expression network nanlysis，WGCNA）将复杂生物过程的基因共表达网络划分为高度相关的几个特征模块，其代表着几组高度协同变化的基因集，并可将模块与特定的临床特征建立关联，从中寻找发挥关键功能的基因，帮助识别参与特定生物学过程的潜在机制以及探索候选生物标志物（Langfelder and Horvath, 2008）
WGCNA基于两个假设：（1）相似表达模式的基因可能存在共调控、功能相关或处于同一通路，（2）基因网络符合无标度分布。基于这两点，可以将基因网络根据表达相似性划分为不同的模块进而找出枢纽基因。
现在我来详细讲解一下WGCNA的具体细节。

WGCNA实现

1. 数据过滤

对于想要分析的数据，我们要保证缺失值不能过多，如果过多我们必须将这些数据去除。好在我们不用手动去除，只需要交给WGCNA的一个函数即可，这样可以保证数据去除的标准更加客观。
此外，对于多个样本的数据，很有可能某几个样本与其他样本的差异太大，这将严重影响我们后续的分析。因此在前期数据过滤部分，我们首先要进行聚类，查看一下哪些样本不符合我们的预期，提前将这些样本去除。下面是一个简单的聚类结果.
. clustring
其中的红线，将样本中的outlier与其他样本区分开。这个cut-off的选取可以自己设置，将自己认为的outlier排除在外即可。

2. 构建共表达调控网络

2.1 构建基因间的相似度矩阵

一个网络关系通常是由一个邻接矩阵所定义的，在这个矩阵中，每个元素的取值范围都在0-1之间，代表两个节点之间的连接强度。比如矩阵中第 i 行，第 j 列的元素我们称之为 a ，它代表的就是 gene i 与gene j 之间的连接强度。
为了得到这个邻接矩阵，我们需要定义一个中间变量，共表达相似性矩阵。以下图为例：
co-expression
首先，我们定义一个，3样本10基因的表达矩阵。每一行是一个基因，每一列代表一个样本，每个元素代表特定样本中某基因的表达值。
有了该矩阵，我们计算Pearson相关系数，从而得到了一个10X10的相似度矩阵。这个矩阵中每个元素都代表着两个基因之间的相似度，如第三行第四列的元素值为0.99，代表Gene 3和Gene 4的相似度是0.99。

使用皮尔逊相关系数具有一个缺点，对于异常值十分的敏感。为了解决该问题，我们通常采用Jackknife相关系数，即先使用Jackknife方法进行重采样，再进行相关性的分析（先抽样，再相关性分析）

图示中的例子采用的是硬阈值（hard thresholding）的方式构建的邻接矩阵，即相似度超过0.8的才认为这俩基因之间存在连接，设置为1。而不超过0.8的则表示二者没有连接，设置为0。

硬阈值就是我们常用的一刀切的阈值

我们需要注意的是，这里定义相似度矩阵时，使用的公式是|r(Gi, Gj)| 。即，对相关系数进行了绝对值的计算，这种方式得到的是unsigned的，即不知道是正相关还是负相关。

公式不重要，重要的是这里有一个绝对值符号；r(Gi, Gj)的取值范围为-1到1

除此之外，我们还可以使用这个公式|(1+ r(Gi, Gj))/2|进行计算。这样的话，我们就可以将负值（负相关）转换为正值，而原有的正值变得更大。利用这种方式，在后续分析的时候，我们更为关注的是那些正相关的基因。

正相关反映到基因调控中，代表正向调控，gene A上调会导致gene B跟着上调

2.2 构建邻接矩阵的参数 cutoff-value 的选择

最简单的一个方法，就是如上图中的例子一样，选择一个阈值（thresholding），然后一刀切。超过这个阈值的，我们认为这俩基因存在连接，低于这个阈值的我们认为这俩基因不存在链接。但是这种定义比较粗暴，比如，我定义一个cutoff为0.8，那么0.79算不算两个基因上存在连接呢？
上面那种阈值的定义方法，我们称之为硬阈值（hard thresholding）。为了避免硬阈值所带来的问题，WGCNA分析构建共表达网络时，使用软阈值（soft thresholding）。软阈值简单来说，就是将相似度进行一个幂运算，从而使得结果符合无尺度网络的标准。

该方法使强相关性变的更强，弱相关性变的更弱

无尺度网络模型具有一个特点就是对于部分节点的错误，具有很高程度的鲁棒性（robustness，可以简单的理解为容错程度）。这是因为该网络具有很多的冗余链接，部分关键节点的错误并不会导致整个网络的调控出现问题。这个特点，与生物体内复杂的网络十分相似，这也是为什么使用无尺度模型的另一个重要原因。
在无尺度网络中，节点的度与该节点出现的概率是负相关的。比如：具有五个连接的节点，其出现的概率要小于有三个连接的节点出现的概率。如下图：
scale

其中，横坐标表示某个节点连通性的log值，纵坐标表示该节点出现的概率的log值。当我们选择的幂指数 (power)，让数据十分符合无尺度网络时，其 $R^2$ 就会十分的接近于1。在这张图上，分布大致遵循一条直线，这被称为近似无尺度拓扑。

连通性，简单的理解为节点上连接的边的个数

我们在使用WGCNA分析时，通常会产生一个如下所示的图，辅助我们进行幂指数 (power) 的选择：

幂的英文即power，有的教程中，会使用β指代该值。

power
两张图的横坐标，均代表幂指数 (power) 的取值。左图的纵坐标是我们上面所讲的R2，越高越好。而右图的纵坐标代表所有节点的平均连通数，越低越好。通常我们对于幂指数 (power) 的选择主要参考左图，右图用于辅助参考。

无尺度网络中，大部分节点的连通度较低，因此平均连通数越低越好

此外，在所有达到设定阈值的幂指数中，幂指数越小越好

WGCNA软阈值power选择标准是什么？
答：一般选取首次使无标度网络结构 R² 达到0.8的power取值，然而有时候无向网络在power小于15或有向网络power小于30内，没有一个power值可以使无标度网络图谱结构 R² 达到0.8，平均连接度较高如在100以上，可能是由于部分样品与其他样品差别太大。这可能由批次效应、样品异质性或实验条件对表达影响太大等造成。可以通过绘制样品聚类查看分组信息和有无异常样品。如果这确实是由有意义的生物变化引起的，也可以使用经验power值。

经验power值

if (is.na(power)){
	power = ifelse(nSamples<20, ifelse(type == "unsigned", 9, 18),
		ifelse(nSamples<30, ifelse(type == "unsigned", 8, 16),
			ifelse(nSamples<40, ifelse(type == "unsigned", 7, 14),
				ifelse(type == "unsigned", 6, 12))
		)
	)
}

2.3 构建拓扑重叠矩阵 - TOM

WGCNA认为邻接矩阵是不够的，基因间的相似性应该在表达和网络拓扑水平上体现。于是，又利用邻接矩阵生成了一个拓扑重叠矩阵/TOM（topological overlap matrix），此矩阵的构建还能够尽可能地最小化噪音和假阳性的影响。

gene A 与 gene B之间的关系，不仅考虑AB直接的关系，还考虑可能的多个间接关系，这样任何一个关系存在假阳性都不至于对AB之间的关系造成巨大的影响。有点无尺度网络那味了~

比如gene A 与 gene B，邻接矩阵关注的是这俩基因之间直接的共表达关系，而拓扑重叠矩阵还考虑中间因素的影响，不止关注gene A对B的调控，还要关注geneA对B的次级影响。比如gene A调控 gene C，gene C又调控gene B。

那么拓扑重叠矩阵怎么计算得到呢？具体看如下公式：
formular
看的我脑袋疼，相信你们脑袋也疼。所以我们还是举个例子把，这样更方便大家的理解。具体如下图：
fig
邻接矩阵中，A与D的关系是1。而在拓扑重叠矩阵中，A和D的关系，不仅需要考虑A和D之间的虚线，还需要A-B-D以及A-C-D。通过B与C的两条间接关系。
ABBD=1， ACCD=1，因此公式中L=1+1 = 2， a =1。l + a 就是A 到D 的直接和间接共表达的总和。而我们可以知道A的k是3 (ki)，D的k也是3(kj)，因此min{ki ,kj}=3，因此最终w 的结果为（2+1）/（3+1-1）=1，即我们得到的拓扑重叠矩阵中，第i行j列的元素值为1。

2.4 基于TOM矩阵进行聚类

有了拓扑异构矩阵（TOM）之后，我们需要进行聚类，从而对基因分块，得到不同的模块（module）。为了更好的展示聚类结果，我们用1减去TOM，从而得到对应的基因不相似性的矩阵（dissTOM）。
接着利用该不相似性的矩阵，进行聚类，聚类结果示例如下：

WGCNA直接将1-TOM得到的矩阵dissTOM，强制转换为距离矩阵。因此越相关的两个基因，其dissTOM中值越小，距离矩阵中对应的距离越近，在图中表示就是越容易聚到一起。

TOM

在这张图上，树上的每一个叶子节点都是一条小竖线，代表着一个基因。左侧的纵坐标（height）表示两个子节点的距离，如gene A 与 gene B 在0.75处merge到了一起，则这俩基因的距离为0.75，不相似度为0.75。

2.5 识别基因共表达模块

得到聚类结果之后，接下来要做的就是找模块（module），即那些连接十分紧密的基因集（共表达程度很高的一类基因），聚类后的每一个簇代表着一个module。
通常将模块树的分支切分开的方法有两种，一种是指定一个高度，将结果分成若干个模块（即：数据预处理部分的红线）。另一种是使用 dynamic branch cut methods，该方法的优点是可以自动化的将module区分开，不再需要手动指定高度，且该方法更为灵活。
使用Dynamic Branch Cut方法处理得到的结果如下：
在这里插入图片描述

2.6 合并相似模块

Dynamic Branch Cut作为一种自动切分模块的方法，有时可能会识别出来两个表达谱非常相似的模块，此时就需要我们手动将这些高度共表达的模块合并，方便后续分析。首先需要对不同的模块进行聚类：

Dynamic Branch Cut利用形状参数确定分类的个数

在这里插入图片描述
这里选择的高度为0.25，对应的相关性为0.75。表示，我们要将相关性大于等于0.75的模块进行合并，最终合并结果如下：

认真看，可以观察到，紫色和黄色两个色块合并到了一起，绿色和红色两种module合并到了一起。

3 模块功能验证

层次聚类的一个缺点是很难确定究竟应该将整个数据集分成几类，尽管我们采用多种分支剪切/模块识别的方法，但是数据的分类问题，是一个开放式的探索问题，并没有一个固定解。即：合理即可。
那么怎么划分模块才是比较合理的呢？通常来说一个共表达模块反应着一个真实的生物学信号，比如通路，比如某个生物学过程。或者反应着噪音，如组织污染，技术偏差，假阳性等。
为了判断一个模块的划分是否合理，通常我们会对一个模块进行GO富集分析，判断是否富集在某一个或者某几个通路。

4 模块与外部关联性状

在WGCNA中，我们将划分得到的模块，与外部性状相关联，来查看每个模块分别与什么性状高度相关，从而找到我们感兴趣性状地对应模块。这些性状可以是各种指标，比如小鼠的重量，尾长，毛色等；也可以是对应植物的高矮，果实多少，颗粒饱满度等。

但是我们知道计算显著性，通常是两个向量之间作比较。但是一个模块是一个矩阵，矩阵和向量之间，计算相关性，我们该如何分析呢？

WGCNA选择使用PCA的方法，计算模块矩阵的主成分，并将其中的PC1定义为eigengenes。

有了eigengennes 我们就可以计算模块与外部性状的相关性了。可视化结果如下：
在这里插入图片描述
热图中的每个元素代表指定模块和对应性状的显著性，我们重点关注那些高度显著的元素，即红色的方块。

模块与指定性状的相关性我们已经分析得到了，但是一个模块通常有很多基因。我们如何找到最相关的基因呢？

WGCNA通过GS和MM给出了一个合理的解决办法。

首先我们将模块中的每个基因与eigengenes进行相关性分析，得到的结果我们称之为module membership（MM）。当结果越接近于0，则我们认为该基因与其所在的模块越不相关。而结果越接近于1或者-1，则我们认为该基因与模块基因高度相关，正负表示其是正相关还是负相关。
此外，module membership与模块内的连通性是高度相关的。其中，高连通性的hub genes倾向于有更高的module membership 值。有了这个因果关系，我们不需要通过cytoscape绘制网络图，也可以找到hub genes了。
模块内的基因，我们不仅想知道它与模块的相关性，还想知道它与模块对应的性状的相关性。于是WGCNA定义了一个新的变量GS (gene significance)，即计算模块内gene 与对应性状的显著性。

下图展示了 module membership的相关性散点图：
在这里插入图片描述从上图我们可以看到，GS和MM还是很相关的，这表明与性状高度相关的基因，通常也是该性状对应模块内比较重要的基因。因此当我们选择感兴趣的重要基因时，推荐选取散点图右上角部分的基因。