分享

WGCNA将共表达基因与表型数据相关联

 生信修炼手册 2019-12-24

单纯的共表达基因集合的结果并不能与我们的实验设计相关联,对于识别到的几十个共表达基因集合,一一进行富集分析去挖掘其功能,看上去如此的盲目,没有目的性,所以我们需要对共表达基因集进一步挖掘,常规的做法就是分析其中与性状相关的共表达基因,然后针对这些基因通过富集分析来研究其功能。

在WGCNA中,通过相关性分析将表型数据和共表达基因关联起来。这种方法要求提供每个样本对应的表型数据的值,利用这个值与module的第一主成分值进行相关性分析,根据相关性分析的结果。识别与表型相关联的modules。

表型数据示例如下

sampleweight_glength_cmab_fat
F2_29036.99.92.53
F2_29148.510.72.9
F2_29245.710.41.04
F2_29350.310.90.91

第一列为样本,其他列代表不同的表型,尽量不要有空值,早进行相关性分析时,空值会被剔除,所以太多的空值会影响相关性分析的结果。

在识别modules的过程中,会根据module的第一主成分,即ME值合并modules, 合并之后的modules需要重新计算对应的ME值,然后用ME值与对应的表型数据的值进行相关性分析,代码如下

# 重新计算ME值 MEs0 <- moduleEigengenes(datExpr, moduleColors)$eigengenes MEs <- orderMEs(MEs0) # 计算ME与表型之间的相关系数和p值 moduleTraitCor <- cor(MEs, datTraits, use = "p"); moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples); # 用热图展示相关性的结果 # 每个单元格标记相关系数和p值 textMatrix <- paste(  signif(moduleTraitCor, 2),  "\n(", signif(moduleTraitPvalue, 1), ")", sep = "") dim(textMatrix) <- dim(moduleTraitCor) labeledHeatmap(  Matrix = moduleTraitCor,  xLabels = names(datTraits),  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"))

可视化的结果如下

在该图中,每一行代表一个module, 每一列代表一种表型,每个单元格的颜色由对应的相关系数进行映射,数值从从-1到1,颜色由绿色过渡到白色,然后过渡到红色。这里在运行时,会有一个有趣的小提示,因为红绿色盲的原因,不推荐采用绿色到红色的颜色渐变,建议采用蓝色到红色的渐变,只需要把greenWhiteRed替换为blueWhiteRed即可,效果图如下

上述只是基本用法,适用于样本属于同一组的情况。设想一下,在组间差异非常大的情况下, 不同分组条件下modules与表型数据的相关性结果肯定也会不同,所以对于样本具有不同分组的数据,需要不同分组分开分析,WGCNA当然也支持这样的分析,不同分组的表达量保存在不同文件中,然后构建一个list对象,长度和分组个数相同,每个元素对应一个分组条件下的表达量数据

# 样本分为male和female两组,分开读取 femData = read.csv("LiverFemale3600.csv") maleData = read.csv("LiverMale3600.csv") # 分组个数 nSets = 2; setLabels = c("Female liver", "Male liver") shortLabels = c("Female", "Male") # 构建总的表达量,长度为nSets的list multiExpr = vector(mode = "list", length = nSets) # 每个元素对应一个分组下的表达量数据 multiExpr[[1]] = list(data = as.data.frame(t(femData[-c(1:8)]))); names(multiExpr[[1]]$data) = femData$substanceBXH; rownames(multiExpr[[1]]$data) = names(femData)[-c(1:8)]; multiExpr[[2]] = list(data = as.data.frame(t(maleData[-c(1:8)]))); names(multiExpr[[2]]$data) = maleData$substanceBXH; rownames(multiExpr[[2]]$data) = names(maleData)[-c(1:8)];

通过上述方式合并不同分组对应的表达量数据,然后一起识别modules, 不考虑分组,所有样本一起识别到的module称为consensus modules, 在后续与表型数据进行相关性分析时,通过循环,对每一组单独进行分析,代码如下

moduleTraitCor = list() moduleTraitPvalue = list() for (set in 1:nSets) {  moduleTraitCor[[set]] = cor(    consMEs[[set]]$data,    Traits[[set]]$data,    use = "p") }

for循环中的代码和一开始提到的基本用法一致,所以对于每个group, 都可以产生上述的相关性结果的热图,除此之外,还可以分析在不同分组中,共表达的趋势是否一致,如果表达趋势不同,一个为正相关,一个为父相关,则用NA表示, 可以得到如下所示的热图

在该图中,只有在两组中共表达趋势相同的modules才会有颜色填充。

所谓的与表型数据关联,其实就是一个相关性分析,最后可以根据相关性的分析结果,筛选与某种表型显著相关的modules。更多细节请参考官方文档。

·end·

    转藏 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多