jamesbang

V1

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

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

## 1写在前面

OK，今天的教程是超大数据集的网络构建与模块识别，当然如果你的电脑或者服务器够好的话，可以跳过这期。🥳

## 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 = TRUE);ylim[2, col] = max(ylim[2, col], powerTables[[set]]\$data[, plotCols[col]], na.rm = TRUE);}}``

### 5.4 可视化一下吧

``sizeGrWindow(8, 6)par(mfcol = c(2,2));par(mar = c(4.2, 4.2 , 2.2, 0.5))cex1 = 0.7;for (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) ;}``

## 6Block-wise的方式构建网络

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

## 7对比一下

### 7.1 一步法结果

``load(file = "./Consensus-NetworkConstruction-auto.RData")bwLabels <-  matchLabels(bnet\$colors, moduleLabels, pThreshold = 1e-7)bwColors <-  labels2colors(bwLabels)table(bwLabels)``

### 7.2 block分别可视化

``sizeGrWindow(12,6)layout(matrix(c(1:4), 2, 2), heights = c(0.8, 0.2), widths = c(1,1))nBlocks = length(bnet\$dendrograms)for (block in 1:nBlocks){  plotDendroAndColors(bnet\$dendrograms[[block]],                       moduleColors[bnet\$blockGenes[[block]]],                      "Module colors",                      main = paste("Gene dendrogram and module colors in block", block),                      dendroLabels = F, hang = 0.03,                      addGuide = T, guideHang = 0.05,                      setLayout = F)  }``

### 7.3 对比一下结果

``sizeGrWindow(12,9)plotDendroAndColors(consTree,                    cbind(moduleColors, bwColors),                    c("Single block", "Blockwise"),                    dendroLabels = F, hang = 0.03,                    addGuide = T, guideHang = 0.05,                    main = "Single block consensus gene dendrogram and module colors")``

📍 往期精彩

##### jamesbang
V1

wx🔍: Grassssss 卷起来了