跟风学一波WGCNA

听说WGCNA是一个看上去比较厉害的转录组分析方法，近期又多次看到了相关的内容，所以跟风学习一下。对于我这种业余学生信的来说确实复杂了些，断断续续看了很久，也只是勉强按流程跟着跑一遍。这篇只有简单的原理和介绍，没有代码。

一、WGCNA简介

对于两组数据，比如实验组和对照组的测序数据，如果我们要找到与实验组条件相关的通路或者基因，可以做差异分析然后进行GO分析、KEGG分析，或者直接使用GSEA分析。那么如果我有多个分组怎么找出每个组的特异分子呢？假设我有10个分组共60个样品，如果两两比较需要比45次，这样不够高效。实际上临床性状可能更加复杂，有多种性状混合。

加权基因共表达网络分析（Weighted Gene Co-Expression Network Analysis，WGCNA），是一种适合进行多样本复杂数据分析的工具，通过计算基因间表达关系，鉴定表达模式相似的基因集合，解析基因集合与样品表型之间的联系，绘制基因集合中基因之间的调控网络并鉴定关键调控基因。WGCNA将复杂生物过程的基因共表达网络划分为高度相关的几个特征模块，其代表着几组高度协同变化的基因集，并可将模块与特定的临床特征建立关联，从中寻找发挥关键功能的基因，帮助识别参与特定生物学过程的潜在机制以及探索候选生物标志物。

相比于只关注差异表达的基因，WGCNA利用数千或近万个变化最大的基因或全部基因的信息识别感兴趣的基因集，并与表型进行显著性关联分析。一是充分利用了信息，二是把数千个基因与表型的关联转换为数个基因集与表型的关联，免去了多重假设检验校正的问题。

二、WGCNA原理

（1）基本假设

WGCNA基于两个假设：（1）相似表达模式的基因可能存在共调控、功能相关或处于同一通路，（2）基因网络符合无尺度分布。

【一个网络图主要包含节点（node）和连接线（edge）两种成分，如果两个节点之间存在相互关系，就画一条线连接起来。网络可以分为有方向的和没有方向的，对于基因网络而言，一般我们预测到的都是没有方向性的。度（degree）代表着一个节点和多少个其他节点的有连接，如果有方向就会分为in-degree和out-degree。随机网络（random network）是指在网络中的所有点之间的关联性都差不多，可以看做是随机连接而成。无尺度网络（scale-free network）具有严重的异质性，其各节点之间的连接状况（度数）具有严重的不均匀分布性：网络中少数称之为Hub点的节点拥有极其多的连接，而大多数节点只有很少量的连接。少数Hub点对无标度网络的运行起着主导的作用。从广义上说，无标度网络的无标度性是描述大量复杂系统整体上严重不均匀分布的一种内在性质。】

无尺度网络服从幂律分布

（2）计算两个基因的相关系数

传统方法上，描述两个基因间的关联程度，可通过计算表达值间的Pearson、Spearman等相关系数获得。为了构建关联网络，通常指定一个筛选阈值，如相关系数大于0.8以上，作为两个基因间具有强关联程度的依据。随后在网络中，节点表示基因，若两个基因间的相关系数大于指定的阈值，则会以边连接以表示两者间可能存在相互作用。但是基于固定阈值法的缺点在于，阈值是人为定义的，将会忽略很多潜在关联。例如，0.79就是不相关吗？同时，这种一刀切的方法也会丢失基因的变化趋势信息，将难以在网络中描述相关性的强弱关系。

经典的共表达网络使用共表达相关性构建网络。相关性测量的值介于-1（完全负相关）到1（完全正相关）之间。在无方向网络（unsigned）中，使用的是绝对相关系数，两个负相关的基因也会被认为共表达，这就导致负相关的基因也会被group到一起。有方向的网络（signed）采取的是把相关性值scale到0和1之间，其中0-0.5代表负相关，0.5-1代表正相关。共表达网络定义为无向的、加权的基因网络。

【Pearson相关系数是用于表示相关性大小的最常用指标，数值介于-1~1之间，越接近0相关性越低，越接近-1或1相关性越高。正负号表明相关方向，正号为正相关、负号为负相关。适用于两个正态分布的连续变量。两个变量X和Y之间的皮尔森相关系数是这样计算的：先按基因进行z-score标准化 (x-μ)/SD，然后用下列公式计算两个变量的协方差/标准差乘积。】

取绝对值

（3）对相关系数加权，获得邻接矩阵

为了解决这些问题，提出了加权的思想，WGCNA的做法是对基因表达值之间的相关系数取β次幂。直接结果是把基因间相关性的强弱的差异放大，使强弱关系更为分明（弱化弱连接，强化强连接），有利于后续聚类（模块）识别（使得相关性数值更符合无尺度网络特征，更具有生物意义）。幂次β就是软阈值（power）。对于基因i和j，相关系数为sij，取β次幂后获得aij = sij^β，aij表示网络中基因节点i和j之间的连接强度。因此，最终网络中基因间的相关强度，取决于β次幂的选择，其取值的高低直接影响模块的构建和模块内基因的划分。β值根据接近无标度网络（scale-free network）的最低值确定，pickSoftThreshold这个函数所做的就是确定合适的power，power值的选择即为确定β值的过程。计算所有基因间表达水平的关联强度（加权相关性值），获得邻接矩阵（Adjacency matrix）。

（4）计算拓扑重叠矩阵

为了保证信息的充分利用，进一步计算了拓扑重叠矩阵（Topological overlap matrix，TOM），以降低噪音和假相关，获得的新的相关性矩阵。其思想是，任何两个基因的相关性不仅仅由它们的表达相似性直接决定，并且如果基因i和j有很多相同的邻接基因，意味着两基因间有相似的表达模式，TOM ij就比较高。因此用TOM矩阵用于反映基因间表达的相似性。

（5）模块识别

通过基因之间的加权相关系数（TOM矩阵）构建分层聚类树（平均连接层级聚类），聚类树的不同分支为模式相似的基因，代表不同的基因模块。将没有被分配到特定模块的基因群定义为灰色模块。这样就可以将几万个基因通过基因表达模式被分成了几十个模块，是一个提取归纳信息的过程。模块（Module）是指相关性非常高的基因集。在无向网络中，模块内是高度相关的基因。在有向网络中，模块内是高度正相关的基因。

把表达模式接近的基因聚在一起，就形成了一个个模块。模块内的基因整合成一个eigengene来代表这个模块。

Module eigengene是给定模型的第一主成分，代表整个模型的基因表达谱，用一个向量代替了一个矩阵，方便后期计算。把每个模块里面的基因进行分析，整合得到一个值，用这个值代替整个模块。表达矩阵原本的成千上万行，就变成了每个模块的eigengene的行，矩阵就变的简单了。

Connectivity（连接度）类似于网络中“度”（degree）的概念，每个基因的连接度是与其相连的基因的边属性之和。在共表达网络中，连接度衡量一个基因与所有其他网络基因的相关性。Intramodular connectivity（模块内连接度）衡量给定基因相对于特定模块的基因的连接或共表达程度。

Module Membership(MM) 对于每个基因，我们通过将基因表达谱与某模块的eigengene相关性来定义Module Membership。例如，基因i与蓝色模块eigengene的相关性MM blue(i)，如果接近0则基因i不是蓝色模块的一部分，接近1或-1则它与蓝色模块基因高度正相关或负相关。Gene significance(GS)是基因i与临床性状的相关性，GS越高则基因i的生物学意义越大。关键基因（Hub gene）是连接度最多或连接多个模块的基因。

（6）术语

（7）整体流程

WGCNA旨在寻找协同表达的基因模块，并探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。

整体分析流程为：

1. 构建基因共表达网络：计算加权的表达相关性（TOM矩阵），即前面的1-4步

2. 识别基因集（模块）：平均连接层级聚类，根据设定标准切分聚类结果，获得不同的基因模块，用聚类树的分枝和不同颜色表示

得到模块之后的分析有：

3. 根据TOM矩阵，绘制相关性图

4. 模块的功能富集：利用 WGCNA 算法得到的基因共表达模块，需要进一步的功能注释以发现其生物学意义

5. 计算基因模块与表型的相关性，鉴定性状相关的模块

6. 找到模块的核心基因，选择感兴趣的关键基因

7. 根据模型中已知基因的功能推测未知基因的功能

亲自手写的

三、没有具体流程只有图

（1）样本聚类树和临床trait heatmap

输入的基因表达矩阵，如果是芯片数据，那么常规的归一化矩阵即可；如果是转录组数据，最好是FPKM值或者其它归一化好的表达量，并进行log2(x+1)转换，也可以使用DESeq2的rlog转换的数据。可以进行一些前处理，比如处理缺失值，按表达水平过滤低表达的基因等，避免它们的影响，同时减少后续运算资源消耗（只要能保留感兴趣的重要基因就可以）。

第一次求相关性的时候，其实跟性状是没有关系的；得到模块和eigengene后第二次相关性才跟性状有关。常规的，性状信息应该是数值，比如患者的年龄、生存时间。如果是分类变量，需要使用1或者0来表示。超过二分类的信息，可以使用model.matrix()制造一个设计矩阵添加进去，这时候得到的相关模块代表该亚组特异的模块。除了添加常见的各种临床信息之外，理论上讲，只要是一列数值就可以添加到性状中，比如基因A的表达量，这时候算出来的跟这个“性状”相关的模块就代表和基因A相关的模块。也可以添加甲基化信息、突变信息、SNP信息等，任何能获得的信息都可以添加。

（2）确定软阈值power

上述图形的横轴均代表权重参数β（软阈值），左图纵轴代表对应的网络中log(k)与log(p(k))相关系数的平方，即无标度拓扑拟合指数R^2（scale free topology fitting index R2，即统计结果中的SFT.R.sq列数值）。相关系数的平方越高，说明该网络越逼近无尺度分布，R^2一般要到0.85以上。之所以有负值是因为又乘以了slope列数值的负方向，仅关注正值即可。右图的纵轴代表对应的基因模块中所有基因邻接函数的均值（mean.k.网络的平均连接度）。综合考虑这些指标选择合适的软阈powers值，使R^2尽可能大但也要保证连通度不能太低。可以选择R^2大于0.85或平均连接度小于100时的power值。最佳的beta值就是sft$powerEstimate。如果实际分析中的纵坐标达不到0.85，软阈值可选取经验值。

软阈值的经验值

（3）构建网络

基于基因表达值矩阵计算基因间表达的相似度，获得邻接度矩阵，并进一步计算拓扑重叠矩阵TOM（Pearson相关系数→邻接矩阵软阈值→拓扑重叠矩阵）。TOM矩阵中的值就反映了两两基因间共表达的相似度，取值0-1，基因间共表达相似度越高，值越趋于1。然后推出基因间的相异性系数（1-TOM），并据此得到基因间的系统聚类树。然后按照混合动态剪切树的标准，设置每个基因模块最少的基因数目为30。再次分析，依次计算每个模块的特征向量值，然后对模块进行聚类分析，将距离较近的模块合并为新的模块（基因按照相异性1-TOM聚类→混合动态剪切树识别模块→得到特征向量→模块聚类与合并）。经过这一步，把输入的表达矩阵的几千个基因组归类成了几十个模块。