【WGCNA】WGCNA学习（一）

其实我一直没用过WGCNA，因为分析网络的方法有很多，但是大家好像都更爱用这个。最近帮人分析的几组数据，他们指名要用WGCNA分析，所以就学习一下。

=======WGCNA简介=========

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

 

WGCNA适合于复杂的转录组数据，研究不同器官/组织类型和不同阶段的发育调控、生物和非生物胁迫的不同时间点响应机制。

======WGCNA中的几个概念======

共表达网络：定义为加权基因网络。点代表基因，边代表基因表达相关性。加权是指对相关性值进行冥次运算(冥次的值也就是软阈值 (power,

pickSoftThreshold这个函数所做的就是确定合适的power))。无向网络的边属性计算方式为abs(cor(genex, geney)) ^ power；有向网络的边属性计算方式为(1+cor(genex, geney)/2) ^ power; sign hybrid的边属性计算方式为cor(genex,geney)^power if cor>0 else 0。这种处理方式强化了强相关，弱化了弱相关或负相关，使得相关性数值更符合无标度网络特征，更具有生物意义。如果没有合适的power，一般是由于部分样品与其它样品因为某种原因差别太大导致的，可根据具体问题移除部分样品或查看后面的经验值。

 

Module(模块)：高度內连的基因集。在无向网络中，模块内是高度相关的基因。在有向网络中，模块内是高度正相关的基因。把基因聚类成模块后，可以对每个模块进行三个层次的分析：1. 功能富集分析查看其功能特征是否与研究目的相符；2. 模块与性状进行关联分析，找出与关注性状相关度最高的模块；3. 模块与样本进行关联分析，找到样品特异高表达的模块。

 

Connectivity (连接度)：类似于网络中"度"(degree)的概念。每个基因的连接度是与其相连的基因的边属性之和。

 

Module eigengene E: 给定模型的第一主成分，代表整个模型的基因表达谱。

 

Intramodular connectivity: 给定基因与给定模型内其他基因的关联度，判断基因所属关系。

 

Module membership: 给定基因表达谱与给定模型的eigengene的相关性。

 

Hub gene: 关键基因 (连接度最多或连接多个模块的基因)。

 

Adjacency matrix(邻接矩阵)：基因和基因之间的加权相关性值构成的矩阵。

 

TOM (Topological overlap matrix)：把邻接矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得的新距离矩阵，这个信息可拿来构建网络或绘制TOM图。

基本分析流程如下：

构建基因共表达网络：使用加权的表达相关性。

识别基因集：基于加权相关性，进行层级聚类分析，并根据设定标准切分聚类结果，获得不同的基因模块，用聚类树的分枝和不同颜色表示。

如果有表型信息，计算基因模块与表型的相关性，鉴定性状相关的模块。

研究模型之间的关系，从系统层面查看不同模型的互作网络。

从关键模型中选择感兴趣的驱动基因，或根据模型中已知基因的功能推测未知基因的功能。导出TOM矩阵，绘制相关性图。

1. 构建基因关系网络

1.1 计算基因间相关关系

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

 

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

 

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

 

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

 

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

 

aij=power(Sij, β)=|Sij|β

 

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

1.2 无尺度网络

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

 

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

 

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

 

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

 

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

 

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

 

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

1.3 确定参数β

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

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

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

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

1.4 计算基因间表达关系

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

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

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

 

建立TOM矩阵：

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

2 构建基因模块

2.1 层次聚类树

基因模块的划分基于基因间的连接稀疏性，将TOM矩阵（Similarity）转化为相异度矩阵（Dissimilarity）

利用基于TOM值的相异度

层次聚类建树

 

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

2.2 动态混合剪切法

第一步：识别满足设定条件的初级模块

1.满足模块预定义的最低基因数目

2.距离集群过远的基因，即使与集群处于同一分支，也去除

3.每个集群与其他周围的集群显著不同

4.处在树分支尖端的每个群集的核心基因紧密相连

第二步：测试步骤

将未分配的基因进行测试，如果足够接近某个初级群集，则分配进去

通常WGCNA使用动态混合剪切法建树

2.3 建树过程的参数

模块最少基因数目（minModuleSize）

 

合并模块的最小距离（mincutHeight）计算模块的特征值，利用模块特征值建树，合并距离很近的模块（如Height小于0.2）

 

模块特征值（Epigengene）

 

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

 

3 筛选基因模块

3.1 表达模式分析

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

 

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

 

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

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

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

 

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

 

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

 

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

3.3 富集分析

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

4.4 依据目标基因筛选模块

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

5 鉴定关键基因

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

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

 

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

 

Module Membership MM,or Epigengene-basedconnectivity 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值）来绘制网络图。

5.2 特定功能基因分析

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

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

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

5.3 目标基因关联分析

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

 

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

 

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

===WGCNA安装===

source("https://bioconductor.org/biocLite.R")

biocLite(c("AnnotationDbi","impute","GO.db", "preprocessCore"))

site="https://mirrors.tuna.tsinghua.edu.cn/CRAN"

install.packages(c("WGCNA","stringr", "reshape2"), repos=site)

本文使用 文章同步助手 同步