WGCNA学习：WGCNA分析原理

2020-08-10 本文已影响0人周小钊

本教程根据PlantTech的WGCNA课程编写，课程还是不错的，所以将该课程给大家分享一下。

WGCNA笔记第一弹

1.WGCNA简介

WGCNA(Weighted Gene Co-Expression Network Analysis, 加权基因共表达网络分析)，鉴定表达模式相似的基因集合（module）。解析基因集合与样品表型之间的联系，绘制基因集合中基因之间的调控网络并鉴定关键调控基因。

WGCNA适合于复杂的转录组数据

研究不同器官/组织类型和不同阶段的发育调控、生物和非生物胁迫的不同时间点响应机制

2.主要内容

一、WGCNA原理
- 1.构建基因关系网络
- 2.构建基因模块
- 3.筛选关键基因
- 4.鉴定关键基因
二、WGCNA应用
- 1.应用场景-实验设计、经典文献
- 2.实例分析-单个材料、两个材料
- 3.实例分析-结合表型
三、WGCNA实战

3.WGCNA原理

4.构建基因关系网络

4.1 计算基因间相关关系

基因间相似性（similarity）：根据基因在不同样品中的表达情况，计算任意两个基因间的相关关系。用Pearson相关系数

基因共表达矩阵： S=[S_ij]

S_ij 表示基因i和基因j的Pearson相关系数。

软阈值：通过加权函数将相关系数变换，形成邻接矩阵（Adjacency Matrix）,矩阵中元素连续化。

邻接函数：power函数（幂指数函数）

a_ij=power(S_ij, β)=|S_ij|^β

需要确定邻接函数的参数β，依据无尺度网络原则，即基因表达网络符合无尺度网络的幂函数分布

4.2 无尺度网络

图片.png

网络图的点指图中的每一个节点，度指与该点的连接数

随机网络（Random network），每个节点的度相对平均

无尺度网络（Scale-free network），少数节点具有明显高于一般点的度，这些点被称为hub，由少数hub与其他节点关联，最终构成整个网络

无尺度网络的幂率分布：节点连接数为k的节点数h，k与h成反比，负相关

尺度：随机网络中每个节点的连接数符合泊松分布，大部分节点的连接数居中，中值称为随机网络的尺度。

无尺度网络符合幂率分布，大多数点只有很少的连接，少数点有很多的连接

基因相关关系，幂函数处理后，少数强相关性不受影响或者影响较小，而相关性弱的取n次幂后，相关性明显下降。

4.3 确定关键参数β

寻找合适的β，使得基因表达关系符合无尺度网络，度数高的节点少，度数低的节点多。

节点度数k与具有该度数节点的个数h服从幂律分布

具体计算度数为k的节点个数的对数值log（k），与该节点出现的概率对数（log(p(k))）呈现负相关，一般会设置相关系数大于0.8

为了检测设置的参数β是否满足无尺度网络，对log₁₀(p(k))和log₁₀(k)作图，同时为更好评估，对两者之间的相关系数做平方，即R²。如果模型R²接近1，则两者之间为很好的线性关系。

4.4 计算基因间表达关系

评估基因间表达关系：直接关系

生物体内基因间的关系：直接关系+间接关系

TOM：用拓扑重叠（topological overlap measure，TOM）来计算基因之间关联程度，除了分析两个基因之间的关系，还考虑这两个基因与其他基因之间的连接。这样更具有生物学意义。

建立TOM矩阵 $（TOM）\Omega=[\omega_{ij}]=[TOM_{ij}]$

$TOM_{ij}=\frac{\sum_{\mu}a_{i\mu}a_{\mu j}+a_{ij}}{min(k_i,k_j)+1-a_{ij}}$

TOM公式中，计算i与j之间的关系，不仅考虑了i和j的直接关系，还考虑了第三个基因μ的间接关系

5 构建基因模块

5.1 层次聚类树

基因模块的划分基于基因间的连接稀疏性，将TOM矩阵（Similarity）转化为相异度矩阵（Dissimilarity）: $d^\omega_{ij}=1-\omega_{ij}$

利用基于TOM值的相异度 $d^\omega_{ ij}$ 层次聚类建树

建树方法：动态剪切树和静态剪切树

5.2 动态混合剪切法

第一步：识别满足设定条件的初级模块
- 1.满足模块预定义的最低基因数目
- 2.距离集群过远的基因，即使与集群处于同一分支，也去除
- 3.每个集群与其他周围的集群显著不同
- 4.处在树分支尖端的每个群集的核心基因紧密相连
第二步：测试步骤
- 将未分配的基因进行测试，如果足够接近某个初级群集，则分配进去
- 通常WGCNA使用动态混合剪切法建树

5.3 建树过程的参数

模块最少基因数目（minModuleSize）
合并模块的最小距离（mincutHeight）计算模块的特征值，利用模块特征值建树，合并距离很近的模块（如Height小于0.2）
模块特征值（Epigengene）
模块内所有基因进行主成分分析（PCA），第一主成分的值即为Epigengene。它代表该模块内基因表达的整体水平

6 筛选基因模块

主要介绍四种方法

6.1 表达模式分析

模块表达模式分析：模块在各个样品中的丰度

模块特征值（Epigengene）：模块内所有基因进行主成分分析（PCA），第一主成分的值即为Epigengene。它代表该模块内基因表达的整体水平。

如果某模块在样品中特征值正或负表达较高，说明模块与这个样品关系紧密

6.2 模块与表型性状关联分析

模块显著性值（Module significance，MS）：模块内所有基因的基因显著性值的平均值。

基因显著性值（Gene significance, GS）:基因表达水平与因变量水平的相关系数。用T检验计算每个基因在不同表型样品组间的差异表达显著性检验P值（Pearson相关系数），通常将P值取以10底对数值定义为基因显著性GS

计算各模块与一表型性状的MS值，如一个模块的MS值显著高于其他模块，则这一模块与该性状存在关联关系

模块特征值显著性（Epigengene significance, ES）:模块特征值与某一性状的相关系数，筛选与性状关联度最高的模块

6.3 富集分析

对各个模块都进行GO和KEGG富集分析，找出与我们研究性状相关通路相关性最强的模块进行深入挖掘

6.4 依据目标基因筛选模块

依据研究目的、前期研究结果和已发表文献，有重点关注的目标基因，可直接筛选目标基因所在的基因模块重点进一步分析

7 鉴定关键基因

7.1 模块内部基因连接度分析

Connectivity(degree)-连接度：与某个基因连接的所有其他基因的总和，即描述一个基因与其他所有基因的关联程度，一般用K值表示。

Intramodular connectivity KIM-模块内部连接度IC:某个模块中的基因与该模块中其他基因的关联程度（共表达程度）。可用来衡量模块身份（module membership,MM).

Module Membership MM,or Epigengene-based connectivity KME:模块身份，用一个基因在所有样本中的表达语与某个模块特征值的表达谱的相关性，来衡量这个基因在这个模块中的身份。

KME值接近0,说明这个基因不是该模块的成员：KME接近1或者－1,说明这个基因与该模块密切相关（正相关或者负相关）。
可以对所有基因计算相对某个模块的KME值，并不一定要是该模块的成员。
KME与KIM高度相关。某个模块中KIM值高的hub基因一定与该模块的KME也很高。
KME与KIM的区别：IC衡量基因在特定模块中的身份，MM衡量基因在全局网络中的位置。

筛选关键基因：
TOM值（模块调控系表中的weight值）大于阈值（默认是0.15)的两个基因才认为是相关的，然后计算每个基因的连接度。即先筛选有足够强度的关系，然后计算连接度。

模块内部高连接度的基因，模块内排名前30或者10%(KME或KIM).
筛选关键基因：将该基因模块身份MM相对于基因显著性GS做散点图，选择右上角MM和GS均高的基因进一步分析。
基因显著性值（Gene significance,GS)因变量水平的相关系数。衡量基因与表型性状的关联程度，GS越高，说明与表型越相关，越具有生物学意义。GS可以为正值或负值（正相关或负相关）
Cytoscape中一般用weight值（TOM值）来绘制网络图。

7.2 特定功能基因分析

高连通性的基因一般位于调控网络的上游；低连通性的基因一般位于调控网络的下游。

调控网络上游一般是调控因子，如转录因子；下游一般是功能性的酶或蛋白分子

重点关注具有调控功能的基因，典型的为转录因子，这些基因往往是关键基因

7.3 目标基因关联分析

依据研究目的，选取跟目标基因关系紧密的基因，如筛选与目标基因的TOM值排名前10，或者TOM值大于0.2的基因

可准确筛选与目标基因存在上下游调控关系的候选基因

当目标基因连接度不高时，可筛选与目标基因TOM值很高，且自身连接度也很高的基因

参考文献：

Langfelder, P., Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9, 559 (2008). https://doi.org/10.1186/1471-2105-9-559

转载请注明周小钊的博客>>WGCNA学习：WGCNA分析原理