WGCNA将共表达基因与表型数据相关联
欢迎关注”生信修炼手册”!
单纯的共表达基因集合的结果并不能与我们的实验设计相关联，对于识别到的几十个共表达基因集合，一一进行富集分析去挖掘其功能，看上去如此的盲目，没有目的性，所以我们需要对共表达基因集进一步挖掘，常规的做法就是分析其中与性状相关的共表达基因，然后针对这些基因通过富集分析来研究其功能。
在WGCNA中，通过相关性分析将表型数据和共表达基因关联起来。这种方法要求提供每个样本对应的表型数据的值，利用这个值与module的第一主成分值进行相关性分析，根据相关性分析的结果。识别与表型相关联的modules。
表型数据示例如下
sample | weight_g | length_cm | ab_fat |
|---|---|---|---|
F2_290 | 36.9 | 9.9 | 2.53 |
F2_291 | 48.5 | 10.7 | 2.9 |
F2_292 | 45.7 | 10.4 | 1.04 |
F2_293 | 50.3 | 10.9 | 0.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·
—如果喜欢，快分享给你的朋友们吧—
腾讯云开发者
