用老分析还被夸做出创新点?这个百搭的分析,让你的课题idea源源不断!
WGCNA应用实战三
大家好,我是风风。今天我们要继续我们的WGCNA系列,这也是我们WGCNA系列的最后一次推文了。在上次的推文中,我们一起看了多数据集的联合分析——由Jeremy Miller提出的对来自多个芯片数据集的数据进行荟萃分析,在最后利用多个数据集,对识别的模块进行了Preservation分析,今天我们接着上次的分析,看看保守性分析后面的内容应该怎么做。
WGCNA系列传送门
kME是一个在WGCNA中常用的元素值,因为它可以用来测量每个基因和每个ME(Module Epigenes)之间的相关性,因此即使是最初没有分配到一个模块的基因也可以被在网络间的被进行比较,我们来看看具体的操作:
rm(list = ls())
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[1:3, 1:3]
datExprA2[1:3, 1:3]
datExprB1[1:3, 1:3]
datExprB2[1:3, 1:3]
length(probesI)
length(probesA)
length(genesI)
length(genesA)
datExprA1g[1:3, 1:3]
datExprA2g[1:3, 1:3]
ME_1A[1:3, 1:3]
head(modulesA1)
geneModuleMembership1 <- signedKME(t(datExprA1g), ME_1A)
colnames(geneModuleMembership1) <- paste("PC",colorsA1,".cor",sep="");
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))))
write.csv(kMEtable1, "kMEtable1.csv", row.names = FALSE)
PCs2A <- moduleEigengenes(t(datExprA2g),
colors = modulesA1)
ME_2A <- PCs2A$eigengenes
geneModuleMembership2 <- signedKME(t(datExprA2g),
ME_2A)
colnames(geneModuleMembership1) <- paste("PC", colorsA1, ".cor", sep="");
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)
write.csv(kMEtable2, "kMEtable2.csv", row.names=FALSE)
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")
}
topGenesKME <- NULL
for (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
softPower <- 10
source("tutorialFunctions.R")
for (co in colorsA1[colorsA1!="grey"])
visantPrepOverall(modulesA1, co, t(datExprA1g), rownames(datExprA1g), 500, softPower, TRUE)
for (co in colorsA1[colorsA1!="grey"])
visantPrepOverall(modulesA1, co, t(datExprA2g), rownames(datExprA2g), 500, softPower, TRUE)
除了上面说的使用外部来源来注释模块之外,在模块分析和比较(以及评估差异表达)中几乎总是可以使用某种表型信息(和我们前面的基础推文一致)。例如,在评估药物治疗效果的实验中,样本至少会被指定为治疗组或对照组。同样,在评估对照和疾病之间变化的实验中,将提供有关疾病状态的表型信息。此外,还经常提供诸如年龄、性别和大脑区域等基本信息。为了演示如何执行这样的分析,我们首先需要有表型信息:
# 构建表型信息(如果有自己的表型信息,直接读取输入即可)
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.01651317
var <- 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中更为显著。
到目前为止,在荟萃分析中,我们已经比较了这样的网络:其中来自一个网络的基因的模块标识符被指定为第二个网络中相同基因的模块标识符。然而,另一种常见情况是我们经常得到两个具有不同模块标识基因的网络。在这种情况下,如何比较两个网络的模块呢?我们可以使用下面的方法:
datExprB2g <- (collapseRows(datExprB2,
genesA,
probesA))[[1]]
GeneAB <- sort(intersect(rownames(datExprB2g), Gene))
head(GeneAB)
datExprB2g <- datExprB2g[GeneAB,]
adjacencyB2 <- adjacency(t(datExprB2g),power=softPower,type="signed");
diag(adjacencyB2) <- 0
dissTOMB2 <- 1-TOMsimilarity(adjacencyB2, TOMType="signed")
geneTreeB2 <- flashClust(as.dist(dissTOMB2), method="average")
mColorh <- NULL
for (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));
}
plotDendroAndColors(geneTreeB2, mColorh, paste("dpSplt =",0:3), main = "",dendroLabels=FALSE);
modulesB2 <- mColorh[,1]
colorsB2 <- names(table(modulesB2))
modulesB2_new = matchModules(Gene, modulesA1, GeneAB, modulesB2)
enrichmentsB2A1 = userListEnrichment(GeneAB,modulesB2,"kMEtable1.csv","A1","enrichmentB2_A1.csv")
enrichmentsB2A1$sigOverlaps
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”获取代码和数据,我们下期见!
单细胞分析专栏传送门
碎碎念专栏传送门(完结)
风之美图系列传送门(完结)
—
END—
撰文丨风 风
排版丨四金兄
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
关键词
基因
方法
内容
网络
结果
最新评论
推荐文章
作者最新文章
你可能感兴趣的文章
Copyright Disclaimer: The copyright of contents (including texts, images, videos and audios) posted above belong to the User who shared or the third-party website which the User shared from. If you found your copyright have been infringed, please send a DMCA takedown notice to [email protected]. For more detail of the source, please click on the button "Read Original Post" below. For other communications, please send to [email protected].
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。