正文
正文部分，我想尽可能的从原理出发，讲解WGCNA背后的大致逻辑。我尽可能的在这部分不介绍代码，以免喧宾夺主般的干扰我们对于WGCNA原理的思考。
先简单给个目录，WGCNA背后逻辑得讲解将从一下几个部分进行：
- 数据过滤
- 构建共表达网络
- 构建相似度矩阵
- 构建邻接矩阵的参数选择
- 构建拓扑重叠矩阵（TOM）
- 利用TOM进行聚类
- 识别基因模块
- 合并相似模块
- 模块功能的验证
- 模块与外部性转的关联
- R中的WGCNA可视化
- Cytoscape的网络可视化
数据过滤
对于想要分析的数据，我们要保证缺失值不能过多，如果过多我们必须将这些数据去除。好在我们不用手动去除，只需要交给WGCNA的一个函数即可，这样可以保证数据去除的标准更加客观。
此外，对于多个样本的数据，很有可能某几个样本与其他样本的差异太大，这将严重影响我们后续的分析。因此在前期数据过滤部分，我们首先要进行聚类，查看一下哪些样本不符合我们的预期，提前将这些样本去除。下面是一个简单的聚类结果
我一直觉得聚类可以单独写一个教程具体讲解一下，起码了解一下常见的聚类方法和原理
首先，我们定义一个，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跟着上调
构建邻接矩阵的参数选择
最简单的一个方法，就是如上图中的例子一样，选择一个阈值（thresholding），然后一刀切。超过这个阈值的，我们认为这俩基因存在连接，低于这个阈值的我们认为这俩基因不存在链接。但是这种定义比较粗暴，比如，我定义一个cutoff为0.8，那么0.79算不算两个基因上存在连接呢？
上面那种阈值的定义方法，我们称之为硬阈值（hard thresholding）。为了避免硬阈值所带来的问题，WGCNA分析构建共表达网络时，使用软阈值（soft thresholding）。软阈值简单来说，就是将相似度进行一个幂运算，从而使得结果符合无尺度网络的标准。
该方法让原本值大的更大，小的更小。
无尺度网络模型具有一个特点就是对于部分节点的错误，具有很高程度的鲁棒性 （robustness，可以简单的理解为容错程度）。这是因为该网络具有很多的冗余链接，部分关键节点的错误并不会导致整个网络的调控出现问题。这个特点，与生物体内复杂的网络十分相似，这也是为什么使用无尺度模型的另一个重要原因。
在无尺度网络中，节点的度与该节点出现的概率是负相关的。比如：具有五个连接的节点，其出现的概率要小于有三个连接的节点出现的概率。如下图：
image.png
看的我脑袋疼，相信你们脑袋也疼。所以我们还是举个例子把，这样更方便大家的理解。具体如下图：
image
邻接矩阵中，A与D的关系是1。而在拓扑重叠矩阵中，A和D的关系，不仅需要考虑A和D之间的虚线，还需要A-B-D以及A-C-D。（通过B与C的两条间接关系）
ABBC=11,=1 ACCD=11=1，因此公式中l =1+1 = 2， a =1。l + a 就是A 到D 的直接和间接共表达的总和。而我们可以知道A的k是3 (k )，D的k也是3，因此min{k ,k }=3，因此最终w 的结果为（2+1）/（3+1-1）=1，即我们得到的拓扑重叠矩阵中，第i行j列的元素值为1。
利用TOM进行聚类
有了拓扑异构矩阵（TOM）之后，我们需要进行聚类，从而对基因分块，得到不同的模块（module）。为了更好的展示聚类结果，我们用1减去TOM，从而得到对应的基因不相似性的矩阵（dissTOM）。
接着利用该不相似性的矩阵，进行聚类，聚类结果示例如下：
WGCNA直接将1-TOM得到的矩阵dissTOM，强制转换为距离矩阵。因此越相关的两个基因，其dissTOM中值越小，距离矩阵中对应的距离越近，在图中表示就是越容易聚到一起。
image
在这张图上，树上的每一个叶子节点都是一条小竖线，代表着一个基因。左侧的纵坐标（height）表示两个子节点的距离，如gene A 与 gene B 在0.75处merge到了一起，则这俩基因的距离为0.75，不相似度为0.75。
识别基因模块
得到聚类结果之后，接下来要做的就是找模块（module），即那些连接十分紧密的基因集（共表达程度很高的一类基因），聚类后的每一个簇代表着一个module。
通常将模块树的分支切分开的方法有两种，一种是指定一个高度，将结果分成若干个模块（即：数据预处理部分的红线）。另一种是使用 dynamic branch cut methods，该方法的优点是可以自动化的将module区分开，不再需要手动指定高度，且该方法更为灵活。
使用Dynamic Branch Cut方法处理得到的结果如下：
!](https://cdn.nlark.com/yuque/0/2020/png/1111451/1607994362237-dd3d8d10-6387-4ba7-a026-9b2a3c6cef13.png#align=left&display=inline&height=1115&margin=%5Bobject%20Object%5D&name=image.png&originHeight=1115&originWidth=1499&size=240716&status=done&style=none&width=1499)
合并相似模块
Dynamic Branch Cut作为一种自动切分模块的方法，有时可能会识别出来两个表达谱非常相似的模块，此时就需要我们手动将这些高度共表达的模块合并，方便后续分析。首先需要对不同的模块进行聚类：
Dynamic Branch Cut利用形状参数确定分类的个数
image
这里选择的高度为0.75，对应的相关性为0.75。表示，我们要将相关性大于等于0.75的模块进行合并，最终合并结果如下：
image.png
认真看，可以观察到，紫色和黄色两个色块合并到了一起，绿色和红色两种module合并到了一起。
模块功能的验证
层次聚类的一个缺点是很难确定究竟应该将整个数据集分成几类，尽管我们采用多种分支剪切/模块识别的方法，但是数据的分类问题，是一个开放式的探索问题，并没有一个固定解。即：合理即可。
那么怎么划分模块才是比较合理的呢？通常来说一个共表达模块反应着一个真实的生物学信号，比如通路，比如某个生物学过程。或者反应着噪音，如组织污染，技术偏差，假阳性等。
为了判断一个模块的划分是否合理，通常我们会对一个模块进行GO富集分析，判断是否富集在某一个或者某几个通路。
模块与外部性状关联
在WGCNA中，我们将划分得到的模块，与外部性状相关联，来查看每个模块分别与什么性状高度相关，从而找到我们感兴趣性状的对应模块。这些性状可以是各种指标，比如小鼠的重量，尾长，毛色等；也可以是对应植物的高矮，果实多少，颗粒饱满度等。
但是我们知道计算显著性，通常是两个向量之间作比较。但是一个模块是一个矩阵，矩阵和向量之间，计算相关性，我们该如何分析呢？WGCNA选择使用PCA的方法，计算模块矩阵的主成分，并将其中的PC1定义为eigengenes。
有了eigengennes 我们就可以计算模块与外部性状的相关性了。可视化结果如下：
image
热图中的每个元素代表指定模块和对应性状的显著性，我们重点关注那些高度显著的元素，即红色的方块。
模块与指定性状的相关性我们已经分析得到了，但是一个模块通常有很多基因。我们如何找到最相关的基因呢？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的相关性散点图：
image
从上图我们可以看到，GS和MM还是很相关的，这表明与性状高度相关的基因，通常也是该性状对应模块内比较重要的基因。因此当我们选择感兴趣的重要基因时，推荐选取散点图右上角部分的基因。
R中的WGCNA可视化
我们通常使用热图展示我们得到的权重网络，其中每一行和每一列分别代表着一个基因。使用热图可以给我们展示哪些基因之间更为相似，通常来说都是一个模块内的基因相似度更高。如下图：
热图中方块的颜色越深（红）表示共表达相关性越高，越浅（黄）表示相关性越弱
image
除了看基因之间的共表达程度，我们还可以看看不同模块之间的关系。这时候，就需要我们之前计算得到的eigengenes了，利用eigengenes我们可以得到不同模块的聚类，可视化不同模块的相关性，如下图：
image
热图中，添加了weight这一外部信息，让我们看到该属性最相关的module分别是哪些。
Cytoscape的网络可视化
最后一步，也是大家最为喜闻乐见的一步，就是将我们得到的模块结果保存，并载入到cytoscape中进行网络的可视化。每个顶点代表一个基因，每条边代表两个基因之间存在共表达关系。
Cytoscape如果讲起来，那将又是一堆字，还是再单独开个教程讲解吧
image
完结撒花，原理/思路部分告一段落，下一部分就是WGCNA的实战部分了。实战部分将会跟着原理部分一步一步走下来，分别讲述对应的代码含义和结果。如果不想等待，而且不反感英语，强烈建议看看这个英文教程：https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/
参考资料
- https://abego.cn/2019/05/31/what-is-wgcna-when-to-use/
- https://zh.wikipedia.org/wiki/%E5%9F%BA%E5%9B%A0%E5%85%B1%E8%A1%A8%E8%BE%BE%E7%BD%91%E7%BB%9C
- https://zh.wikipedia.org/wiki/%E5%8A%A0%E6%9D%83%E7%9B%B8%E5%85%B3%E7%BD%91%E7%BB%9C%E5%88%86%E6%9E%90
- https://zhuanlan.zhihu.com/p/36132370
- http://www.wuchangsong.com/?p=198
- http://www.sci666.net/58128.html
- https://blog.csdn.net/weixin_44377608/article/details/90489878
- https://www.jianshu.com/p/fe4d3c77786e
- http://www.sci666.net/58731.html
- [http://www.biotrainee.com/thread-646-1-1.html](
网友评论