TCGA 数据分析实战 —— WGCNA

加权基因共表达网络分析（ WGCNA , Weighted gene co-expression network analysis ）是一种用来描述不同基因在样本中的表达关联模式的系统生物学方法。

通过将表达高度相关的基因聚集成不同的模块，并探究不同模块与样本表型之间的关联。还可以探究模块内的关键基因的功能，作为潜在的生物标志物或治疗靶点进行后续分析

WGCNA 模块识别算法大致包含以下几个步骤：

输入数据的格式要符合行为样本，列为基因的矩阵格式，因为计算的是基因之间的相关性，所以数据可以是标准化的表达值或者是 read counts 。

探针集或基因可以通过平均表达量或方差（如中位数或绝对中位差）进行过滤，因为低表达或无变化的基因通常代表噪音。

注意 ：并不推荐使用差异基因作为输入矩阵，通过差异表达基因过滤将会导致一个（或几个高度相关的）基因聚成一个模块，同时，也破坏了无标度拓扑的假设，所以通过无标度拓扑拟合来选择软阈值的将会失败。

主要是过滤一些离群或异常的样本，可以对样本数据进行聚类，如果存在异常样本，则其在聚类图中会显示出离群现象，可考虑将其剔除。

首先，对基因的表达量进行 0-1 标准化，即

其中， 为样本方差

然后，使用 pearson 计算基因之间的相关性

两个基因的共表达相似性表示为

然后将基因之间的相似度转换为邻接值，对于非加权网络，计算方式为

其中 为硬阈值，大于等于该阈值表示这两个基因之间存在连接，而低于阈值则认为两个基因没有连接。它们并不能反映共表达信息的连续性质，因此可能导致信息损失。例如，阈值为 08 ，那 079 是不是应该也有一定的相关性呢？

在介绍软阈值之前，我们先引出两个图论的概念：

度表示为节点所连接的边的数量

无标度网络具有很好的鲁棒性，网络中某些节点的错误并不会导致整个网络的瘫痪，具有很多的代偿连接。而这一特点，与生物体中的复杂生化网络非常类似，只有少数的基因执行着关键性的功能，而大多数的基因执行较为单一的功能。

无标度网络中，节点 d 的度为 k 的概率满足幂律分布

通过对数变换，变为

从这个公式可以看出，节点的度数与其出现的概率是负相关的，通过计算各个节点的度数 k 与该度数 k 在所有节点度数中的占比的 pearson 相关性，我们可以得到关于无标度网络的适应系数。该系数越接近 1 则越像无标度网络，越接近 0 则越像随机网络。

所以，对于加权网络，其邻接值的计算方式为：

当软阈值 时，会让相关系数小的更小，而大的更大。

可以根据适应系数来筛选软阈值

光有邻接矩阵是不够的，基因间的相似性应该要同时体现在其表达和网络拓扑水平，为了能能够尽可能地最小化噪音和假阳性的影响，因此引入了拓扑重叠矩阵

这个概念的主要表达的是，两个基因 a 和 b 之间的相关性，不光考虑两个基因的表达相关性，还需要考虑一些 A 和 B 共有的表达相关基因 u ，如果 u 足够多，则说明 A 与 B 的网络重叠性强，应该被聚成一类

换个说法，两个人之间的亲密度不仅与他们两人之间有关，还与他们的共同好友有关，共同好友越多，说明他们两人之间应该越亲密

计算公式为：

其中， 分别为 i 和 j 的度数

表示的是两个基因的相似性，转换成距离度量就是 ，并使用该值来进行聚类，并分割模块

我们以 TCGA 的乳腺癌数据作为示例，来完整的做一遍 WGCNA 分析

先安装模块

获取 50 个样本的 FPKM 数据， WGCNA 最少需要 15 个样本， 20 个以上的样本会更好，样本越多越好，这里为了方便，我们只挑了 50 个样本

过滤基因，取绝对中位差 top 5000 的基因

过滤异常样本

确定软阈值的时候，需要选择网络类型，不同的网络类型，其计算邻接值的方法是不一样的。

默认为 unsigned

我在 RStudio 中使用 enableWGCNAThreads() 会引发下面的错误

所以，我改用了 allowWGCNAThreads() ，就可以运行了

绘制软阈值曲线

其中横坐标为软阈值的梯度，第一幅图的纵坐标为无标度网络适应系数，越大越好；第二幅图的纵坐标为节点的平均连通度，越小越好。

查看系统给我们推荐的软阈值

与我们从图上看到的结果是一致的，如果出现了异常的值，也就是说在有效的 power 梯度范围内（无向网络在 power 小于 15 ，有向网络 power 小于 30 ），无法使适应系数的值超过 08 ，且平均连接度在 100 以上

可能是由于部分样品与其他样品差别较大。这可能是由于批次效应、样品异质性或实验条件对表达影响太大等因素造成的。

可以对样本绘制聚类图来查看有无异常样品，如果这确实是由于生物学差异引起的，也可以使用下面的经验 power 值。

一步法构建网络，我们使用上面推荐的软阈值 5

查看各模块的基因数量

可以使用 labels2colors 函数将数值转换为颜色名称

使用 plotDendroAndColors 函数来展示各个模块的层次聚类结果

其中，无法聚类到模块中的基因会标示为灰色，如果灰色区域较多，可能由于样本中基因共表达趋势不明显，可能需要调整基因过滤的方法。

展示模块之间的相关性

展示 TOM 矩阵，为了节省时间，我们只使用第一个聚类分支

或者更换一种配色

颜色越深表示基因表达的相关性更高，我们可以看到，模块内的基因之间具有较高的共表达，而模块之间的表达相关性较低

将整个网络全部导出成 Cytoscape 输入文件

保存网络

也可以提取某一模块的基因

获取到基因之后，可以进行富集分析找到相关的生物学通路

我们可以分析各网络模块与样本表型之间的关系，从而找到与我们感兴趣表型相关的模块。

样本表型可以是各种指标，比如肿瘤分期分级、已知的分类亚型、药物响应等，并计算模块与这些表型之间是否具有显著相关性

但是模块是一个矩阵，无法直接计算矩阵和向量之间的相关性，需要转换为向量之间的相关性。

而 WGCNA 选择使用 PCA 的方法对数据降维，并将第一主成分定义为 eigengenes ，然后计算 eigengenes 与表型之间的相关性

先获取并处理临床数据

计算模块与 ER 状态的相关性

如果使用的是其他相关性方法，则可以使用 bicorAndPvalue 函数来计算显著性

绘制相关性图

可以看到有些模块的相关性挺高的，而且也具有显著性。我们计算出模块与表型之间相关性之后，可以挑选最相关的那些模块来进行后续分析。但是，模块本身可能还包含很多的基因，还需要进一步识别关键基因基因。

如何寻找关键基因呢？我们可以计算所有基因与模块之间的相关性，也可以计算基因与表型之间的相关性。如果存在一些基因，既与表型显著相关又跟某个模块显著相关，那么这些基因可能就是非常重要的关键基因了

从上图中，我们可以看到 paleturquoise 具有较高的相关性，且具有显著性，我们就来尝试找找这个模块的关键基因

计算基因与模块的相关性

再计算基因与表型的相关性

展示模块内基因与模块和表型之间的相关性

从图中我们可以看出，基因与表型的相关性和基因与模块的相关性还是有一定的线性趋势的，这说明与表型高度相关的基因，通常也是该表型对应模块内比较重要的基因。

因此，当我们要选择关键基因时，推荐选取散点图中右上角部分的基因，即两个相关性均较大的基因

我们可以导出这个模块的网络

转自“ 医学统计园 ”微信公众号。

读入clinicaljson文件

计算文件长度n，在这里n为348

初始化变量

利用一个for循环由json文件中提取信息

将提取的信息做成一个dataFrame

需要。

当我们使用TCGA表达数据，数据以gene-samplematrix形式整理好了，ICGC也是一个很大的癌症数据库，存放了几十种癌症的数据。我们需要将系数矩阵修改成gene-sample。

DataTable dt = new DataTable();

dtColumnsAdd(new DataColumn("PreRevDate0", typeof(decimal)));

DataColumn col = new DataColumn();

colColumnName = "PreRevDate1";

colExpression = "ABS(ConvertToInt32(PreRevDate0))";

colDataType = typeof(decimal);

dtColumnsAdd(col);

DataRow dr = dtNewRow();

dr["PreRevDate0"] = -1;

dtRowsAdd(dr);

答案就在 TCGA barcode ，样本标签描述了样本类型，是正常的还是异常的。还是对照组。比如胶质瘤RNAseq的bar code ，有174个样本类似于这个：

TCGA-06-0681-11A-41R-A36H-07

TCGA-06-0649-01B-01R-1849-01

第四个字段：11A和01B描述的就是样本类型，1-9是肿瘤，10-19是正常，20-29是对照。A 和 B 我也不知道啥意思。由于TCGA barcode 字段宽度是严格的。因此用substr就可提取

names=colnames(RNAseq_dat)

a=asnumeric(substr(names,14,15))

table(a)

可以看见数据中有5个是正常组织样本

----------------------

Xena 网站（网页链接)有整理好的TCGA数据，包括数据集和样本表格。样本表格数据详细，包含生存期，肿瘤分期分级，突变，亚型等等。

以上就是关于TCGA 数据分析实战 —— WGCNA全部的内容，包括:TCGA 数据分析实战 —— WGCNA、R语言提取TCGA数据库clinical.json中的临床信息、icgc数据库的数据需要log fpkm+1嘛等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！