jamesbang

V1

2023/05/11阅读：9主题：雁栖湖

# 🤠 WGCNA | 不止一个组的WGCNA怎么分析嘞！？~（二）（共识网络分析-第二步-构建网络与模块-一步法）

## 2用到的包

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

## 3示例数据

``load("./Consensus-dataInput.RData")``

## 4提取数据集数量

``nSets <-  checkSets(multiExpr)\$nSets``

## 5构建网络与模块识别（一步法）

### 5.1 创建power向量

``powers <-  c(seq(4,10,by=1), seq(12,20, by=2))``

### 5.2 计算power

``powerTables <-  vector(mode = "list", length = nSets) for (set in 1:nSets){   powerTables[[set]] = list(data = pickSoftThreshold(multiExpr[[set]]\$data,                                                      powerVector=powers,                                                      verbose = 2)[[2]])  collectGarbage()}``

### 5.3 可视化参数设置

``  colors = c("black", "red")  plotCols = c(2,5,6,7)   colNames = c("Scale Free Topology Model Fit", "Mean connectivity", "Median connectivity", "Max connectivity")  ylim = matrix(NA, nrow = 2, ncol = 4);  for (set in 1:nSets) {for (col in 1:length(plotCols))  {    ylim[1, col] = min(ylim[1, col], powerTables[[set]]\$data[, plotCols[col]], na.rm = T);     ylim[2, col] = max(ylim[2, col], powerTables[[set]]\$data[, plotCols[col]], na.rm = T);   } }``

### 5.4 可视化

ok，可以正式可视化了。😗

``sizeGrWindow(8, 6)par(mfcol = c(2,2))par(mar = c(4.2, 4.2 , 2.2, 0.5))cex1 = 0.7for (col in 1:length(plotCols)) for (set in 1:nSets){if (set==1){plot(powerTables[[set]]\$data[,1], -sign(powerTables[[set]]\$data[,3])*powerTables[[set]]\$data[,2],xlab="Soft Threshold (power)",ylab=colNames[col],type="n", ylim = ylim[, col],main = colNames[col]);addGrid();}if (col==1){text(powerTables[[set]]\$data[,1], -sign(powerTables[[set]]\$data[,3])*powerTables[[set]]\$data[,2],labels=powers,cex=cex1,col=colors[set]);} elsetext(powerTables[[set]]\$data[,1], powerTables[[set]]\$data[,plotCols[col]],labels=powers,cex=cex1,col=colors[set]);if (col==1){legend("bottomright", legend = setLabels, col = colors, pch = 20) ;} elselegend("topright", legend = setLabels, col = colors, pch = 20) ;}``

### 5.5 识别模块

``net <-  blockwiseConsensusModules(multiExpr, power = 6, minModuleSize = 30,                                   deepSplit = 2, pamRespectsDendro = F,                                   mergeCutHeight = 0.25,                                  # maxBlockSize = 20000,                                  numericLabels = T, minKMEtoStay = 0,                                   saveTOMs = T,                                  verbose = 5)``

### 5.6 查看网络

``names(net)``

### 5.7 转换label

``consMEs <- net\$multiMEsmoduleLabels <- net\$colorsmoduleColors <-  labels2colors(moduleLabels) consTree <-  net\$dendrograms[[1]]``

## 6可视化网络及模块结果

``sizeGrWindow(8,6)plotDendroAndColors(consTree, moduleColors, "Module colors",                     dendroLabels = F, hang = 0.03,                    addGuide = T, guideHang = 0.05,                    main = "Consensus gene dendrogram and module colors")``

## 7Save一下

``save(consMEs, moduleLabels, moduleColors, consTree,      file = "./Consensus-NetworkConstruction-auto.RData")``

📍 往期精彩

##### jamesbang
V1

wx🔍: Grassssss 卷起来了