WGCNA应用实战二
大家好，我是风风。今天我们要继续我们的WGCNA系列。在之前的推文中，我们已经一起做了对单个数据集进行了WGCNA分析并完整地走完了单个WGCNA的流程和针对多数据集联合分析的共识性分析，同时，在构建网络的过程中，我们也使用了“自动构建拓扑异构网络”法和“逐步法构建拓扑异构网络”两种方法，那么我们WGCNA的常用知识就算基本完结了，接下来无非就是利用WGCNA进行实战的内容。我们将通过三次不同角度的实战，利用我们自己的数据进行分析，看看如何将WGCNA应用到我们的文章之中。
在上一次的推文中，我们尝试了将WGCNA分析和免疫浸润分析的结果联合起来，寻找表达数据和免疫特征之间的关系，从而去进行相应的文章分析，这是经典的WGCNA联系表型的方法，把免疫浸润分析的结果替换成其他的表型也一样可以分析。

在我们前面的基础内容中，曾经使用共识性网络，来进行WGCNA联合分析，今天我们一样来看一下多数据集的联合分析——由Jeremy Miller提出的对来自多个芯片数据集的数据进行荟萃分析，Jeremy Miller提出的这种方法也被WGCNA的作者将链接放在了WGCNA的官方文档，因此权威性和准确性都很强。作者认为，对于芯片数据而言，分析的内容主要分为两种：其一是差异表达分析，也就是我们“挑圈联靠”中的“挑”，找到差异表达的分子，A基因在实验组中高表达，B基因实验组中低表达这一类内容；其二便是共表达分析，注意啊，不是常规的相关性分析，而是共表达分析（熟熟悉咱们“生信全书”的学员应该能够回想起来，WGCNA的课程内容是放在了“挑圈联靠”中“圈”部分的内容中，而相关性分析又是在“联”部分的内容中）。共表达的过程并不是仅仅进行相关性分析，在前面WGCNA系列的第一个推文中，我们就给大家展示过了作者关于WGCNA的步骤流程图了。

那么，当同样的内容有多个芯片数据，或者说相似的研究中的多个芯片数据，可能会出现的一个问题，即“如果A组和B组都进行了芯片数据的分析研究并报告了一些结果，这些结果之间的兼容性怎么样？”对于差异表达分析的结果，可以使用RRA分析，而对于WGCNA的分析的结果，我们一样也可以达到类似的目的，我们来看看具体的操作步骤。
WGCNA系列传送门
仅凭一个分析就能发3.5+SCI！万字长文教你彻底学会WGCNA！不分享说不过去吧？
比另一半更爱你！被你们问的最多的生信分析！一个分析就是一篇3+SCI！高级！给力！
被问爆的WGCNA教程，高分秘籍，有手就会！（附代码）
决定了！把免疫做到新高度就用这个分析！|（附代码 )
数据处理
在进行相应的分析之前，我们同样先对数据进行预处理，使用的数据和相关代码放在公众号后台，回复关键词即可获取：
rm(list = ls())# BiocManager::install("WGCNA") 安装R包# BiocManager::install("Hmisc") # BiocManager::install("flashClust") # BiocManager::install("qvalue") # BiocManager::install("dynamicTreeCut") # BiocManager::install("impute") options(stringsAsFactors = FALSE)library(WGCNA)library(Hmisc)library(flashClust)library(dynamicTreeCut)library(impute)library(ggplot2)# 加载数据load("metaAnalysisData.RData")# 查看数据# datExprA1和datExprA22是来自Illumina Human ref-12平台的两个数据集datExprA1[1:3, 1:3]## C.1385.CA3.M.81.4.O C.3795.CA1.F.72.4.V C.1761.CA1b.M.86.4.O## ILMN_2055271 162.8675 216.7612 185.6294## ILMN_1653355 180.3644 249.8743 148.2917## ILMN_2359168 1246.7213 443.4179 1787.8173datExprA2[1:3, 1:3]## A.4147.CA3.M.66.4.V A.704.CA1.M.81.4.O A.4125.CA1.F.61.4.V## ILMN_2055271 148.5686 183.1183 213.0718## ILMN_1653355 144.7133 197.0489 183.0671## ILMN_2359168 803.9261 500.7845 591.5554# datExprB1和datExprB2是来自Illumina Human ref-12和Affymetrix HG-U133A的数据集datExprB1[1:3, 1:3]## C.1385.CA3.M.81.4.O C.3795.CA1.F.72.4.V C.1761.CA1b.M.86.4.O## ILMN_2055271 162.8675 216.7612 185.6294## ILMN_1653355 180.3644 249.8743 148.2917## ILMN_2359168 1246.7213 443.4179 1787.8173datExprB2[1:3, 1:3]## GSM21203.cel GSM21204.cel GSM21208.cel## 1007_s_at 2590.52064 2071.477139 2110.995464## 1053_at 36.68135 6.354923 1.658427## 117_at 32.07228 47.248082 89.273107# probesI和probesA是用于Human ref-12和HG-U133A平台的探针length(probesI)## [1] 21696length(probesA)## [1] 20198# genesI和genesA是probesI和probesA对应的基因symbol名称length(genesI)## [1] 21696length(genesA)## [1] 20198# 接下来我们需要对数据集进行预处理，以确保数据之间具有可比性。这里我们需要注意有两种情况：其一是如果我们的两份甚至多份数据集都来自同一平台，那么我们就不再需要对数据进行额外的处理，直接进行注释后就可以进行分析了；其二是如果我们感兴趣的数据集来自不同的平台，比如其中一个数据集来自Illumina平台，另一个来自Affymetrix平台，如果是这种情况，则我们需要以某种方式匹配探针，我们可以根据基因symbol为每个数据集中的每个探针确定对应的symbol名称，然后使用基因symbol代替探针ID作为每个数据集中的标识符，然后再对两个数据集进行交集处理使两个数据集所包含的基因symbol一致，我们可以使用以下代码达到目的：source("collapseRows_NEW.R") datExprB1g <- (collapseRows(datExprB1, genesI, probesI))[[1]] datExprB2g <- (collapseRows(datExprB2, genesA, probesA))[[1]]## Warning: row names of input data and probes not identical...## ... Attempting to proceed anyway. Check results carefully.# collapseRows，collapse是折叠的意思，collapseRows就是对行进行折叠处理，简单地理解就是注释，对datExprB1和datExprB2的行名根据genesI/probesI和genesA/probesA进行匹配注释，从而将探针ID转换为基因symbolhead(datExprB1g)[1:6, 1:6]## C.1385.CA3.M.81.4.O C.3795.CA1.F.72.4.V C.1761.CA1b.M.86.4.O## 1-Mar 195.7141 190.9397 210.2140## 1-Sep 189.3278 153.0674 164.3864## 10-Sep 189.7120 151.4790 179.8286## 11-Sep 1226.4431 1064.2611 1177.9998## 12-Sep 297.6668 1180.4090 341.5574## 15-Sep 2907.7309 1887.1917 2986.0669## C.1013.CA3.M.90.4.O C.3740.CA1.F.88.4.V C.1716.CA3.M.90.4.O## 1-Mar 191.1562 230.7255 220.8315## 1-Sep 162.0791 174.7665 189.8200## 10-Sep 194.2755 206.6645 210.5338## 11-Sep 988.9611 1790.5023 1618.3355## 12-Sep 308.3526 352.7494 445.1497## 15-Sep 2217.0365 3300.1172 2918.7352head(datExprB2g)[1:6, 1:6]## GSM21203.cel GSM21204.cel GSM21208.cel GSM21209.cel GSM21210.cel## A1CF 195.79707 186.38870 239.23413 203.23650 121.78074## A2BP1 870.26767 934.05645 2406.66220 851.55378 1126.94997## A2M 913.00166 1361.38048 1263.90340 1151.67477 1823.09367## A4GALT 16.26832 66.94742 42.68259 47.57441 15.61746## A4GNT 128.87320 102.63451 98.61871 125.92344 76.22605## AACS 557.08256 496.86176 503.28366 461.01426 405.07454## GSM21211.cel## A1CF 118.99223## A2BP1 950.84910## A2M 1076.33419## A4GALT 37.58282## A4GNT 90.73498## AACS 421.07370# 接着我们处理datExprA1和datExprA2数据集，使用intersect()函数进行取交集处理commonProbesA <- intersect (rownames(datExprA1), rownames(datExprA2)) head(commonProbesA)## [1] "ILMN_2055271" "ILMN_1653355" "ILMN_2359168" "ILMN_1787689" "ILMN_1745607"## [6] "ILMN_2136495"# 同样提取交集行名对应的表达矩阵datExprA1p <- datExprA1[commonProbesA, ] dim(datExprA1p)## [1] 21696 32datExprA2p <- datExprA2[commonProbesA, ] dim(datExprA2p)## [1] 21696 31# 数据之间是否有可比性，其中最重要的一点是两份数据集中是否包含的元素相同，就像我们构建模型使用的训练集和验证集一样，如果构建模型的分子只在训练集中存在，而在验证集中不存在，那么也不存在验证的可能性，因此取交集是关键前提步骤# 接着我们同样处理datExprB1g和datExprB2g数据集，使用intersect()函数进行取交集处理commonGenesB <- intersect (rownames(datExprB1g), rownames(datExprB2g)) head(commonGenesB)## [1] "A2BP1" "A2M" "A4GALT" "A4GNT" "AACS" "AAK1"datExprB1g <- datExprB1g[commonGenesB, ] dim(datExprB1g)## [1] 8057 32datExprB2g <- datExprB2g[commonGenesB, ] dim(datExprB2g)## [1] 8057 26# 这样的话，用于比较的数据集中的每一行就都对应于相同的探针/基因了
数据集间的可比性分析
那究竟要如何快速评估两个数据集之间的可比性呢？利用WGCNA的方法就是将两个数据集之间的平均基因表达和总体连通性的度量关联起来。这些属性的相关性越高，我们就越有可能在分析的后续阶段找到两个数据集之间的相似性。具体操作如下：
# 联合免疫浸润和WGCNA分析softPower<- 10 # 设置软阈值rankExprA1<- rank(rowMeans(datExprA1p)) rankExprA2<- rank(rowMeans(datExprA2p)) random5000<- sample(commonProbesA,5000) rankConnA1<- rank(softConnectivity(t(datExprA1p[random5000,]),type="signed",power=softPower))## softConnectivity: FYI: connecitivty of genes with less than 11 valid samples will be returned as NA.## ..calculating connectivities..rankConnA2<- rank(softConnectivity(t(datExprA2p[random5000,]),type="signed",power=softPower))## softConnectivity: FYI: connecitivty of genes with less than 11 valid samples will be returned as NA.## ..calculating connectivities..rankExprB1<- rank(rowMeans(datExprB1g)) rankExprB2<- rank(rowMeans(datExprB2g)) random5000<- sample(commonGenesB,5000) # 构建连接度网络rankConnB1<- rank(softConnectivity(t(datExprB1g[random5000,]),type="signed",power=softPower))## softConnectivity: FYI: connecitivty of genes with less than 11 valid samples will be returned as NA.## ..calculating connectivities..rankConnB2<- rank(softConnectivity(t(datExprB2g[random5000,]),type="signed",power=softPower))## softConnectivity: FYI: connecitivty of genes with less than 9 valid samples will be returned as NA.## ..calculating connectivities..# pdf("generalNetworkProperties.pdf", height=10, width=9) par(mfrow=c(2,2)) verboseScatterplot(rankExprA1,rankExprA2,col = "steelblue",xlab="Ranked Expression (A1)", ylab="Ranked Expression (A2)") verboseScatterplot(rankConnA1,rankConnA2,col = "steelblue",xlab="Ranked Connectivity (A1)", ylab="Ranked Connectivity (A2)") verboseScatterplot(rankExprB1,rankExprB2,col = "steelblue",xlab="Ranked Expression (B1)", ylab="Ranked Expression (B2)") verboseScatterplot(rankConnB1,rankConnB2,col = "steelblue",xlab="Ranked Connectivity (B1)", ylab="Ranked Connectivity (B2)")
# dev.off()# 针对这个结果图片，可以给出如下解读：1）如果相关性为正且p值具有显著统计学意义，说明数据集之间的结果可以比较；2）A的相关性和p-values比b的相关性和p-values要好，这是因为两个A数据集是来源于同一个平台运行的。因此我们可以得出结论：来自不同平台的数据集不如来自同一平台的数据集具有可比性，但它们之间仍然具有可比性。
进行WGCNA分析
接下来我们使用来自同一平台的A1和A2进行后续的分析。为了计算方便，我们首先选择数据集中表达最多的5000个探针，然后对相同基因名的只保留其中一个探针：
# 使用rank函数排序后提取表达量最大的前5000个探针keepGenesExpr <- rank(-rowMeans(datExprA1p)) <= 5000datExprA1g <- datExprA1p[keepGenesExpr, ] datExprA2g <- datExprA2p[keepGenesExpr, ] # 进行匹配注释keepGenesDups <- (collapseRows(datExprA1g, genesI, probesI))[[2]]## Warning: row names of input data and probes not identical...## ... Attempting to proceed anyway. Check results carefully.head(keepGenesDups)## group selectedRowID ## 11-Sep "11-Sep""ILMN_1788778"## 15-Sep "15-Sep""ILMN_2407082"## 2-Mar "2-Mar""ILMN_1703142"## 2-Sep "2-Sep""ILMN_2365711"## 3-Sep "3-Sep""ILMN_1746673"## 4-Mar "4-Mar""ILMN_1789991"datExprA1g <- datExprA1g[keepGenesDups[,2], ] datExprA2g <- datExprA2g[keepGenesDups[,2], ] rownames(datExprA1g) <- rownames(datExprA2g) <- keepGenesDups[,1] # 接下来，我们将计算运行WGCNA所需的所有值，这些都是我们前面WGCNA分析使用过的步骤了# 对数据集A1进行WGCNA分析adjacencyA1 <- adjacency(t(datExprA1g), power = softPower, type = "signed")diag(adjacencyA1) = 0dissTOMA1 <- 1-TOMsimilarity(adjacencyA1, TOMType = "signed")## ..connectivity..## ..matrix multiplication (system BLAS)..## ..normalization..## ..done.geneTreeA1 <- flashClust(as.dist(dissTOMA1),                        method = "average") # 使用同样的方法对数据集A2进行分析adjacencyA2 <- adjacency(t(datExprA2g), power = softPower,                        type = "signed")diag(adjacencyA2) <- 0dissTOMA2 <- 1-TOMsimilarity(adjacencyA2, TOMType = "signed")## ..connectivity..## ..matrix multiplication (system BLAS)..## ..normalization..## ..done.geneTreeA2 <- flashClust(as.dist(dissTOMA2), method = "average") save.image("tutorial.RData")  # 运行时间较长，保留分析结果方便后续分析# 进行可视化# pdf("dendrogram.pdf", height = 6, width = 16) par(mfrow=c(1,2)) plot(geneTreeA1, xlab="", sub="", main="Gene clustering on TOM-based dissimilarity (A1)", labels=FALSE, hang=0.04); plot(geneTreeA2, xlab="", sub="", main="Gene clustering on TOM-based dissimilarity (A2)", labels=FALSE,hang=0.04)
# dev.off() # 从聚类图的结果来看，这两份数据的聚类结果都很好，合理推测后续进行WGCNA的分析也会很好。# 接下来我们将根据数据集A1来确定模块mColorh=NULL for(ds in 0:3){ tree = cutreeHybrid(dendro = geneTreeA1, pamStage=FALSE, minClusterSize = (30-3*ds), cutHeight = 0.99, deepSplit = ds, distM = dissTOMA1) mColorh = cbind(mColorh,labels2colors(tree$labels))}## ..done.## ..done.## ..done.## ..done.# pdf("Module_choices.pdf", height=10, width=25) plotDendroAndColors(geneTreeA1,mColorh,paste("dpSplt =", 0:3), main = "",dendroLabels=FALSE);
# dev.off() # 根据结果，这里我们选择dpSplt = 0的结果（此处主要根据自己研究目的进行选择，比较自由，没有统一标准）modulesA1 =  mColorh[, 1] # 接下来，我们计算可视化的主要数据。第一个数据是模块特征基因(ME)，代表模块中所有基因最大方差百分比的单一值。换句话说，如果我们显示X模块的ME发挥了某一类功能或者表型，那么X模块中的大多数基因也很有可能发挥了某一类同样的功能或者表型。PCs1A<- moduleEigengenes(t(datExprA1g),colors=modulesA1) ME_1A<- PCs1A$eigengenes distPC1A<- 1-abs(cor(ME_1A,use = "p")) distPC1A<- ifelse(is.na(distPC1A), 0,distPC1A)pcTree1A<- hclust(as.dist(distPC1A),method = "a") MDS_1A<- cmdscale(as.dist(distPC1A), 2) colorsA1 <- names(table(modulesA1)) save.image("tutorial.RData") # pdf("ModuleEigengeneVisualizations.pdf", height=6, width=6) par(mfrow=c(1,1), mar=c(0, 3, 1, 1) + 0.1, cex=1) plot(pcTree1A,xlab="",ylab="",main="",sub="")
plot(MDS_1A,col= colorsA1, main="MDS plot",cex=2, pch=19)
ordergenes = geneTreeA1$order plotMat(scale(log(datExprA1g[ordergenes,])),rlabels= modulesA1[ordergenes], clabels= colnames(datExprA1g),rcols = modulesA1[ordergenes])
for (which.module in names(table(modulesA1))){ ME = ME_1A[, paste("ME",which.module, sep="")] barplot(ME, col=which.module, main="", cex.main=2, ylab="eigengene expression", xlab="array sample") }
# dev.off()# 第一个聚类图是模块特征基因的树状图:同一分支中的模块具有相对相似的表达，第二个图是模块特征基因的多维缩放(MDS)图，组合在一起的模块具有相对相似的表达式。第三个图是所有基因的热图，每一列代表一个样本(x轴)，可以根据表型信息或聚类对其进行自定义排序。剩下的图是每一个模块的ME相对表达图。# 现在我们有了所有的WGCNA变量和模块定义，我们可以开始评估网络A1中的模块在网络A2中的保守性。作为定性评估，我们将A1中的模块施加到数据集A2的网络上，然后绘制得到的网络，也就是用A2的数据来验证A1的网络：# pdf("Final_modules.pdf", height=8, width=12) plotDendroAndColors(geneTreeA1, modulesA1,"Modules", dendroLabels=FALSE, hang=0.03, addGuide=TRUE, guideHang=0.05,main="Gene dendrogram and module colors (A1)")
plotDendroAndColors(geneTreeA2, modulesA1, "Modules", dendroLabels=FALSE, hang=0.03, addGuide=TRUE, guideHang=0.05,main="Gene dendrogram and module colors (A2)")
# dev.off()# 为了量化这个保守性分析的结果，我们利用了WGCNA包中内置的另一个函数:modulePreservation，这俄格函数使用多种策略评估一项研究中的模块在另一项研究中的保守性情况，并输出单个z分数：multiExpr = list(A1=list(data=t(datExprA1g)), A2=list(data=t(datExprA2g))) multiColor = list(A1 = modulesA1) mp=modulePreservation(multiExpr, multiColor, referenceNetworks = 1, verbose = 3, networkType = "signed", nPermutations = 10, maxGoldModuleSize = 100, maxModuleSize = 400)## ..checking data for excessive amounts of missing data..## Flagging genes and samples with too many missing values...## ..step 1## Flagging genes and samples with too many missing values...## ..step 1## ..unassigned 'module' name: grey ## ..all network sample 'module' name: gold## ..calculating observed preservation values## ..calculating permutation Z scores## ..Working with set 1 as reference set## ....working with set 2 as test set## ......working on permutation 1## ......working on permutation 2## ......working on permutation 3## ......working on permutation 4## ......working on permutation 5## ......working on permutation 6## ......working on permutation 7## ......working on permutation 8## ......working on permutation 9## ......working on permutation 10stats = mp$preservation$Z$ref.A1$inColumnsAlsoPresentIn.A2 stats[order(-stats[,2]),c(1:2)]## moduleSize Zsummary.pres## blue 400 47.700836## brown 400 47.638466## black 188 37.550882## red 210 29.901600## green 325 27.197690## turquoise 400 24.520431## pink 107 17.395392## yellow 400 15.374599## magenta 70 13.565823## grey 400 11.877713## gold 100 5.808094# 一般情况下，“Zsummary”的值越高，数据集之间的模块保守性越好: 5 < Z < 10为中度保守性，Z > 10为高度保守性
好了，今天的内容到此结束，剩下的内容我们下次继续。
后台回复“feng 64”获取代码和数据，我们下期见！
学到了！别人发了CNS的数据，我还能发吗？教你一招，做出新意！（附代码）
终于找到发高分生信文章的财富密码了！（附代码！）
不到一刻钟！高分SCI高逼格套路，我先替你们感受了！
性价比最高的2区纯生信套路，学会够评副高用了！
单细胞分析专栏传送门
大年初一，单细胞分析从nature communication开始
8+单细胞文章从解读到套路总结全部送给你
30分多的单细胞分析文章是怎么做的？
不用R，点点点就能get的单细胞公共数据下载，你“学废”了吗
单细胞都有两个“对象”，你一个都没有
8成生信审稿人必问！单细胞第一张可视化图，你漏了吗？
我从Science中，“偷”到了这个实用的聚类分析技能！真舍不得分享！（附代码）
不骗你！学了这个R包，搞定九成高大上的单细胞分析！包治百病！一包封神！小白赶紧收起来！
Cell的生信分析能有多高级？单细胞研究必备分析！不会做么？简单，我来教你！（附代码）
放大招！单细胞筛分子标志物，让你的课题逼格up! up! up! 你一定要试试！（附详细代码）
来了来了，不就是5+单细胞美图复现吗，我给就是了！
拆篇Nature 又怎样！分析+范文+代码，喂到嘴边的方法，赶紧收藏！
全网首发！单细胞路在何方？北大张泽民课题组万字综述给你指明方向！
单细胞分析绕不过去的坎，三种方法帮你解决！总有一款适合你！
全网首发！万字长文详解Smart-seq2超详细单细胞教程！没别的，就是教你发CNS级别SCI！
相见恨晚！不同数据，不同注释方法，来对号入座！
疯了！想白嫖10+单细胞SCI？只要解决这一步，你真的可以！
碎碎念专栏传送门（完结）
拨开数据迷云雾，得见缺失值真面目
师兄用临床数据一年发五篇文章，原来靠这个诀窍！（建议收藏！）
单个组学玩腻了，来试试20分文章喜欢的多组学联合聚类吧！
顶刊JCO的ROC曲线绘制秘诀，让你的临床研究再上一个level！
风之美图系列传送门（完结）
手把手教你复现一篇CNS级别美图！附代码，建议收藏！
CNS级别的美图复现起来难不难？小白也能快速上手（附代码）（二）
小白课堂 | 学个R包，复现一张CNS级别美图（附代码）
学一个R包给文章加一张CNS级别的美图（四）
学一个R包给文章加一张CNS级别的美图（五）
学一个R包给文章加一张CNS级别的美图（六）
拆和叠？原来你和高阶论文的图片就差这点
可繁可简，可盐可甜，这才是upset的真面目吧！
高分文章不讲xx，竟然用这种图片
荡气回肠一柄剑，看T细胞舍生取义
Cancer Cell的一张美图送给你，新年也发个CC
嘀嘀嘀，美图系列末班车已到站，来收拾最后的图片下车
—
END
—

撰文丨风   风
排版丨四金兄
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~

继续阅读
阅读原文
关键词
网络
数据集
内容
方法
探针