太惊艳 !多数据集联合分析,这样做,能让你的SCI分数更上一层楼!
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系列传送门
在进行相应的分析之前,我们同样先对数据进行预处理,使用的数据和相关代码放在公众号后台,回复关键词即可获取:
rm(list = ls())
"WGCNA") 安装R包 BiocManager::install(
"Hmisc") BiocManager::install(
"flashClust") BiocManager::install(
"qvalue") BiocManager::install(
"dynamicTreeCut") BiocManager::install(
"impute") BiocManager::install(
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.8173
datExprA2[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.8173
datExprB2[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] 21696
length(probesA)
# [1] 20198
genesI和genesA是probesI和probesA对应的基因symbol名称
length(genesI)
# [1] 21696
length(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转换为基因symbol
head(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.7352
head(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 32
datExprA2p <- 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 32
datExprB2g <- 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)
c(2,2)) =
verboseScatterplot(rankExprA1,
rankExprA2,
col = "steelblue",
xlab="Ranked Expression (A1)",
ylab="Ranked Expression (A2)")
verboseScatterplot(rankConnA1,
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,
col = "steelblue",
xlab="Ranked Connectivity (B1)",
ylab="Ranked Connectivity (B2)")
# dev.off()
# 针对这个结果图片,可以给出如下解读:1)如果相关性为正且p值具有显著统计学意义,说明数据集之间的结果可以比较;2)A的相关性和p-values比b的相关性和p-values要好,这是因为两个A数据集是来源于同一个平台运行的。因此我们可以得出结论:来自不同平台的数据集不如来自同一平台的数据集具有可比性,但它们之间仍然具有可比性。
接下来我们使用来自同一平台的A1和A2进行后续的分析。为了计算方便,我们首先选择数据集中表达最多的5000个探针,然后对相同基因名的只保留其中一个探针:
keepGenesExpr <- rank(-rowMeans(datExprA1p)) <= 5000
datExprA1g <- datExprA1p[keepGenesExpr, ]
datExprA2g <- datExprA2p[keepGenesExpr, ]
keepGenesDups <- (collapseRows(datExprA1g,
genesI,
probesI))[[2]]
head(keepGenesDups)
datExprA1g <- datExprA1g[keepGenesDups[,2], ]
datExprA2g <- datExprA2g[keepGenesDups[,2], ]
rownames(datExprA1g) <- rownames(datExprA2g) <- keepGenesDups[,1]
adjacencyA1 <- adjacency(t(datExprA1g),
power = softPower,
type = "signed")
diag(adjacencyA1) = 0
dissTOMA1 <- 1-TOMsimilarity(adjacencyA1,
TOMType = "signed")
geneTreeA1 <- flashClust(as.dist(dissTOMA1),
method = "average")
adjacencyA2 <- adjacency(t(datExprA2g),
power = softPower,
type = "signed")
diag(adjacencyA2) <- 0
dissTOMA2 <- 1-TOMsimilarity(adjacencyA2,
TOMType = "signed")
geneTreeA2 <- flashClust(as.dist(dissTOMA2),
method = "average")
save.image("tutorial.RData")
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)
mColorh,
", 0:3), =
main = "",
dendroLabels=FALSE);
# dev.off()
# 根据结果,这里我们选择dpSplt = 0的结果(此处主要根据自己研究目的进行选择,比较自由,没有统一标准)
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,
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)
c(1,1), mar=c(0, 3, 1, 1) + 0.1, cex=1) =
xlab="",ylab="",main="",sub="")
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 10
stats = 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”获取代码和数据,我们下期见!
单细胞分析专栏传送门
碎碎念专栏传送门(完结)
风之美图系列传送门(完结)
—
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]。