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

原创

庐州月光 2022-06-21 06:12:43 ©著作权

文章标签 数据 for循环颜色渐变 文章分类 OpenStack 云计算

©著作权归作者所有：来自51CTO博客作者庐州月光的原创作品，请联系作者获取转载授权，否则将追究法律责任

欢迎关注”生信修炼手册”!

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

在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"))

可视化的结果如下

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

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

WGCNA将共表达基因与表型数据相关联_for循环_02

上述只是基本用法，适用于样本属于同一组的情况。设想一下，在组间差异非常大的情况下，不同分组条件下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表示, 可以得到如下所示的热图

WGCNA将共表达基因与表型数据相关联_颜色渐变_03