正文
正文部分，我想尽可能的从原理出发，讲解WGCNA背后的大致逻辑。我尽可能的在这部分不介绍代码，以免喧宾夺主般的干扰我们对于WGCNA原理的思考。
先简单给个目录，WGCNA背后逻辑得讲解将从一下几个部分进行：
- 数据过滤
- 构建共表达网络
- 构建相似度矩阵
- 构建邻接矩阵的参数选择
- 构建拓扑重叠矩阵（TOM）
- 利用TOM进行聚类
- 识别基因模块
- 合并相似模块
- 模块功能的验证
- 模块与外部性转的关联
- R中的WGCNA可视化
- Cytoscape的网络可视化
数据过滤
对于想要分析的数据，我们要保证缺失值不能过多，如果过多我们必须将这些数据去除。好在我们不用手动去除，只需要交给WGCNA的一个函数即可，这样可以保证数据去除的标准更加客观。
此外，对于多个样本的数据，很有可能某几个样本与其他样本的差异太大，这将严重影响我们后续的分析。因此在前期数据过滤部分，我们首先要进行聚类，查看一下哪些样本不符合我们的预期，提前将这些样本去除。下面是一个简单的聚类结果
我一直觉得聚类可以单独写一个教程具体讲解一下，起码了解一下常见的聚类方法和原理
其中的红线，将样本中的outlier与其他样本区分开。这个cut-off的选取可以自己设置，将自己认为的outlier排除在外即可。
我相信即使不懂聚类的原理，也能看出来左上的那个分支所代表的样本，和其他样本相差甚远
构建共表达调控网络
这一部分是WGCNA的重头戏，不得不好好说道一下……
构建相似度矩阵
一个网络关系通常是由一个邻接矩阵所定义的，在这个矩阵中，每个元素的取值范围都在0-1之间，代表两个节点之间的连接强度。比如矩阵中第i行，第j列的元素我们称之为a ，它代表的就是gene i 与gene j之间的连接强度。
为了得到这个邻接矩阵，我们需要定义一个中间变量，共表达相似性矩阵。以下图为例：
首先，我们定义一个，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，可以简单的理解为容错程度）。这是因为该网络具有很多的冗余链接，部分关键节点的错误并不会导致整个网络的调控出现问题。这个特点，与生物体内复杂的网络十分相似，这也是为什么使用无尺度模型的另一个重要原因。
在无尺度网络中，节点的度与该节点出现的概率是负相关的。比如：具有五个连接的节点，其出现的概率要小于有三个连接的节点出现的概率。如下图：
其中，横坐标表示某个节点连通性的log值，纵坐标表示该节点出现的概率的log值。当我们选择的幂指数 (power)，让数据十分符合无尺度网络时，其R 就会十分的接近于1。在这张图上，分布大致遵循一条直线，这被称为近似无尺度拓扑。
连通性，简单的理解为节点上连接的边的个数
我们在使用WGCNA分析时，通常会产生一个如下所示的图辅助我们进行幂指数 (power) 的选择：
幂的英文即power，有的教程中，会使用β指代该值。
两张图的横坐标，均代表幂指数 (power) 的取值。左图的纵坐标是我们上面所讲的R ，越高越好。而右图的纵坐标代表所有节点的平均连通数，越低越好。通常我们对于幂指数 (power) 的选择主要参考左图，右图用于辅助参考。
无尺度网络中，大部分节点的连通度较低，因此平均连通数越低越好。
此外，在所有达到设定阈值的幂指数中，幂指数越小越好，不然我们直接选择1000，或者一个非常大的值它不香么？
构建拓扑重叠矩阵(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。
既然考虑了次级调控，为什么不考虑更次一级的调控呢？这里将其称之为三级调控。如，A不仅自己调控B，还通过调控C，C调控D，D又调控B的方式，间接的调控B。两个原因，一个可能是计算量太大，太耗费计算资源了。另一个可能是，三级调控对基因B的影响基本可以忽略不计了，这样也就可以解释为什么没有四级调控，五级调控……
那么拓扑重叠矩阵怎么计算得到呢？具体看如下公式：
看的我脑袋疼，相信你们脑袋也疼。所以我们还是举个例子把，这样更方便大家的理解。具体如下图：
邻接矩阵中，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中值越小，距离矩阵中对应的距离越近，在图中表示就是越容易聚到一起。
在这张图上，树上的每一个叶子节点都是一条小竖线，代表着一个基因。左侧的纵坐标（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利用形状参数确定分类的个数
这里选择的高度为0.75，对应的相关性为0.75。表示，我们要将相关性大于等于0.75的模块进行合并，最终合并结果如下：
认真看，可以观察到，紫色和黄色两个色块合并到了一起，绿色和红色两种module合并到了一起。
模块功能的验证
层次聚类的一个缺点是很难确定究竟应该将整个数据集分成几类，尽管我们采用多种分支剪切/模块识别的方法，但是数据的分类问题，是一个开放式的探索问题，并没有一个固定解。即：合理即可。
那么怎么划分模块才是比较合理的呢？通常来说一个共表达模块反应着一个真实的生物学信号，比如通路，比如某个生物学过程。或者反应着噪音，如组织污染，技术偏差，假阳性等。
为了判断一个模块的划分是否合理，通常我们会对一个模块进行GO富集分析，判断是否富集在某一个或者某几个通路。
模块与外部性状关联
在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还是很相关的，这表明与性状高度相关的基因，通常也是该性状对应模块内比较重要的基因。因此当我们选择感兴趣的重要基因时，推荐选取散点图右上角部分的基因。
R中的WGCNA可视化
我们通常使用热图展示我们得到的权重网络，其中每一行和每一列分别代表着一个基因。使用热图可以给我们展示哪些基因之间更为相似，通常来说都是一个模块内的基因相似度更高。如下图：
热图中方块的颜色越深（红）表示共表达相关性越高，越浅（黄）表示相关性越弱
除了看基因之间的共表达程度，我们还可以看看不同模块之间的关系。这时候，就需要我们之前计算得到的eigengenes了，利用eigengenes我们可以得到不同模块的聚类，可视化不同模块的相关性，如下图：
热图中，添加了weight这一外部信息，让我们看到该属性最相关的module分别是哪些。
Cytoscape的网络可视化
最后一步，也是大家最为喜闻乐见的一步，就是将我们得到的模块结果保存，并载入到cytoscape中进行网络的可视化。每个顶点代表一个基因，每条边代表两个基因之间存在共表达关系。
Cytoscape如果讲起来，那将又是一堆字，还是再单独开个教程讲解吧
完结撒花，原理/思路部分告一段落，下一部分就是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](
