基于TCGA肿瘤数据进行基因共表达网络分析

WGCNA概述
WGCNA相关概念
WGCNA基本分析流程
WGCNA应用示例
- 输入数据和参数选择
- 安装WGCNA
- WGCNA示例
- - 1. 数据输入
  - 2. 构建加权相关性邻接矩阵
  - 3. 对相关系数邻接矩阵加权优化（软阈值法）
  - 4. 计算拓扑重叠矩阵
参考资料

WGCNA概述

WGCNA（Weighted Gene Coexpression Network Analysis）是一个基于基因表达数据构建基因共表达网络的方法。WGCNA和差异基因分析的区别在于:

差异基因分析主要针对样本之间的差异，
WGCNA主要针对基因之间的关系。

WGCNA原文

WGCNA begins with the understanding that the information captured by microarray experiments is far richer than a list of
differentially expressed genes. Rather, microarray data are more completely represented by considering the relationships
between measured transcripts, which can be assessed by pair-wise correlations between gene expression profiles. In most microarray data analyses, however, these relationships go essentially unexplored. WGCNA starts from the level of
thousands of genes, identifies clinically interesting gene modules, and finally uses intramodular connectivity, gene
significance (e.g. based on the correlation of a gene expression profile with a sample trait) to identify key genes
in the disease pathways for further validation.

WGCNA 从数千个基因的层面开始，识别临床上感兴趣的基因模块，最后使用模块内连接、基因显著性（例如基于基因表达谱与样本特征的相关性）来识别疾病通路中的关键基因，以进一步验证。

WGCNA通过分析基因之间的关联性，将基因区分为多个模块。最后通过这些模块和样本表型之间的关联性分析，寻找特定表型的分子特征。WGCNA根据需要计算的步骤可大致分成以下几个步骤：

数据前处理
构建加权相关性邻接矩阵
计算拓扑重叠矩阵
对基因进行层次聚类，划分模块
总之，WGCNA方法旨在寻找协同表达的基因模块（module），并探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。适用于复杂的数据模式，一般可应用的研究方向有：不同器官或组织类型发育调控，统一组织不同发育调控，非生物胁迫不同时间点应答、病原菌侵染后不同时间点应答。

WGCNA alleviates the multiple testing problem inherent in microarray data analysis. Instead of relating thousands of genes to
a microarray sample trait, it focuses on the relationship between a few (typically less than 10) modules and the sample trait. Toward this end, it calculates the eigengene significance (correlation between sample trait and eigengene) and the corresponding
p-value for each module. The module definition does not make use of a priori defined gene sets. Instead, modules are constructed
from the expression data by using hierarchical clustering. Although it is advisable to relate the resulting modules to gene ontology information to assess their biological plausibility,
it is not required.

WGCNA 缓解了微阵列数据分析中固有的多重测试问题。它不是将数千个基因与微阵列样本特征相关联，而是关注少数（通常少于 10 个）模块与样本特征之间的关系。为此，它计算每个模块的特征基因显著性（样本性状和特征基因之间的相关性）和相应的 p 值。模块定义不使用先验定义的基因集。相反，模块是通过使用层次聚类从表达式数据构建的。尽管建议将生成的模块与基因ontology信息相关联以评估其生物学合理性，但这不是必需的。

加权基因共表达网络（WGCNA）是用来描述不同样品之间基因关联模式的系统生物学方法，可以用来鉴定高度协同变化的基因集，并根据基因集的内连性和基因集与表型之间的关联鉴定候补生物标记基因或治疗靶点。
与只关注差异表达的基因相比，WGCNA利用数千或近万个变化最大的基因或全部基因的信息识别感兴趣的基因集，并与表型进行显著性关联分析。一是充分利用了信息，二是把数千个基因与表型的关联转换为数个基因集与表型的关联，免去了多重假设检验校正的问题。

从方法上将，WGCNA分为表达量聚类分析和表型关联两部分，主要包括基因之间的相关系数计算、基因模块的确定、共表达网络、模块与性状关联四大步骤。

第一步：计算任意两个基因之间的相关系数（Pearson Coefficient）。一般需要设置阈值来进行筛选，高于阈值的则认为是相似的。但是如果，如果将阈值设为0.8，则很难说明0.8与0.79具有显著差别。因此，WGCNA分析时采用相关系数加权值，即对基因相关系数取N次幂，使得网络中基因之间的连接服从无尺度网络分布（scale-free networks），这种算法具有生物学意义。
第二步：通过基因之间的相关系数构建分层聚类书，聚类树的不同分支代表不同的基因模块，不同颜色代表不同的模块。基于基因的加权相关系数，将基因按照表达模式进行分类，将模式相似的基因归为一个模块。这样就可以将众多基因通过基因表达模式分成几十个模块，这也是一个提取归纳信息的过程。

Because the modules may correspond to biological pathways, focusing the analysis on intramodular hub genes (or the module
eigengenes) amounts to a biologically motivated data reduction scheme.
由于模块可能对应于生物通路，因此将分析重点放在模块内hub基因（或模块特征基因）上，这相当于一种生物驱动的数据减少方案。
Because the expression profiles of intramodular hub genes are highly correlated, typically dozens of candidate biomarkers result.
由于模块内hub基因的表达谱高度相关，通常会产生数十种候选生物标志物。
Although these candidates are statistically equivalent, they may differ in terms of biological plausibility or clinical utility.
尽管这些候选基因在统计上是等效的，但它们在生物学合理性或临床效用方面可能有所不同。
Gene ontology information can be useful for further prioritizing intramodular hub genes.
基因ontology信息可用于进一步确定模块内hub基因的优先级。

WGCNA相关概念

共表达网络。定义为加权基因网络。点表示基因，边表示基因表达相关性。加权是指对相关性值进行幂次运算（幂次的值也就是软阈值soft thresholding（power：β≥1\beta \ge 1β≥1, pickSoftThreshold）这个函数所作的就是确定合适的power）。
- 无向网络(unsigned network) 的边属性计算方式为abs(cor(genex, geney))^power：ai,j=∣cor(xi,xj)∣βa_{i,j}=|cor(x_i,x_j)|^\betaai,j=∣cor(xi,xj)∣β
- 有向网络(signed co-expression network) 的边属性计算方式为(1+cor(genex,geney)/2)^power：ai,j=∣(1+cor(xi,xj))/2∣βa_{i,j}=|(1+cor(x_i,x_j))/2|^\betaai,j=∣(1+cor(xi,xj))/2∣β
- sign hybrid的边属性计算方式为cor(genex, geney)^power if cor>0 else 0。

We define co-expression networks as undirected, weighted gene networks. The nodes of such a network correspond to gene expression
profiles, and edges between genes are determined by the pairwise correlations between gene expressions.
我们将共表达网络定义为无向、加权的基因网络。 这种网络的节点对应于基因表达谱，基因之间的边由基因表达之间的成对相关性决定。By raising the absolute value of the correlation to a power(soft thresholding), the weighted gene co-expression network
construction emphasizes high correlations at the expense of low correlations.
通过将相关性的绝对值提高到幂（软阈值），加权基因共表达网络构建以牺牲低相关性为代价强调高相关性。

这种处理方式强化了强相关，弱化了弱相关或负相关，使得相关性数值更符合无标度网络特征，更具有生物意义。如果没有合适的power，一般是由于部分样品与其他样品因为某种原因差别太大导致的，可根据具体问题移除部分样品或查看后面的经验值。

Module（模块）：高度内联的基因集。在无向网络中，模块内是高度相关的基因。在有向网络中，模块内是高度正相关的基因。

Modules are clusters of highly interconnected genes. In an unsigned coexpression network, modules correspond to clusters of genes
with high absolute correlations. In a signed network, modules correspond to positively correlated genes

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

Connectivity（连接度）：类似于网络中”度(degree)“的概念。每个基因的连接都是与其相连的基因的边属性之和。
Module eigengene E：给定模型的第一主成分，代表整个模型的基因表达谱。即用一个向量代替一个矩阵，方便后期计算。
Intramodular connectivity：给定基因与给定模型内其他基因的关联度，判断基因所属关系。
Module membership：给定基因表达谱与给定模型的eigengene的相关性。
Hub gene：关键基因（连接度最多或连接多个模块的基因）。
Adjacency matrix（邻接矩阵）：基因和基因之间的加权相关性值构成的矩阵。
TOM（Topological overlap matrix）：把临界矩阵转换为拓扑重叠矩阵，以降低噪声和假相关，获得的新距离矩阵，这个信息可以拿来构建网络或绘制TOM图。

WGCNA基本分析流程

构建基因共表达网络：使用加权的表达相关性
识别基因集：基于加权相关性，进行层次聚类分析，并根据设定标准切分聚类结果，获得不同的基因模块，用聚类树的分支和不同颜色表示。
如果有表型信息，计算基因模块与表型的相关性，鉴定性状相关的模块。
研究模型之间的关系，从系统层面查看不同模型的互作网络。
从关键模型中选择感兴趣的驱动基因，或根据模型中已知基因的功能来推测未知基因的功能。
导出TOM矩阵，绘制相关性图。

WGCNA应用示例

R包WGCNA是用于计算各种加权关联分析的功能集合，可用于网络构建,基因筛选,基因簇鉴定,拓扑特征计算,数据模拟和可视化等。

输入数据和参数选择

注意事项：

WGCNA本质是基于相关系数的网络分析方法，适用于多样品数据模式，一般要求样本数多于15个。样本数多于20时效果更好，样本越多，结果越稳定。
基因表达矩阵：常规表达矩阵即可，即基因在行，样品在列，进入分析前做一个转置。RPKM、FPKM或其它标准化方法影响不大，推荐使用Deseq2的varianceStabilizingTransformation或log2(x+1)对标准化后的数据做个转换。如果数据来自不同的批次，需要先移除批次效应 。如果数据存在系统偏移，需要做下quantile normalization。
性状矩阵：用于关联分析的性状必须是数值型特征（如下面示例中的Height, Weight,Diameter）。如果是区域或分类变量，需要转换为0-1矩阵的形式(1表示属于此组或有此属性，0表示不属于此组或无此属性，如样品分组信息WT, KO, OE)。

ID  WT  KO  OE Height Weight Diameter
samp1   1   0   0   1   2   3
samp2   1   0   0   2   4   6
samp3   0   1   0   10  20  50
samp4   0   1   0   15  30  80
samp5   0   0   1   NA  9   8
samp6   0   0   1   4   8   7

推荐使用Signed network和Robust correlation (bicor)。(这个根据自己的需要选择)
无向网络在power小于15或有向网络power小于30内，没有一个power值可以使无标度网络图谱结构R×2R\times 2R×2达到0.8或者平均连接度降到100以下，可能是由于部分样品与其他样品差别太大造成的。这可能由批次效应、样品异质性或实验条件对表达影响太大等造成，可以通过绘制样品聚类查看分组信息、关联批次信息、处理信息和有无异常样品 (可以使用热图简化，增加行或列属性)。如果这确实是由有意义的生物变化引起的，也可以使用后面程序中的经验power值。

安装WGCNA

if (!requireNamespace("BiocManager", quietly = TRUE))install.packages("BiocManager")
BiocManager::install(c("WGCNA","stringr","reshape2"))

WGCNA示例

1. 数据输入

WGCNA需要的目标数据是储存有若干样本以及其对应的各基因表达量的数据集，经过数据前处理将其转化为合适的数据结构，以进行后续的分析。
我们需要的样本为由m个样本、n个基因构成的m x n矩阵，矩阵内的元素为某样本中某基因的基因表达量值。获取数据后先检查其数据是否符合格式，例如如果样本和基因行列对调了，将其矩阵转置即可。
如果不确定样本数据是否可靠，例如是否有异常值，可以对样本数据进行一次聚类分析，若有某样本有异常值，则其在聚类分析后可视化显离群现象，可考虑将其剔除。

> library(WGCNA)
载入需要的程辑包：dynamicTreeCut
载入需要的程辑包：fastcluster载入程辑包：‘fastcluster’The following object is masked from ‘package:stats’:hclust载入程辑包：‘WGCNA’The following object is masked from ‘package:stats’:corWarning message:
程辑包‘WGCNA’是用R版本4.1.2 来建造的
> library(reshape2)
> library(stringr)
> options(stringsAsFactors = FALSE)
> enableWGCNAThreads() # 打开多线程
Allowing parallel execution with up to 15 working processes.

2. 构建加权相关性邻接矩阵

计算基因表达量的相关系数
对于m×nm\times nm×n的表达矩阵，取每两两基因在所有样本中的表达量（即取任意两列），得两个列向量。由于基因表达量数据参差不齐，有时偏差较大，因此对列向量进行以下方法得标准化，以便后续计算：
scale(x)u=xu−mean(x)var(x)scale(x)_u=\frac{x_u-mean(x)}{\sqrt{var(x)}}scale(x)u=var(x)xu−mean(x)
其中xux_uxu为原向量，scale(x)uscale(x)_uscale(x)u为标准化后得向量。

对两个标准化后的列向量以Pearson相关进行线性相关性计算。Pearson相关系数适用于两组连续变量，并且：
a) 两个变量分别服从正态分布，通常用t检验检查相关系数的显著性；
b) 两个变量的标准差不为0。
Pearson相关系数（连续变量）的定义为：两个连续变量(X,Y)的pearson相关性系数P(x,y)等于它们之间的协方差cov(X,Y)除以它们各自标准差的乘积(σX,σY)。系数的取值总是在-1.0到1.0之间，接近0的变量被成为无相关性，接近1或者-1被称为具有强相关性。

3. 对相关系数邻接矩阵加权优化（软阈值法）

先介绍两个概念：随机网络和无标度网络。在数学的图论中，这里提到的网络其实就是图的概念。

随机网络的概念是：在网络中所有点之间的关联性都差不多，可以看做是随机连接而形成的。
无标度网络的概念是：在网络中有极少数的中心点与网络里大部分的点连接，而大部分的点只与那些极少数的中心点连接。
度的概念：度可以理解为网络中某个结点相连的线段数量。

事实上，无标度网络更符合我们认知的事实：极少数的基因执行着生物体中更多的功能，而大多数的外围基因执行着较少的，较为边缘的功能。在受到外因伤害时，无标度网络下的基因会有更强的抵抗力，只要中心基因不受影响，伤害可以忽略不计；而随机网络下的基因在受到外因伤害时，其受损程度和伤害程度完全成正比。因此这也是符合生物进化论中适者生存的理论的。

无标度网络的定义：假设在网络中任取一个结点，取到一个度数为ddd的结点的概率为P(d)P(d)P(d)，且P(d)P(d)P(d)与ddd的概率分布满足幂律分布：
P(d)≈d−αP(d)\approx d^{-\alpha}P(d)≈d−α
由于ddd和P(d)P(d)P(d)取值范围均为0到正无穷，因此可以取对数变换得：
logP(d)≈−αlogdlogP(d)\approx -\alpha log dlogP(d)≈−αlogd
可以看出，如果网络完全符合无标度网络的特征，则其概率密度函数取对数后的logP(d)logP(d)logP(d)与logdlog dlogd线性相关。基于这个原理，我们对各个结点的度数 d 及该度数在所有结点度数中的占比进行pearson线性相关分析，得到关于无标度网络的适应系数（scale-free fitting index R2）。这个系数可以评定构建的网络接近无标度网络的程度，越接近1则越像无标度网络，越接近0则越像随机网络。

基于以上理论，就产生了给基因的相关系数加上一个指数 β 的想法。当β＞1β＞1β＞1时，不同的相关系数之间的差距会越来越大。例如基因A和基因C的相关系数为0.3，基因A和B为0.6，当给他们加上β=2的次方时，相关系数变为0.09和0.36，放大了其相关系数差异的显著程度。这个β可以根据前文提到的无标度网络适应系数进行选择，例如我们要构建无标度网络适应系数大于0.9的网络，就计算各个β值加权后网络的无标度网络适应系数，超过0.9的最小β是最优选择。
Aij=(∣cor(xi,xj)∣)βA_{ij}=(|cor(x_i, x_j)|)^\betaAij=(∣cor(xi,xj)∣)β
在R语言的WGCNA包中pickSoftThreshold函数已将此功能实现并可视化，可供手动选择β的值。
至此，加权的基因相关性矩阵构建完毕。

4. 计算拓扑重叠矩阵

如果直接以上一步得到的两两基因间的加权相关系数进行模块识别和划分准确性较差，且对数据的利用度也不高。因此引入拓扑重叠矩阵的概念。
这个概念可以通过一个社交网络的例子解释：假设A、B同属于同一个公司，且A与B不认识，互相无关联，但A与B都各有该公司的很多共同好友，则可以认为A与B的社交网络重叠性强，因此认为A与B在同一模块（所属的公司）。
在WGCNA中应用类似：如果基因 i 与基因 j 之间可以通过若干中间关联的基因 u 而形成联系，则这个联系的强弱会一定程度上被考虑进化分模块的依据中。在WGCNA中，拓扑重叠矩阵的计算方法为：
TOMij=∑u≠i,j(ai,u∗au,j)−aijmin{∑u≠iaiu,∑u≠jauj}+1−aijTOM_{ij}=\frac{\sum_{u\ne i,j} (a_{i,u} * a_{u,j})-a_{ij}}{min\{\sum_{u\ne i} a_{iu}, \sum_{u\ne j} a_{uj}\}+1-a_{ij}}TOMij=min{∑u=iaiu,∑u=jauj}+1−aij∑u=i,j(ai,u∗au,j)−aij
这里可以算作是WGCNA的精髓了。以 u 为遍历因子，遍历基因列表里除$ i，j 以外的所有基因，并将u与以外的所有基因，并将 u 与以外的所有基因，并将u与i，j$之间的相关系数以上述公式运算，巧妙地用其他基因作为桥梁将 i,ji, ji,j 基因在网络中表现出的重叠程度数据化。如果i,ji, ji,j在网络中连接情况越相似，则TOMijTOM_{ij}TOMij越大，i和j越有可能在同一个表达模块中。
至此，拓扑重叠矩阵计算完毕。

分析流程如下：

7. 从TCGA数据库下载要研究的癌症的基因、miRNA表达数据
8. 分别对基因、miRNA进行表达差异分析
9. 分别针对基因和miRNA进行WGCNA分析，获得聚类模块。并分析与病理学分期相关联的基因模块
10. 验证基因模块与临床特征相关的基因
11. 选择模块中的hub gene和hub miRNA进行生存分析
12. 从病理分期相关联的模块中验证的基因和miRNA做代谢通路分析，构建miRNA-基因的调控网络

下面以TCGA乳腺癌基因数据为例，运用WGCNA对乳腺癌各个亚型做共表达模块。我们首先从TCGA数据库里下载乳腺癌的RNA-seq数据和临床病理资料。

参考资料

[1] WGCNA: an R package for weighted correlation network analysis
[2] WGCNA分析，简单全面的最新教程
[3] Pearson，Spearman，Kendall三大相关系数简单介绍
[4] WGCNA总结笔记：深入研究算法及其实现原理（1）

【生信分析】基于TCGA肿瘤数据进行基因共表达网络分析相关推荐

基因共表达网络分析java,RNA-seq数据的基因共表达网络分析
目录 1 背景 BackgroundTypes of biological networks Motivation for using co-expression networks Network i ...
肿瘤/非肿瘤/单基因/单细胞/非编码：史上最全生信分析攻略！！！
解读生信之美,探讨每篇文献背后的逻辑非肿瘤专栏:条条大路通罗马 1.4+非肿瘤生信分析+铁死亡/焦亡/自噬/代谢/免疫的万能钥匙短评:适合一些热门机制如铁死亡/焦亡/自噬等在非肿瘤疾病中的研究 2 ...
这是入门生信，学习生信分析思路和数据可视化的首选？
封面来源:https://www.zhihu.com/question/304747766 常规转录组是我们最常接触到的一种高通量测序数据类型,其实验方法成熟,花费较低,是大部分CNS必备的技术,以后 ...
图形化开放式生信分析系统开发 - 9 Illumina测序仪测序数据自动拆分
前文链接: 图形化开放式生信分析系统开发 - 1 需求分析及技术实现图形化开放式生信分析系统开发 - 2 样本信息处理图形化开放式生信分析系统开发 - 3 生信分析流程的进化图形化开放式生信分析 ...
免费生信课程|多组学数据整合分析之转录组和蛋白质组分析
搜索"基因组Genome",轻松关注不迷路生科云网址:https://www.bioincloud.tech/ 01 课程简介多组学技术是结合两种或两种以上组学研究方法,如基因 ...
iMeta | Sangerbox: 交互式整合临床生信分析平台
点击蓝字关注我们 Sangerbox:交互.友好的整合临床生信分析平台 https://doi.org/10.1002/imt2.36 ●2022年7月8日,哈尔滨医科大学宋湘团队联合杭州慕谷科技有 ...
生信分析矫正P值_TCGA数据分析系列（二）：数据库之GEPIA2
所谓工欲善其事,必先利其器,从今天开始,我们来介绍TCGA数据库的使用.今天我们来介绍一款非常容易上手的数据库:GEPIA2 链接点这里话不多说,直接进入正题 GEP1A2数据来源如图所示,GEP ...
生信分析流程构建的几大流派
导言构建生信分析流程是生物信息学从业人员必备的技能之一,对该项能力的评估常常是各大公司招录人员的参考项目之一. 在进行 ngsjs 项目时,我做了一张示意图来表示一些高通量测序数据分析项目重现性的要 ...
福建农林大学朱方捷组招聘讲师/副教授/助理——生信分析方向
福建农林大学转录系统生物学课题组招聘(生信分析.组培转化) 工作地点: 福建农林大学海峡联合研究院薪金: 18-30万招聘岗位: 讲师/副教授/助理实验室网址: http://hbmcsysbi ...
找不到r低版本_R 语言与数据挖掘直播班开始招生，生信分析帮你发高分文章
数据单薄很难支撑文章内容?数据有了不知道怎么处理作出高级的图片?这个时候需要的是生信分析--深度的数据挖掘和分析处理,可以帮助临床医生通过数据处理得到自己想要的信息,更快速地发文章. 学习哪种生信分析 ...

【生信分析】基于TCGA肿瘤数据进行基因共表达网络分析