wgcna_2

knitr::opts_chunk$set(echo = TRUE)

5.导入性状数据

rm(list = ls())
options(stringsAsFactors = F)

load("wgcna-network.Rdata")
datTraits <- read.delim("sample_data/ClinicalTraits.txt", row.names=1, quote="")

head(rownames(datTraits));head(rownames(datExpr))

可以发现两个文件的rownames不一致。

#确保两个文件rownames一致
traitRows = match(rownames(datExpr), rownames(datTraits))
datTraits = datTraits[traitRows, ]

head(rownames(datTraits));head(rownames(datExpr))

6.模块与性状关系

#获得样本数目和基因数目
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
  
#计算模块特征向量 MEs
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
  
#计算模块（MEs）与性状之间的相关性
moduleTraitCor = cor(MEs, datTraits, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

#相关性 heatmap
sizeGrWindow(10,6)

#连接相关性和 pvalue
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)


par(mar = c(6, 8.5, 3, 3))
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")
)

7. 基因与性状关系（GS）& 基因与模块关系（MM）
   计算 ##

modNames = substring(names(MEs), 3)
traitNames = names(datTraits)

7.1计算module membership (MM): 基因（TMP）与模块（MEs）的相关性

geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"))
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership),nSamples)) 
names(geneModuleMembership) = paste("MM", modNames, sep="")
names(MMPvalue) = paste("p.MM", modNames, sep="")

7.2计算Gene Significance (GS): 基因（TMP）与性状的相关性

geneTraitSignificance = as.data.frame(cor(datExpr, datTraits, use = "p"))
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance),nSamples))
names(geneTraitSignificance) = paste("GS.", traitNames, sep="")
names(GSPvalue) = paste("p.GS.", traitNames, sep="")

整合文件

geneInfo<-cbind(geneModuleMembership, MMPvalue, geneTraitSignificance, GSPvalue)
write.table(geneInfo, file = "geneInfo.txt", 
            sep = "\t", 
            quote = F)

save(datExpr, datTraits, MEs, geneModuleMembership, geneTraitSignificance, 
     file = "trait_analysis.RData")