jamesbang

V1

2023/02/15阅读：31主题：雁栖湖

# 🤩 WGCNA | 值得你深入学习的生信分析方法！~（网状分析-第三步-模块与特征分析）

## 2用到的包

``rm(list = ls())library(WGCNA)library(tidyverse)``

## 3示例数据

``load("FemaleLiver-01-dataInput.RData")load("FemaleLiver-02-networkConstruction-auto.RData")``

## 4模块与外部特征关联

### 4.1 量化模块与特征之间的关系

``nGenes <-  ncol(datExpr)nSamples <-  nrow(datExpr)MEs0 <-  moduleEigengenes(datExpr, moduleColors)\$eigengenesMEs <-  orderMEs(MEs0)moduleTraitCor <- cor(MEs, datTraits, use = "p")moduleTraitPvalue <-  corPvalueStudent(moduleTraitCor, nSamples)``

``sizeGrWindow(10,6)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"))``

### 4.2 计算Gene Significance 和 Module Membership

1️⃣ 接着我们将`Gene Significance``GS`） 定义为量化`基因``traits`之间相关性的绝对值。

2️⃣ `Module Membership``MM`）定义为模块的`eigengene`与基因表达谱之间的相关性。

``weight <-  as.data.frame(datTraits\$weight_g);names(weight) <-  "weight"modNames <-  substring(names(MEs), 3)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="")geneTraitSignificance <-  as.data.frame(cor(datExpr, weight, use = "p"))GSPvalue <-  as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))names(geneTraitSignificance) <-  paste("GS.", names(weight), sep="")names(GSPvalue) <-  paste("p.GS.", names(weight), sep="")``

### 4.3 模块内部分析

``module <-  "magenta"column <-  match(module, modNames)moduleGenes <-  moduleColors==modulesizeGrWindow(7, 7)par(mfrow = c(1,1))verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),                   abs(geneTraitSignificance[moduleGenes, 1]),xlab = paste("Module Membership in", module, "module"),ylab = "Gene significance for body weight",main = paste("Module membership vs. gene significance\n"),cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)``

### 4.4 批量输出

``modNames <-  substring(names(MEs), 3)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="")traitNames <- names(datTraits)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="")for (trait in traitNames){  traitColumn = match(trait,traitNames)    for (module2 in modNames){    column = match(module2, modNames)    moduleGenes = moduleColors==module2    if (nrow(geneModuleMembership[moduleGenes,]) > 1){      pdf(file = paste0("./module_", trait, "_", module,".pdf"),          width=7,height=7)            par(mfrow = c(1,1))            verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),                         abs(geneTraitSignificance[moduleGenes, traitColumn]),                         xlab = paste("Module Membership in", module, "module"),                         ylab = paste("Gene significance for ",trait),                         main = paste("Module membership vs. gene significance\n"),                         cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)      dev.off()    }  }}``

## 5结果汇总输出

### 5.1 读入并整理注释文件

``annot <-  read.csv(file = "./FemaleLiver-Data/GeneAnnotation.csv");dim(annot)names(annot)probes <-  names(datExpr)probes2annot <-  match(probes, annot\$substanceBXH)sum(is.na(probes2annot))``

### 5.2 整理并输出结果文件

``geneInfo0 <-  data.frame(substanceBXH = probes,                         geneSymbol = annot\$gene_symbol[probes2annot],                         LocusLinkID = annot\$LocusLinkID[probes2annot],                         moduleColor = moduleColors,                         geneTraitSignificance,                         GSPvalue)modOrder <-  order(-abs(cor(MEs, weight, use = "p")))for (mod in 1:ncol(geneModuleMembership)){oldNames = names(geneInfo0)geneInfo0 = data.frame(geneInfo0, geneModuleMembership[, modOrder[mod]],MMPvalue[, modOrder[mod]]);names(geneInfo0) = c(oldNames, paste("MM.", modNames[modOrder[mod]], sep=""),paste("p.MM.", modNames[modOrder[mod]], sep=""))}geneOrder <-  order(geneInfo0\$moduleColor, -abs(geneInfo0\$GS.weight));geneInfo <-  geneInfo0[geneOrder, ]write.csv(geneInfo, file = "geneInfo.csv")DT::datatable(geneInfo)``

## 6如何引用

📍
```Langfelder, P., Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9, 559 (2008). https://doi.org/10.1186/1471-2105-9-559 ```

📍 往期精彩

##### jamesbang
V1

wx🔍: Grassssss 卷起来了