用老分析还被夸做出创新点？这个百搭的分析，让你的课题idea源源不断！

WGCNA应用实战三

大家好，我是风风。今天我们要继续我们的WGCNA系列，这也是我们WGCNA系列的最后一次推文了。在上次的推文中，我们一起看了多数据集的联合分析——由Jeremy Miller提出的对来自多个芯片数据集的数据进行荟萃分析，在最后利用多个数据集，对识别的模块进行了Preservation分析，今天我们接着上次的分析，看看保守性分析后面的内容应该怎么做。

WGCNA系列传送门

仅凭一个分析就能发3.5+SCI！万字长文教你彻底学会WGCNA！不分享说不过去吧？

比另一半更爱你！被你们问的最多的生信分析！一个分析就是一篇3+SCI！高级！给力！

被问爆的WGCNA教程，高分秘籍，有手就会！（附代码）

决定了！把免疫做到新高度就用这个分析！|（附代码 )

太惊艳！多数据集联合分析，这样做，能让你的SCI分数更上一层楼！

识别Module Membership（MM/kME）

kME是一个在WGCNA中常用的元素值，因为它可以用来测量每个基因和每个ME（Module Epigenes）之间的相关性，因此即使是最初没有分配到一个模块的基因也可以被在网络间的被进行比较，我们来看看具体的操作：

rm(list = ls())# BiocManager::install("WGCNA") 安装R包# BiocManager::install("Hmisc") # BiocManager::install("flashClust") # BiocManager::install("qvalue") # BiocManager::install("dynamicTreeCut") # BiocManager::install("impute") # 加载R包options(stringsAsFactors = FALSE)library(WGCNA)library(Hmisc)library(flashClust)library(qvalue)library(dynamicTreeCut)library(impute)library(ggplot2)source("collapseRows_NEW.R") # 加载数据load("metaAnalysisData.RData")# 加载前面运行的结果load("step.Rda")# 查看数据# 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# datExprA1g、datExprA2g、ME_1A、colorsA1、modulesA1都是前面运行的结果datExprA1g[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## 11-Sep 1226.443 1064.261 1178.000## 15-Sep 2907.731 1887.192 2986.067## 2-Mar 1047.459 1906.995 1072.987datExprA2g[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## 11-Sep 1720.986 1465.163 1896.664## 15-Sep 3578.331 2488.593 2088.214## 2-Mar 1097.692 1780.784 1716.164ME_1A[1:3, 1:3]## MEblack MEblue MEbrown## C.1385.CA3.M.81.4.O 0.12461289 0.1949176 0.2090233## C.3795.CA1.F.72.4.V -0.14704184 -0.2249308 -0.1999316## C.1761.CA1b.M.86.4.O 0.07249517 0.2964738 0.2158585head(modulesA1)## [1] "yellow""pink""turquoise""green""black""red"# 首先是A1数据集# 计算kME值geneModuleMembership1 <- signedKME(t(datExprA1g), ME_1A) colnames(geneModuleMembership1) <- paste("PC",colorsA1,".cor",sep=""); # 计算p值MMPvalue1 <- corPvalueStudent(as.matrix(geneModuleMembership1), dim(datExprA1g)[[2]])colnames(MMPvalue1) <- paste("PC", colorsA1,".pval", sep="")Gene <- rownames(datExprA1g) # 构建数据框kMEtable1 <- cbind(Gene, Gene, modulesA1) for (i in1:length(colorsA1)) kMEtable1 <- cbind(kMEtable1, geneModuleMembership1[,i], MMPvalue1[,i]) colnames(kMEtable1) <- c("PSID","Gene","Module", sort(c(colnames(geneModuleMembership1), colnames(MMPvalue1)))) # 保存A1的结果write.csv(kMEtable1, "kMEtable1.csv", row.names = FALSE)# 接着是A2数据集，使用A1的模块分配确定kME值# 计算A2的MEsPCs2A <- moduleEigengenes(t(datExprA2g), colors = modulesA1) ME_2A <- PCs2A$eigengenes # 计算相关性geneModuleMembership2 <- signedKME(t(datExprA2g), ME_2A) colnames(geneModuleMembership1) <- paste("PC", colorsA1, ".cor", sep=""); # 计算P值MMPvalue2 <- corPvalueStudent(as.matrix(geneModuleMembership2),dim(datExprA2g)[[2]]); colnames(MMPvalue2) <- paste("PC",colorsA1,".pval",sep=""); kMEtable2 <- cbind(Gene,Gene,modulesA1) for (i in1:length(colorsA1)) kMEtable2 <- cbind(kMEtable2, geneModuleMembership2[,i], MMPvalue2[,i]) colnames(kMEtable2) <- colnames(kMEtable1) # 保存A2的结果 write.csv(kMEtable2, "kMEtable2.csv", row.names=FALSE)# 现在我们已经有了两个网络的KME值，我们可以来比较下生成的网络

# 首先是绘制A1中每个基因的KME值与A2中每个基因的相应KME值，具有高度相关性的点的模块是高度稳定的模块：

# 利用所有基因创建这些曲线图的结果图片# pdf("all_kMEtable2_vs_kMEtable1.pdf", height = 8, width = 8) for (c in1:length(colorsA1)){ verboseScatterplot(geneModuleMembership2[,c], geneModuleMembership1[,c], main=colorsA1[c], xlab="kME in A2", ylab="kME in A1") }

# dev.off() # 利用最初分配给给定模块的基因子集创建这些曲线图的结果图片# pdf("inModule_kMEtable2_vs_kMEtable1.pdf", height = 8, width = 8) for (cin1:length(colorsA1)){ inMod = modulesA1== colorsA1[c] verboseScatterplot(geneModuleMembership2[inMod,c],geneModuleMembership1[inMod,c],main=colorsA1[c], xlab="kME in A2",ylab="kME in A1") }

# dev.off() # save.image("tutorial.RData") # 这样我们就绘制了所有模块的散点图

# 一般我们认为，使用所有基因的散点图可以包含所有正相关和负相关的基因，但通常也包含大量干扰信息。而仅使用模块内基因是评估HUB基因保守性的直观方法：如果这些基因显示集合间的相关性，那么位于散点图右上角的基因很可能是数据集之间的公共HUB基因。(HUB基因是与MES和模块内高度连通性显著相关的基因

# 我们还可以通过确定哪些基因在两个网络中都有极高的KME值，来确定哪些基因是这两个网络中的关键HUB基因：

topGenesKME <- NULLfor (c in1:length(colorsA1)){ kMErank1 <- rank(-geneModuleMembership1[,c]) kMErank2 <- rank(-geneModuleMembership2[,c]) maxKMErank <- rank(apply(cbind(kMErank1,kMErank2+.00001), 1, max)) topGenesKME <- cbind(topGenesKME,Gene[maxKMErank<=10]) }colnames(topGenesKME) <- colorsA1 topGenesKME## black blue brown green grey magenta ## [1,] "BSN""ACTR10""CADPS""AGT""ATP13A2""AHSA1"## [2,] "C19orf4""ATP5C1""LRFN5""CYBRD1""DNM1""DNAJA1"## [3,] "CACNA1E""ATP5O""MAP7D2""DPYSL3""EEF1A2""DNAJB1"## [4,] "DNM1""COX6B1""NAPB""FERMT2""HSPA12A""DNAJB6"## [5,] "EEF1A2""KIAA1279""NECAP1""HEPACAM""PACSIN1""FKBP4"## [6,] "ELOVL4""LOC646630""NSF""IL17D""PDXP""HSPA1A"## [7,] "FLJ33996""PRDX3""PGM2L1""LOC219854""PRRT3""HSPA1B"## [8,] "LOC649970""TBCA""PREPL""MAPRE1""RAB6B""HSPB1"## [9,] "LPHN1""TMEM85""SLC9A6""NTRK2""SCAMP5""HSPD1"## [10,] "THY1""TXN""WDR7""RHOBTB3""TMEM130""STIP1"## pink red turquoise yellow ## [1,] "15-Sep""CALN1""ANKRD30B""C5orf4"## [2,] "C1orf19""COL4A3BP""C9orf130""CLIC4"## [3,] "DAZAP2""FAM130A2""CEP27""CNTN2"## [4,] "DYNC1I2""GFOD1""EIF2C2""CREB5"## [5,] "ESD""NUDCD3""IL10""CTNNA1"## [6,] "GABARAPL2""PODXL2""LOC653086""HRASLS3"## [7,] "KHDRBS1""PPFIA4""SDHALP1""LASS2"## [8,] "POLR2G""RUSC2""SLC4A5""PLEKHH1"## [9,] "SF3B14""SLC4A3""ZNF652""RDX"## [10,] "SFT2D1""SNPH""ZNF786""VAMP3"

# 这样我们就识别出了模块中的关键HUB基因，这些基因代表了两个网络中基于KME的每个模块的前10个基因

# 我们还可以输出每一个模块的结果用于cytoscape构建可视化网络softPower <- 10# 设置软阈值source("tutorialFunctions.R") for (co in colorsA1[colorsA1!="grey"]) visantPrepOverall(modulesA1, co, t(datExprA1g), rownames(datExprA1g), 500, softPower, TRUE)## black_connectivityOverall.csv written.## 25 0.49 80## black_visantOverall.csv written.## blue_connectivityOverall.csv written.## 201 0.4075 10## blue_visantOverall.csv written.## brown_connectivityOverall.csv written.## 201 0.5325 28## brown_visantOverall.csv written.## green_connectivityOverall.csv written.## 4 0.28 54## green_visantOverall.csv written.## magenta_connectivityOverall.csv written.## 19 0.43 14## magenta_visantOverall.csv written.## pink_connectivityOverall.csv written.## 9 0.33 54## pink_visantOverall.csv written.## red_connectivityOverall.csv written.## 4 0.28 64## red_visantOverall.csv written.## turquoise_connectivityOverall.csv written.## 301 0.45875 2## turquoise_visantOverall.csv written.## yellow_connectivityOverall.csv written.## 101 0.345 32## yellow_visantOverall.csv written.

# 文件夹下将会产生很多文件，每一个模块都有两个关键文件：1) 模块名_connectivityOverall.csv：这个文件包含给定模块中所有基因的列表，这些基因从模块内连接性(KIN)最高到最低排序，以及它们的平均表达。因此，排在首位的基因应该是中枢HUB基因。2) 模块名_visantOverall.csv：该文件包含cytoscape可视化网络需要的所有输入，也可以使用文档推荐的visant软件，操作也很简单，直接将文件导入软件之中即可

# 接着同样对第二个数据集运行相同的操作：for (co in colorsA1[colorsA1!="grey"]) visantPrepOverall(modulesA1, co, t(datExprA2g), rownames(datExprA2g), 500, softPower, TRUE)## black_connectivityOverall.csv written.## 30 0.54 100## black_visantOverall.csv written.## blue_connectivityOverall.csv written.## 201 0.3525 30## blue_visantOverall.csv written.## brown_connectivityOverall.csv written.## 201 0.5425 42## brown_visantOverall.csv written.## green_connectivityOverall.csv written.## 2 0.26 36## green_visantOverall.csv written.## magenta_connectivityOverall.csv written.## 15 0.39 56## magenta_visantOverall.csv written.## pink_connectivityOverall.csv written.## 3 0.27 44## pink_visantOverall.csv written.## red_connectivityOverall.csv written.## 101 0.355 50## red_visantOverall.csv written.## turquoise_connectivityOverall.csv written.## 24 0.48 64## turquoise_visantOverall.csv written.## yellow_connectivityOverall.csv written.## 101 0.195 54## yellow_visantOverall.csv written.

# 此外，有了模块中的特定基因后，我们还可以进行相应的富集分析的注释，可以利用下面这些网站对模块基因进行注释，分析模块基因参与的生物学功能

# ChiliBot: http://www.chilibot.net/ # WebGestalt: http://bioinfo.vanderbilt.edu/webgestalt/ # Ingenuity: http://www.ingenuity.com/ # Galaxy: http://main.g2.bx.psu.edu/ # GSEA: http://www.broadinstitute.org/gsea/index.jsp # UGET: http://genome.ucla.edu/projects/UGET # STRING: http://string.embl.de/

联合表型信息确定关键表达基因和模块

除了上面说的使用外部来源来注释模块之外，在模块分析和比较(以及评估差异表达)中几乎总是可以使用某种表型信息（和我们前面的基础推文一致）。例如，在评估药物治疗效果的实验中，样本至少会被指定为治疗组或对照组。同样，在评估对照和疾病之间变化的实验中，将提供有关疾病状态的表型信息。此外，还经常提供诸如年龄、性别和大脑区域等基本信息。为了演示如何执行这样的分析，我们首先需要有表型信息：

# 构建表型信息（如果有自己的表型信息，直接读取输入即可）region = rep("CA1",32); region[c(1,4,6,11,12,15,16,17,22,24,25,26,28,29,31,32)] = "CA3"age = c(81,72,86,90,88,90,90,74,83,73,73,70,85,85,75,90,72,70,90,84,75,85,80, 86,85,84,81,88,80,90,83,74) # 接着我们将在数据集A1中找到在年龄和大脑区域中表现出最高差异表达的基因：# 联合表型寻找差异基因# region (Top 3)var <- list(region == "CA1", region == "CA3") datReg <- t(apply(datExprA1g, 1, t.test.l)) colnames(datReg) <- c("MeanCA1","MeanCA3","SD_CA1","SD_CA3","PvalRegion") datReg[order(datReg[,5])[1:3],]## MeanCA1 MeanCA3 SD_CA1 SD_CA3 PvalRegion## NRIP3 10509.884 22481.709 2951.5061 6239.7370 6.739892e-07## SAT1 4656.685 2964.674 885.3605 658.5044 1.335158e-06## KCNH3 3814.777 1840.010 1049.7946 749.2281 1.493591e-06# age (Top 3)var <- age datAge <- t(apply(datExprA1g,1,cor.test.l)) colnames(datAge) <- c("CorrAge","PvalAge") datAge[order(datAge[,2])[1:3],]## CorrAge PvalAge## DDX42 0.7116889 4.947339e-06## SCD 0.6961764 9.660017e-06## KIAA0907 -0.6818583 1.728760e-05

# 我们可以使用类似的策略来查找所有与区域或年龄相关的模块，使用模块eigengene作为每个模块中表达水平基因的代表

var <- list(region == "CA1", region == "CA3") datRegM <- t(apply(t(ME_1A), 1, t.test.l)) colnames(datRegM) <- c("MeanCA1","MeanCA3","SD_CA1","SD_CA3","PvalRegion") datRegM[datRegM[,5]<0.02,]## MeanCA1 MeanCA3 SD_CA1 SD_CA3 PvalRegion## MEblack -0.07880119 0.07880119 0.1250271 0.1943906 0.01136588## MEgreen 0.07462621 -0.07462621 0.1864653 0.1414803 0.01651317var <- age datAgeM <- t(apply(t(ME_1A),1,cor.test.l)) colnames(datAgeM) <- c("CorrAge","PvalAge") datAgeM[datAgeM[,2]<0.02,]## CorrAge PvalAge## MEmagenta -0.4735782 0.006185179## MEpink 0.4747031 0.006047230

# 总的来说，我们发现了几个与区域和年龄相关的模块：绿色模块中CA1高表达，粉红色模块中CA3高表达，黑色模块显示随着年龄的增加表达增加，洋红和红色显示随着年龄的增加表达减少(p < 0.02)。

# 接着我们可以使用WGCNA库中的函数来可视化这些结果# pdf("RegionAgePlots.pdf", width = 16, height = 4) par(mfrow=c(1,4)) verboseBoxplot(as.numeric(datExprA1g["NRIP3",]), region, notch = FALSE, main="NRIP3 expression -", las=2, xlab="Region", ylab="") verboseScatterplot(age, as.numeric(datExprA1g["DDX42",]), main="DDX42 expression -", las=2, abline=TRUE, xlab="Age", ylab="") verboseBoxplot(as.numeric(ME_1A[,"MEgreen"]), region, notch = FALSE, main="Green ME expr. -", las=2, xlab="Region", ylab="") verboseScatterplot(age, as.numeric(ME_1A[,"MEmagenta"]), main="Magenta ME expr. -", las=2, abline=TRUE, xlab="Age", ylab="")

# dev.off()

# 我们也可以使用这些图片在不同的数据集上直观地比较基因或模块与表型的关系。要做到这一点，我们需要两个数据集中相同特征的信息，然后我们只需画出可比较的图表即可：

region2 <- rep("CA1",31) region2[c(1,4,7,8,9,11,15,18,20,21,22,26,28,29,30)] = "CA3"var <- list(region2=="CA1", region2=="CA3") datReg2 <- t(apply(datExprA2g,1,t.test.l)) colnames(datReg2) <- c("MeanCA1","MeanCA3","SD_CA1","SD_CA3","PvalRegion") datRegM2 <- t(apply(t(ME_2A),1,t.test.l)) colnames(datRegM2) <- c("MeanCA1","MeanCA3","SD_CA1","SD_CA3","PvalRegion") # pdf("RegionAgePlots12.pdf", width = 16, height = 4) par(mfrow=c(1,4)) verboseBoxplot(as.numeric(datExprA1g["NRIP3",]), region, notch = FALSE, main="NRIP3 expression (A1) -", las=2, xlab="Region (A1)", ylab="") verboseBoxplot(as.numeric(datExprA2g["NRIP3",]), region2, notch = FALSE, main="NRIP3 expression (A2) -", las=2, xlab="Region (A2)", ylab="") verboseBoxplot(as.numeric(ME_1A[,"MEgreen"]), region, notch = FALSE, main="Green ME expr. (A1) -", las=2, xlab="Region (A1)", ylab="") verboseBoxplot(as.numeric(ME_2A[,"MEgreen"]), region2, notch = FALSE, main="Green ME expr. (A2) -", las=2, xlab="Region (A2)", ylab="")

# dev.off()

# 从结果中我们发现NRIP3和绿色模块在A1和A2中关于区域表型的趋势比较相似，并且在数据集A1中更为显著。

比较具有不同模块定义的网络

到目前为止，在荟萃分析中，我们已经比较了这样的网络：其中来自一个网络的基因的模块标识符被指定为第二个网络中相同基因的模块标识符。然而，另一种常见情况是我们经常得到两个具有不同模块标识基因的网络。在这种情况下，如何比较两个网络的模块呢？我们可以使用下面的方法：

# 利用datExprB2构建一个拓扑异构网络datExprB2g <- (collapseRows(datExprB2, genesA, probesA))[[1]]## Warning: row names of input data and probes not identical...## ... Attempting to proceed anyway. Check results carefully.GeneAB <- sort(intersect(rownames(datExprB2g), Gene)) head(GeneAB)## [1] "A2BP1""A2M""AARS""AASDHPPT""AATF""ABAT"datExprB2g <- datExprB2g[GeneAB,] adjacencyB2 <- adjacency(t(datExprB2g),power=softPower,type="signed"); diag(adjacencyB2) <- 0dissTOMB2 <- 1-TOMsimilarity(adjacencyB2, TOMType="signed")## ..connectivity..## ..matrix multiplication (system BLAS)..## ..normalization..## ..done.geneTreeB2 <- flashClust(as.dist(dissTOMB2), method="average") mColorh <- NULLfor (ds in0:3){ tree <- cutreeHybrid(dendro = geneTreeB2, pamStage=FALSE, minClusterSize = (30-3*ds), cutHeight = 0.99, deepSplit = ds, distM = dissTOMB2) mColorh <- cbind(mColorh,labels2colors(tree$labels)); }## ..done.## ..done.## ..done.## ..done.# pdf("Module_choices_B2.pdf", height = 10, width = 25)plotDendroAndColors(geneTreeB2, mColorh, paste("dpSplt =",0:3), main = "",dendroLabels=FALSE);

# dev.off() modulesB2 <- mColorh[,1] colorsB2 <- names(table(modulesB2))

# 我们现在可以检查模块重叠：一旦我们有了模块定义，我们就可以检查网络A1中的哪些模块包含大量与B2中的模块重叠的基因。然后，这些将是两个网络中彼此对应的模块，并且来自一个网络的标签可以被重新映射，使得对应的模块被分配相同的颜色：

modulesB2_new = matchModules(Gene, modulesA1, GeneAB, modulesB2)## Old - New labels: black - greenyellow## Old - New labels: blue - tan## Old - New labels: brown - brown## Old - New labels: green - turquoise## Old - New labels: grey - grey## Old - New labels: magenta - green## Old - New labels: pink - salmon## Old - New labels: purple - blue## Old - New labels: red - pink## Old - New labels: turquoise - yellow## Old - New labels: yellow - cyanenrichmentsB2A1 = userListEnrichment(GeneAB,modulesB2,"kMEtable1.csv","A1","enrichmentB2_A1.csv") enrichmentsB2A1$sigOverlaps## InputCategories UserDefinedCategories CorrectedPvalues ## pOut "brown""A1_brown""3.27182333304024e-54"## pOut "turquoise""A1_yellow""2.02712525903512e-44"## pOut "blue""A1_brown""3.62384230587342e-36"## pOut "green""A1_turquoise""4.51672464729032e-27"## pOut "purple""A1_blue""7.98370066904828e-24"## pOut "magenta""A1_green""4.25773848289337e-23"## pOut "blue""A1_blue""1.4748170949109e-22"## pOut "grey""A1_grey""6.93179215442581e-19"## pOut "pink""A1_yellow""3.99346048959251e-16"## pOut "yellow""A1_turquoise""5.65304886894366e-16"## pOut "brown""A1_black""1.7292566856617e-13"## pOut "black""A1_turquoise""1.83080307092811e-13"## pOut "red""A1_pink""7.21374392668378e-11"## pOut "grey""A1_red""2.03648220919647e-07"## pOut "red""A1_blue""9.74459904505315e-07"## pOut "turquoise""A1_pink""4.65400459255205e-05"## pOut "red""A1_yellow""5.61984797719157e-05"## pOut "pink""A1_turquoise""0.000130701279529847"

# 总的来说，这两个网络之间有很多明显重叠的模块。这可以通过使用以下代码一次绘制两个树状图来实现

# pdf("A1_B2_dendrogram_plots.pdf", height = 5, width = 15) # plotDendroAndColors(geneTreeA1, modulesA1, "A1", dendroLabels=FALSE, hang=0.03, # addGuide=TRUE, guideHang=0.05, main="A1 dendrogram") plotDendroAndColors(geneTreeB2, modulesB2_new, "B2_new", dendroLabels=FALSE, hang=0.03, addGuide=TRUE, guideHang=0.05, main="B2 dendrogram")

# dev.off()

# 这种方法适用于当我们分别用两个或者多个数据集构建WGCNA网络之后，对模块的稳定性以及模块进行比较的分析。比较模块注释是确定网络间模块一致性的另一种方法。还可以通过对两个数据集的模块进行富集分析，然后查看哪些模块共享了相同注释来实现。例如，如果A1中的红色模块和B2中的棕色模块都有线粒体基因的过度表示，那我们就可能会认为这两个模块是重叠的模块，即使它们没有大量共同的基因。同样，如果A1和B2中的黄色模块都富含星形胶质细胞基因，这就进一步证明了这些基因具有共同的功能。

好了，今天的内容到此结束，这两次的内容是我个人非常喜欢的WGCNA的内容，我个人觉得，Jeremy Miller的这种方法，不仅提供了一种多数据集之间联合分析的方法，还提供了一种更见稳健地寻找hub基因并进行注释获取相应通路和表型的方法，这种基于多数据集的方法的分析结果更加稳定，实验验证得到阳性结果的机率也更大，此外，这种方法也避免了不同平台数据或者不同来源数据合并后需要进行批间差矫正的问题，这一点我个人觉得非常重要！好了，我们WGCNA系列到此全部结束，下一次我们将开启新的系列。

后台回复“feng 65”获取代码和数据，我们下期见！