WGCNA 分析 生信技能树
史庭昕
2023.8.21
WGCNA 分析
(readme.assets/image-20190928155237950.png)
大家首先可以看到3个教程：
- 2016-WGCNA-HCC-hub-gene.pdf 中文文章范例）
- WGCNA_GBMTutorialHorvath.pdf
- WGCNA_YeastTutorialHorvath.pdf
基本概念
WGCNA其译为加权基因共表达网络分析。该分析方法旨在寻找协同表达的基因模块(module)，并探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。
适用于复杂的数据模式，推荐5组(或者15个样品)以上的数据。一般可应用的研究方向有：不同器官或组织类型发育调控、同一组织不同发育调控、非生物胁迫不同时间点应答、病原菌侵染后不同时间点应答。
基本原理
从方法上来讲，WGCNA分为表达量聚类分析和表型关联两部分，主要包括基因之间相关系数计算、基因模块的确定、共表达网络、模块与性状关联四个步骤。
第一步计算任意两个基因之间的相关系数（Person Coefficient）。为了衡量两个基因是否具有相似表达模式，一般需要设置阈值来筛选，高于阈值的则认为是相似的。但是这样如果将阈值设为0.8，那么很难说明0.8和0.79两个是有显著差别的。因此，WGCNA分析时采用相关系数加权值，即对基因相关系数取N次幂，使得网络中的基因之间的连接服从无尺度网络分布(scale-freenetworks)，这种算法更具生物学意义。
第二步通过基因之间的相关系数构建分层聚类树，聚类树的不同分支代表不同的基因模块，不同颜色代表不同的模块。基于基因的加权相关系数，将基因按照表达模式进行分类，将模式相似的基因归为一个模块。这样就可以将几万个基因通过基因表达模式被分成了几十个模块，是一个提取归纳信息的过程。
WGCNA术语
权重(weghted)
基因之间不仅仅是相关与否，还记录着它们的相关性数值，数值就是基因之间的联系的权重(相关性)。
权重
Module
模块(module)：表达模式相似的基因分为一类，这样的一类基因成为模块；
Eigengene
Eigengene（eigen + gene）：基因和样本构成的矩阵，https://en.wiktionary.org/wiki/eigengene
Adjacency Matrix
邻近矩阵：是图的一种存储形式，用一个一维数组存放图中所有顶点数据；用一个二维数组存放顶点间关系（边或弧）的数据，这个二维数组称为邻接矩阵；在WGCNA分析里面指的是基因与基因之间的相关性系数矩阵。 如果用了阈值来判断基因相关与否，那么这个邻近矩阵就是0/1矩阵，只记录基因相关与否。但是WGCNA没有用阈值来卡基因的相关性，而是记录了所有基因之间的相关性。
Topological Overlap Matrix (TOM)
WGNA认为基因之间的简单的相关性不足以计算共表达，所以它利用上面的邻近矩阵，又计算了一个新的邻近矩阵。一般来说，TOM就是WGCNA分析的最终结果，后续的只是对TOM的下游注释。
下游分析
得到模块之后的分析有：
1.模块的功能富集
2.模块与性状之间的相关性
3.模块与样本间的相关系数
挖掘模块的关键信息：
1.找到模块的核心基因
2.利用关系预测基因功能
代码示例
其中第一步数据准备反而是最复杂的，取决于大家的R语言水平，这个数据GSE48213-wgcna-input.RData我已经保存下来咯，如果大家不会做，又想体验一下这个WGCNA流程，就可以直接load我保存好的数据文件即可。
step1: 输入数据的准备
这里主要是表达矩阵， 如果是芯片数据，那么常规的归一化矩阵即可，如果是转录组数据，最好是RPKM/TPM值或者其它归一化好的表达量。然后就是临床信息或者其它表型，总之就是样本的属性。
为了保证后续脚本的统一性，表达矩阵统一用datExpr标识，临床 信息统一用datTraits标识。(PS: 如果你R语言很差，变量名不要轻易修改)
rm(list = ls())
#options(stringsAsFactors = F)
# 切换工作目录
#setwd('/data3/txshi/chip_seq/peak_analysis/WGCNA/test/my_WGCNA-master/230821Rmd/')
load('/data3/txshi/chip_seq/peak_analysis/WGCNA/test/my_WGCNA-master/GSE48213-wgcna-input.RData')#查看数据结构
fpkm[1:4,1:4]## GSM1172844 GSM1172845 GSM1172846 GSM1172847
## ENSG00000000003 95.21255 95.69868 19.99467 65.6863763
## ENSG00000000005 0.00000 0.00000 0.00000 0.1492021
## ENSG00000000419 453.20831 243.64804 142.05818 200.4131493
## ENSG00000000457 18.10439 26.56661 16.12776 12.0873135head(datTraits)## gsm cellline subtype
## GSM1172844 GSM1172844 184A1 Non-malignant
## GSM1172845 GSM1172845 184B5 Non-malignant
## GSM1172846 GSM1172846 21MT1 Basal
## GSM1172847 GSM1172847 21MT2 Basal
## GSM1172848 GSM1172848 21NT Basal
## GSM1172849 GSM1172849 21PT Basaltable(datTraits$subtype)##
## Basal Claudin-low Luminal Non-malignant unknown
## 14 6 27 5 4library(WGCNA)## Loading required package: dynamicTreeCut## Loading required package: fastcluster##
## Attaching package: 'fastcluster'## The following object is masked from 'package:stats':
##
## hclust## ##
## Attaching package: 'WGCNA'## The following object is masked from 'package:stats':
##
## corRNAseq_voom <- fpkm
## 因为WGCNA针对的是基因进行聚类，而一般我们的聚类是针对样本用hclust即可，所以这个时候需要转置。
WGCNA_matrix = t(RNAseq_voom[order(apply(RNAseq_voom,1,mad), decreasing = T)[1:5000],])
datExpr0 <- WGCNA_matrix ## top 5000 mad genes
datExpr <- datExpr0
## 下面主要是为了防止临床表型与样本名字对不上
sampleNames = rownames(datExpr);
traitRows = match(sampleNames, datTraits$gsm)
rownames(datTraits) = datTraits[traitRows, 1] 上面代码里面的rpkm就是我们的转录组数据的表达矩阵，以rpkm为单位。而datTraits就是所有样本对应的表型信息。需要自己制作，这个是学习WGCNA的基础，本次实例代码都是基于这两个数据。至于如何做出上面代码的两个例子，取决于大家自己的项目，我这里给出自己的代码，仅供参考哈！
# 56 breast cancer cell lines were profiled to identify patterns of gene expression associated with subtype and response to therapeutic compounds.
if(F){
## https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48213
#wget -c ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE48nnn/GSE48213/suppl/GSE48213_RAW.tar
#tar -xf GSE48213_RAW.tar
#gzip -d *.gz
## 首先在GSE48213_RAW目录里面生成tmp.txt文件，使用shell脚本
# awk '{print FILENAME"\t"$0}' GSM*.txt |grep -v EnsEMBL_Gene_ID >tmp.txt
# 其实也可以直接使用R来读取GSE48213_RAW.tar里面的gz文件，这里就不演示了
# 可以参考：https://mp.weixin.qq.com/s/OLc9QmfN0YcT548VAYgOPA 里面的教程
## 然后把tmp.txt导入R语言里面用reshape2处理即可
# 这个 tmp.txt 文件应该是100M左右大小哦。
a=read.table('GSE48213_RAW/tmp.txt',sep = '\t',stringsAsFactors = F)
library(reshape2)
fpkm <- dcast(a,formula = V2~V1)
rownames(fpkm)=fpkm[,1]
fpkm=fpkm[,-1]
colnames(fpkm)=sapply(colnames(fpkm),function(x) strsplit(x,"_")[[1]][1])
library(GEOquery)
a=getGEO('GSE48213')
metadata=pData(a[[1]])[,c(2,10,12)]
datTraits = data.frame(gsm=metadata[,1],
cellline=trimws(sapply(as.character(metadata$characteristics_ch1),function(x) strsplit(x,":")[[1]][2])),
subtype=trimws(sapply(as.character(metadata$characteristics_ch1.2),function(x) strsplit(x,":")[[1]][2]))
)
save(fpkm,datTraits,file = 'GSE48213-wgcna-input.RData')
}很明显，这个数据GSE48213-wgcna-input.RData我已经保存下来咯，如果大家不会做，又想体验一下这个WGCNA流程，那么可以找我email求取这个数据哦。我的邮箱是jmzeng1314@163.com
我给大家演示的示例数据大概是下面这个样子：
> head(datTraits) ## 56 个细胞系的分类信息，表型
gsm cellline subtype
GSM1172844 GSM1172844 184A1 Non-malignant
GSM1172845 GSM1172845 184B5 Non-malignant
GSM1172846 GSM1172846 21MT1 Basal
GSM1172847 GSM1172847 21MT2 Basal
GSM1172848 GSM1172848 21NT Basal
GSM1172849 GSM1172849 21PT Basal
> fpkm[1:4,1:4] ## 56个细胞系的36953个基因的表达矩阵
GSM1172844 GSM1172845 GSM1172846 GSM1172847
ENSG00000000003 95.21255 95.69868 19.99467 65.6863763
ENSG00000000005 0.00000 0.00000 0.00000 0.1492021
ENSG00000000419 453.20831 243.64804 142.05818 200.4131493
ENSG00000000457 18.10439 26.56661 16.12776 12.0873135
> 这个数据集里面的56种细胞系被分成了5组，如果要分开两两做差异分析，有10种组合，也就是说需要做10次差异分析，每个差异分析结果都需要去注释，会比较麻烦，这个时候WGCNA就派上用场啦。当然，如果你一定要去做差异分析，我也给你代码：https://github.com/jmzeng1314/my-R/blob/master/10-RNA-seq-3-groups/hisat2_mm10_htseq.R
实际上多个分组，差异分析策略是非常个性化的， 比如：https://mp.weixin.qq.com/s/hc6JkKxyelc7b1M1MRiHRQ
step2:确定最佳beta值
选择合适“软阀值（soft thresholding power）”beta，同样的，也是使用教程标准代码即可：
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)## pickSoftThreshold: will use block size 5000.
## pickSoftThreshold: calculating connectivity for given powers...
## ..working on genes 1 through 5000 of 5000## Warning: executing %dopar% sequentially: no parallel backend registered## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1 0.0944 -0.904 0.885 1040.000 1.02e+03 1810.00
## 2 2 0.4910 -1.580 0.952 328.000 3.03e+02 866.00
## 3 3 0.7030 -1.860 0.983 128.000 1.08e+02 474.00
## 4 4 0.7920 -2.000 0.991 57.300 4.38e+01 283.00
## 5 5 0.8490 -2.060 0.996 28.400 1.95e+01 179.00
## 6 6 0.8810 -2.090 0.991 15.200 9.45e+00 118.00
## 7 7 0.9040 -2.070 0.990 8.640 4.89e+00 80.60
## 8 8 0.9220 -2.040 0.994 5.170 2.67e+00 56.40
## 9 9 0.9330 -2.030 0.995 3.240 1.54e+00 40.50
## 10 10 0.9350 -2.020 0.989 2.100 9.29e-01 30.00
## 11 12 0.9250 -2.030 0.977 0.971 3.63e-01 17.30
## 12 14 0.9210 -2.020 0.982 0.496 1.56e-01 10.50
## 13 16 0.9250 -1.970 0.992 0.275 7.04e-02 6.61
## 14 18 0.8940 -1.960 0.973 0.163 3.31e-02 4.31
## 15 20 0.9220 -1.820 0.986 0.102 1.63e-02 2.89#设置网络构建参数选择范围，计算无尺度分布拓扑矩阵
# Plot the results:
##sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
关键就是理解pickSoftThreshold函数及其返回的对象，最佳的beta值就是sft$powerEstimate
最佳beta值
参数beta取值默认是1到30，上述图形的横轴均代表权重参数β，左图纵轴代表对应的网络中log(k)与log(p(k))相关系数的平方。相关系数的平方越高，说明该网络越逼近无网路尺度的分布。右图的纵轴代表对应的基因模块中所有基因邻接函数的均值。最佳的beta值就是sft$powerEstimate，已经被保存到变量了，不需要知道具体是什么，后面的代码都用这个即可，在本例子里面是6。
即使你不理解它，也可以使用代码拿到合适“软阀值（soft thresholding power）”beta进行后续分析。
step3：一步法构建共表达矩阵
有了表达矩阵和估计好的最佳beta值，就可以直接构建共表达矩阵了。 datExpr： 这是输入数据矩阵，包含基因表达数据。 power：该参数用于指定构建网络时使用的软阈值功率。本代码中将其设置为 sft$powerEstimate，该值在上一步中计算得出。 maxBlockSize： 该参数用于指定每个区块中包含的最大基因数。本代码将其设置为 6000。 TOMType： 该参数用于指定拓扑重叠矩阵的类型。在本代码中设置为 “无符号”。 minModuleSize（最小模块大小）： 此参数用于指定模块的最小尺寸。在本代码中设置为 30。 reassignThreshold（重新分配阈值）： 该参数指定将基因重新分配给模块的阈值。在本代码中设置为 0。 mergeCutHeight（合并切割高度）： 该参数用于指定合并模块的高度。本代码将其设置为 0.25。 numericLabels： 此参数用于指定模块标签是否应为数字。本代码将其设置为 “true”。 pamRespectsDendro（尊重树枝图）：该参数用于指定 PAM 聚类是否应尊重树枝图。在本代码中设置为 FALSE。 saveTOMs：保存拓扑图： 该参数用于指定是否保存拓扑重叠矩阵。在本代码中设置为 FALSE。
net = blockwiseModules(
datExpr,
power = sft$powerEstimate, #软阈值
maxBlockSize = 6000, #每个区块中的最大基因数
TOMType = "unsigned", #拓扑重叠矩阵的类型
minModuleSize = 30, #最小模块大小
reassignThreshold = 0, #将基因重新分配给模块的阈值
mergeCutHeight = 0.25, #指定合并模块的高度
numericLabels = TRUE, #模块标签是否应为数字
pamRespectsDendro = FALSE, #PAM 聚类是否应尊重树枝图
saveTOMs = F, #是否保存拓扑重叠矩阵
verbose = 3
)## Calculating module eigengenes block-wise from all genes
## Flagging genes and samples with too many missing values...
## ..step 1
## ..Working on block 1 .
## TOM calculation: adjacency..
## ..will not use multithreading.
## Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## ....clustering..
## ....detecting modules..
## ....calculating module eigengenes..
## ....checking kME in modules..
## ..removing 97 genes from module 1 because their KME is too low.
## ..removing 6 genes from module 2 because their KME is too low.
## ..removing 1 genes from module 4 because their KME is too low.
## ..removing 2 genes from module 6 because their KME is too low.
## ..removing 7 genes from module 8 because their KME is too low.
## ..removing 1 genes from module 9 because their KME is too low.
## ..removing 2 genes from module 10 because their KME is too low.
## ..removing 1 genes from module 12 because their KME is too low.
## ..removing 3 genes from module 13 because their KME is too low.
## ..removing 2 genes from module 15 because their KME is too low.
## ..removing 7 genes from module 17 because their KME is too low.
## ..removing 1 genes from module 21 because their KME is too low.
## ..removing 1 genes from module 22 because their KME is too low.
## ..removing 1 genes from module 25 because their KME is too low.
## ..removing 2 genes from module 27 because their KME is too low.
## ..removing 3 genes from module 29 because their KME is too low.
## ..removing 4 genes from module 30 because their KME is too low.
## ..removing 2 genes from module 32 because their KME is too low.
## ..merging modules that are too close..
## mergeCloseModules: Merging modules whose distance is less than 0.25
## Calculating new MEs...table(net$colors) ##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 246 1671 355 305 279 270 241 175 168 124 108 101 88 87 84 84
## 16 17 18 19 20 21 22 23 24 25 26 27
## 77 67 62 61 60 47 44 41 40 40 38 37所有的核心就在这一步，把输入的表达矩阵的几千个基因组归类成了几十个模块。大体思路：计算基因间的邻接性，根据邻接性计算基因间的相似性，然后推出基因间的相异性系数，并据此得到基因间的系统聚类树。然后按照混合动态剪切树的标准，设置每个基因模块最少的基因数目为30。
根据动态剪切法确定基因模块后，再次分析，依次计算每个模块的特征向量值，然后对模块进行聚类分析，将距离较近的模块合并为新的模块。
然后是分布法完成网络构建，仅供有探索精神的同学挑战。
构建加权共表达网络分为两步：
1. 计算邻近值，也是就是两个基因在不样品中表达量的表达相关系数(pearson correlation rho)，
2. 计算topology overlap similarity (TOM)。 WGCNA认为，只通过计算两个基因的表达相关系数构建共表达网络是不足够的。
于是他们用TOM表示两个基因在网络结构上的相似性，即两个基因如果具有相似的邻近基因，这两个基因更倾向于有相互作用。
参考 2.b.3 in https://labs.genetics.ucla.edu/horvath/htdocs/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/FemaleLiver-02-networkConstr-man.pdf
###################################################不运行#####################################################################################
#(1)网络构建 Co-expression similarity and adjacency
# 构建邻接矩阵用于表示基因之间的相似性
adjacency = adjacency(datExpr, power = sft$powerEstimate)
#(2) 邻近矩阵到拓扑矩阵的转换，Turn adjacency into topological overlap
# 邻接矩阵被转换成了拓扑矩阵（TOM），通常用于表示基因之间的拓扑相似性。
TOM = TOMsimilarity(adjacency);## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.# TOMsimilarity函数用于计算TOM。然后，通过计算 1 - TOM得到了一个关于基因之间差异性的矩阵 dissTOM
dissTOM = 1-TOM
# (3) 聚类拓扑矩阵 Call the hierarchical clustering function
# 使用层次聚类（hierarchical clustering）来对基因进行聚类
# as.dist(dissTOM)将上一步中得到的差异性矩阵转换成了一个距离矩阵，然后 hclust函数用于进行层次聚类
# method = "average"表示使用平均链接法进行聚类。聚类的结果将形成一个树状图（dendrogram）。
geneTree = hclust(as.dist(dissTOM), method = "average");
# Plot the resulting clustering tree (dendrogram)
sizeGrWindow(12,9)
## 这个时候的geneTree与一步法的 net$dendrograms[[1]] 性质类似，但是还需要进行进一步处理
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04);
#(4) 聚类分支的修整 dynamicTreeCut
# 对先前得到的聚类结果进行修整。
# 首先，设置了一个最小模块大小（ minModuleSize = 30），用于指定最小模块中基因的数量。
# 然后，使用 cutreeDynamic函数对基因树进行动态切割，得到了不同的模块，这些模块包含具有相似表达模式的基因。
# 最后，通过表格展示了每个模块中包含的基因数量。
#
# We like large modules, so we set the minimum module size relatively high:
minModuleSize = 30;
# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
deepSplit = 2, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize);## ..cutHeight not given, setting it to 0.998 ===> 99% of the (truncated) height range in dendro.
## ..done.table(dynamicMods)## dynamicMods
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 8 1274 214 193 176 163 160 159 146 146 138 135 124 120 116 115
## 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
## 107 101 97 93 93 78 76 75 67 66 66 65 65 63 63 54
## 32 33 34 35 36 37 38 39 40
## 52 50 45 44 43 41 40 36 33#4. 绘画结果展示
# Convert numeric lables into colors
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)## dynamicColors
## black blue brown cyan darkgreen
## 159 214 193 116 76
## darkgrey darkmagenta darkolivegreen darkorange darkred
## 67 45 50 66 78
## darkturquoise green greenyellow grey grey60
## 75 163 135 8 101
## lightcyan lightgreen lightyellow magenta mediumpurple3
## 107 97 93 146 33
## midnightblue orange orangered4 paleturquoise pink
## 115 66 36 54 146
## plum1 purple red royalblue saddlebrown
## 40 138 160 93 63
## salmon sienna3 skyblue skyblue3 steelblue
## 120 44 65 41 63
## tan turquoise violet white yellow
## 124 1274 52 65 176
## yellowgreen
## 43# Plot the dendrogram and colors underneath
#sizeGrWindow(8,6)
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors")
#5. 聚类结果相似模块的融合，Merging of modules whose expression profiles are very similar
#在聚类树中每一leaf是一个短线，代表一个基因，
#不同分之间靠的越近表示有高的共表达基因，将共表达极其相似的modules进行融合
# 计算了模块的特征基因（eigengenes），这些基因代表了模块的表达模式。
# 然后，计算了模块特征基因之间的差异性，将其用于层次聚类，得到了一个新的树状图。接着，通过设置一个相关性阈值（ MEDissThres = 0.25），将具有高相关性的模块合并在一起。
# 最后，通过表格展示了合并后的模块以及它们的特征基因。
# Calculate eigengenes
MEList = moduleEigengenes(datExpr, colors = dynamicColors)
MEs = MEList$eigengenes
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs);
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average");
# Plot the result
#sizeGrWindow(7, 6)
plot(METree, main = "Clustering of module eigengenes",
xlab = "", sub = "")
#选择有75%相关性的进行融合
MEDissThres = 0.25
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")
# Call an automatic merging function
merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)## mergeCloseModules: Merging modules whose distance is less than 0.25
## multiSetMEs: Calculating module MEs.
## Working on set 1 ...
## moduleEigengenes: Calculating 41 module eigengenes in given set.
## multiSetMEs: Calculating module MEs.
## Working on set 1 ...
## moduleEigengenes: Calculating 34 module eigengenes in given set.
## multiSetMEs: Calculating module MEs.
## Working on set 1 ...
## moduleEigengenes: Calculating 33 module eigengenes in given set.
## Calculating new MEs...
## multiSetMEs: Calculating module MEs.
## Working on set 1 ...
## moduleEigengenes: Calculating 33 module eigengenes in given set.# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new merged modules:
mergedMEs = merge$newMEsstep4: 模块可视化
这里用不同的颜色来代表那些所有的模块，其中灰色默认是无法归类于任何模块的那些基因，如果灰色模块里面的基因太多，那么前期对表达矩阵挑选基因的步骤可能就不太合适。
# Convert labels to colors for plotting
mergedColors = labels2colors(net$colors)
table(mergedColors)## mergedColors
## black blue brown cyan darkgreen
## 175 355 305 84 44
## darkgrey darkorange darkred darkturquoise green
## 40 38 47 41 270
## greenyellow grey grey60 lightcyan lightgreen
## 101 246 67 77 62
## lightyellow magenta midnightblue orange pink
## 61 124 84 40 168
## purple red royalblue salmon tan
## 108 241 60 87 88
## turquoise white yellow
## 1671 37 279# Plot the dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
## assign all of the gene to their corresponding module
## hclust for the genes.
基因的模块可视化
这里的重点就是plotDendroAndColors函数，它接受一个聚类的对象，以及该对象里面包含的所有个体所对应的颜色。比如对表达矩阵进行hclust之后，加上表达矩阵里面所有样本的分组信息对应的颜色，也是可以用plotDendroAndColors函数可视化的，比如下面的样品图：
# ###################################################不运行#####################################################################################
# #明确样本数和基因数
# nGenes = ncol(datExpr)
# nSamples = nrow(datExpr)
# #首先针对样本做个系统聚类树
# datExpr_tree<-hclust(dist(datExpr), method = "average")
# par(mar = c(0,5,2,0))
# plot(datExpr_tree, main = "Sample clustering", sub="", xlab="", cex.lab = 2,
# cex.axis = 1, cex.main = 1,cex.lab=1)
# ## 如果这个时候样本是有性状，或者临床表型的，可以加进去看看是否聚类合理
# #针对前面构造的样品矩阵添加对应颜色
# sample_colors <- numbers2colors(as.numeric(factor(datTraits$Tumor.Type)),
# colors = c("white","blue","red","green"),signed = FALSE)
# ## 这个给样品添加对应颜色的代码需要自行修改以适应自己的数据分析项目。
# # sample_colors <- numbers2colors( datTraits ,signed = FALSE)
# ## 如果样品有多种分类情况，而且 datTraits 里面都是分类信息，那么可以直接用上面代码，当然，这样给的颜色不明显，意义不大。
# #构造10个样品的系统聚类树及性状热图
# par(mar = c(1,4,3,1),cex=0.8)
# plotDendroAndColors(datExpr_tree, sample_colors,
# groupLabels = colnames(sample),
# cex.dendroLabels = 0.8,
# marAll = c(1, 4, 3, 1),
# cex.rowText = 0.01,
# main = "Sample dendrogram and trait heatmap")上面给样本进行聚类的代码可以不运行，其实跟WGCNA本身关系不大。
样本的聚类可视化
可以看到这些乳腺癌的细胞系的表达谱聚类情况并不是完全与其分类匹配，所以仅仅是根据样本的分组信息做差异分析并不完全准确。
step5:模块和性状的关系
## step 5 (最重要的) 模块和性状的关系
## 这一步主要是针对于连续变量，如果是分类变量，需要转换成连续变量方可使用
table(datTraits$subtype)##
## Basal Claudin-low Luminal Non-malignant unknown
## 14 6 27 5 4if(T){
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
design=model.matrix(~0+ datTraits$subtype)
colnames(design)=levels(datTraits$subtype)
moduleColors <- labels2colors(net$colors)
# Recalculate MEs with color labels 计算表达矩阵(datExpr)中的基因模块的特征向量（MEs）。
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0); ##不同颜色的模块的ME值矩 (样本vs模块)
# 计算MEs与样本亚型(datTraits$subtype)之间的相关性。
moduleTraitCor = cor(MEs, design , use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
# 使用热图和条形图将这些关系可视化。
sizeGrWindow(10,6)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
png("step5-Module-trait-relationships.png",width = 800,height = 1200,res = 120)
par(mar = c(6, 8.5, 3, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = colnames(design),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = greenWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
zlim = c(-1,1),
main = paste("Module-trait relationships"))
dev.off()
# 除了上面的热图展现形状与基因模块的相关性外
# 还可以是条形图,但是只能是指定某个形状
# 或者自己循环一下批量出图。
#从样本亚型(datTraits$subtype)中提取Luminal亚型。
Luminal = as.data.frame(design[,3]);
names(Luminal) = "Luminal"
y=Luminal
#计算Luminal亚型与表达矩阵(datExpr)之间的相关性。
GS1=as.numeric(cor(y,datExpr, use="p"))
GeneSignificance=abs(GS1)
# Next module significance is defined as average gene significance.
#计算基因模块的重要性。
ModuleSignificance=tapply(GeneSignificance,
moduleColors, mean, na.rm=T)
sizeGrWindow(8,7)
par(mfrow = c(1,1))
# 如果模块太多，下面的展示就不友好
# 不过，我们可以自定义出图。
plotModuleSignificance(GeneSignificance,moduleColors)
}## Warning in greenWhiteRed(50): WGCNA::greenWhiteRed: this palette is not suitable for people
## with green-red color blindness (the most common kind of color blindness).
## Consider using the function blueWhiteRed instead.通过模块与各种表型的相关系数，可以很清楚的挑选自己感兴趣的模块进行下游分析了。这个图就是把moduleTraitCor这个矩阵给用热图可视化一下。
模块和性状的关系
因为一些历史遗留问题，这个图片缺乏X轴的标记。
从上图已经可以看到跟乳腺癌分类相关的基因模块了，包括"Basal" "Claudin-low" "Luminal" "Non-malignant" "unknown" 这5类所对应的不同模块的基因列表。可以看到每一种乳腺癌都有跟它强烈相关的模块，可以作为它的表达signature，模块里面的基因可以拿去做下游分析。我们看到Luminal表型跟棕色的模块相关性高达0.86，而且极其显著的相关，所以值得我们挖掘，这个模块里面的基因是什么，为什么如此的相关呢？
step6:感兴趣性状的模块的具体基因分析
性状跟模块虽然求出了相关性，可以挑选最相关的那些模块来分析，但是模块本身仍然包含非常多的基因，还需进一步的寻找最重要的基因。所有的模块都可以跟基因算出相关系数，所有的连续型性状也可以跟基因的表达值算出相关系数。主要参考资料：PDF document, R script 如果跟性状显著相关基因也跟某个模块显著相关，那么这些基因可能就非常重要。
首先计算模块与基因的相关性矩阵
# names (colors) of the modules
modNames = substring(names(MEs), 3) #提取从第三个字符开始的子字符串
## 算出每个模块跟基因的皮尔森相关系数矩阵
## MEs是每个模块在每个样本里面的值
## datExpr是每个基因在每个样本的表达量
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p")); #使用Pearson相关系数计算datExpr和MEs之间的相关性
#使用学生t分布计算geneModuleMembership中每个相关系数的p值
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
names(geneModuleMembership) = paste("MM", modNames, sep=""); #将"MM"与每个元素的名称进行连接来重命名geneModuleMembership的列。
names(MMPvalue) = paste("p.MM", modNames, sep=""); #将"p.MM"与每个元素的名称进行连接来重命名MMPvalue的列再计算性状与基因的相关性矩阵
## 只有连续型性状才能只有计算
## 这里把是否属于 Luminal 表型这个变量用0,1进行数值化。
Luminal = as.data.frame(design[,3]); #提取了设计对象的第三列，并将其转换为数据框，然后将其赋值给Luminal。
names(Luminal) = "Luminal"
geneTraitSignificance = as.data.frame(cor(datExpr, Luminal, use = "p")); #使用Pearson相关系数计算datExpr和Luminal之间的相关性
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples)); #使用学生t分布计算geneTraitSignificance中每个相关系数的p值
names(geneTraitSignificance) = paste("GS.", names(Luminal), sep=""); #将"GS."与Luminal中每个列名称进行连接来重命名geneTraitSignificance的列
names(GSPvalue) = paste("p.GS.", names(Luminal), sep=""); #将"p.GS."与Luminal中每个列名称进行连接来重命名GSPvalue的列最后把两个相关性矩阵联合起来,指定感兴趣模块进行分析
module = "brown"
column = match(module, modNames);
moduleGenes = moduleColors==module;
sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance for Luminal",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
模块和性状里面的指定基因的相关性比较
可以看到这些基因不仅仅是跟其对应的模块高度相关，而且是跟其对应的性状高度相关，进一步说明了基因值得深度探究。
step7:网络的可视化
主要参考资料：PDF document, R script
首先针对所有基因画热图
# 主要是可视化 TOM矩阵，WGCNA的标准配图
# 然后可视化不同 模块 的相关性 热图
# 不同模块的层次聚类图
# 还有模块诊断，主要是 intramodular connectivity
if(T){
#获取数据集中基因和样本的数量
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
#获取基因的树状图
geneTree = net$dendrograms[[1]];
#根据基因的表达水平计算基因之间的相似性。power 参数用于控制基因之间的连接程度。
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);
plotTOM = dissTOM^7; #将相似性矩阵提高到7次方。
diag(plotTOM) = NA; #将矩阵的对角线元素设置为NA。
#TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
nSelect = 400 #随机从数据集中选择400个基因。
# For reproducibility, we set the random seed设置随机种子以保证可重复性
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select]; #从相似性矩阵中选择一部分基因。
# There’s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
#使用层次聚类重新对选择的基因进行聚类。
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select]; #根据模块成员关系为每个基因分配颜色。
# Open a graphical window
sizeGrWindow(9,9)
# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# the color palette; setting the diagonal to NA also improves the clarity of the plot
#设置新的相似性矩阵以进行绘图。对角线元素设置为NA以提高可读性。
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
png("step7-Network-heatmap.png",width = 800,height = 600)
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")
dev.off()
# Recalculate module eigengenes
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
## 只有连续型性状才能只有计算
## 这里把是否属 Luminal 表型这个变量0,1进行数值化
Luminal = as.data.frame(design[,3]);
names(Luminal) = "Luminal"
# Add the weight to existing module eigengenes
MET = orderMEs(cbind(MEs, Luminal))
# Plot the relationships among the eigengenes and the trait
sizeGrWindow(5,7.5);
par(cex = 0.9)
png("step7-Eigengene-dendrogram.png",width = 800,height = 600)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
= 90)
dev.off()
# Plot the dendrogram
sizeGrWindow(6,6);
par(cex = 1.0)
## 模块的进化树
png("step7-Eigengene-dendrogram-hclust.png",width = 800,height = 600)
plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
plotHeatmaps = FALSE)
dev.off()
# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
par(cex = 1.0)
## 性状与模块热
png("step7-Eigengene-adjacency-heatmap.png",width = 800,height = 600)
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),
plotDendrograms = FALSE, xLabelsAngle = 90)
dev.off()
}## TOM calculation: adjacency..
## ..will not use multithreading.
## Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.## pdf
## 2这个非常消耗计算资源和时间，所以建议选取其中部分基因作图即可，我就没有画，而且根据下面的代码选取部分基因来作图！
然后随机选取部分基因作图
nSelect = 400
# For reproducibility, we set the random seed
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
# There’s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
# Open a graphical window
sizeGrWindow(9,9)
# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# the color palette; setting the diagonal to NA also improves the clarity of the plot
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")
模块热图
这个图凑数的意义居多，如果是把全部基因画上去，可以很清楚的看到各个区块颜色差异。
最后画模块和性状的关系
moduleEigengenes 是 R 语言中的一个函数，用于计算模块特征基因。它的输入是一个表达数据矩阵和一个模块颜色向量，返回一个列表，其中包含以下组件：
eigengenes：模块特征基因的数据框，每一列对应一个特征基因。列名以相应的颜色为前缀，例如 MEturquoise 等。如果 returnValidOnly==FALSE，则计算失败的模块特征基因的所有组成部分都设置为 NA。 averageExpr：如果 align == “along average”，则包含每个模块中的平均归一化表达式的数据框。列名以相应的颜色为前缀，例如 AEturquoise 等。 varExplained：数据框，其中每一列对应一个模块，其中组件 varExplained [PC, module] 给出了主成分 PC 解释的模块 module 的方差。计算是精确的，不受计算的主成分数量影响。在此数据框中记录最多 10 个方差解释值。 nPC：输入 nPC 的副本。 validMEs：布尔向量。每个组件（对应于数据中的列）如果相应的特征基因有效，则为 TRUE；否则为 FALSE。有效特征基因包括主成分和它们的 hubgene 近似值。当 returnValidOnly==FALSE 时，根据定义，所有返回的特征基因都是有效的，并且 validMEs 的条目都为 TRUE。 validColors：输入颜色的副本，其中与无效模块对应的条目设置为灰色（如果给定），否则为 0（如果颜色是数字）；否则为“grey”。 allOK：布尔标志，表示是否已正确计算所有特征基因，无论是作为主成分还是作为 hubgene 平均近似值。 allPC：布尔标志，表示所有返回的特征基因是否都是主成分。 isPC：布尔向量。每个组件（对应于 eigengenes 中的列）如果相应的特征基因是第一个主成分，则为 TRUE；否则为 FALSE 或无效。 isHub：布尔向量。每个组件（对应于 eigengenes 中的列）如果相应的特征基因是 hubgene 近似值，则为 TRUE；否则为 FALSE 或无效。 validAEs：布尔向量。每个组件（对应于 eigengenes 中的列）如果相应模块平均表达式有效，则为 TRUE。 allAEOK：布尔标志，表示所有返回的模块平均表达式是否包含有效数据。请注意，returnValidOnly==TRUE 不意味着 allAEOK==TRUE：如果其相应特征基因已正确计算，则可能返回一些无效平均表达式。
# Recalculate module eigengenes
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
## 只有连续型性状才能只有计算
## 这里把是否属于 Luminal 表型这个变量用0,1进行数值化。
Luminal = as.data.frame(design[,3]);
names(Luminal) = "Luminal"
# Add the weight to existing module eigengenes
MET = orderMEs(cbind(MEs, Luminal))
# Plot the relationships among the eigengenes and the trait
sizeGrWindow(5,7.5);
par(cex = 0.9)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
= 90)
# Plot the dendrogram
sizeGrWindow(6,6);
par(cex = 1.0)
## 模块的聚类图
plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
plotHeatmaps = FALSE)
# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
par(cex = 1.0)
## 性状与模块热图
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),
plotDendrograms = FALSE, xLabelsAngle = 90)
性状与模块热图
step8:提取指定模块的基因名
## step 8
# 主要是关心具体某个模块内部的基因
if(T){
# Select module
module = "brown";
# Select module probes
probes = colnames(datExpr) ## 我们例子里面的probe就是基因，获取数据集中所有基因的名称
inModule = (moduleColors==module); #确定每个基因是否属于所选模块。
modProbes = probes[inModule]; #获取所选模块中的所有基因。
head(modProbes)
# 如果使用WGCNA包自带的热图就很丑。
which.module="brown";
dat=datExpr[,moduleColors==which.module ] #获取所选模块中的所有基因表达数据。
#绘制热图。其中，scale 函数用于对数据进行归一化，nrgcols 参数用于设置颜色数量，rlabels 和 clabels 参数用于显示行和列标签，rcols 参数用于设置颜色，title 参数用于设置标题。
plotMat(t(scale(dat)),nrgcols=30,rlabels=T,
clabels=T,rcols=which.module,
title=which.module )
datExpr[1:4,1:4]
#获取所选模块中的所有基因表达数据，并将其转置。
dat=t(datExpr[,moduleColors==which.module ] )
library(pheatmap)
pheatmap(dat ,show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
n=t(scale(t(log(dat+1)))) # 'scale'可以对log-ratio数值进行归一化。对数据进行归一化和对数变换，并将其转置。
#对数据进行截断，以便更好地显示热图。
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
#绘制带有注释的热图。其中，datTraits$subtype 是样本类型信息，用于为每个样本添加注释。
pheatmap(n,show_colnames =F,show_rownames = F)
group_list=datTraits$subtype
ac=data.frame(g=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac )
# 可以很清晰的看到，所有的形状相关的模块基因
# 其实未必就不是差异表达基因。
}
有了基因信息，下游分析就很简单了。包括GO/KEGG等功能数据库的注释
image-20190928155929933
Step9: 模块的导出
主要模块里面的基因直接的相互作用关系信息可以导出到cytoscape,VisANT等网络可视化软件。
# Recalculate topological overlap
TOM = TOMsimilarityFromExpr(datExpr, power = 6); ## TOM calculation: adjacency..
## ..will not use multithreading.
## Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.# Select module
module = "brown";
# Select module probes
probes = colnames(datExpr) ## 我们例子里面的probe就是基因名
inModule = (moduleColors==module);
modProbes = probes[inModule];
## 也是提取指定模块的基因名
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)
## 模块对应的基因关系矩阵 首先是导出到VisANT
vis = exportNetworkToVisANT(modTOM,
file = paste("VisANTInput-", module, ".txt", sep=""),
weighted = TRUE,
threshold = 0)然后是导出到cytoscape
cyt = exportNetworkToCytoscape(
modTOM,
edgeFile = paste("CytoscapeInput-edges-", paste(module, collapse="-"), ".txt", sep=""),
nodeFile = paste("CytoscapeInput-nodes-", paste(module, collapse="-"), ".txt", sep=""),
weighted = TRUE,
threshold = 0.02,
nodeNames = modProbes,
nodeAttr = moduleColors[inModule]
);如果模块包含的基因太多，网络太复杂，还可以进行筛选，比如：
nTop = 30;
IMConn = softConnectivity(datExpr[, modProbes]);
top = (rank(-IMConn) <= nTop)
filter <- modTOM[top, top]后面就是cytoscape自身的教程了，这里不再赘述，我博客有比较详尽的介绍。