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

# 重新计算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"))

# 样本分为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)];

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