jamesbang

V1

2023/05/17阅读：12主题：雁栖湖

# 🤠 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 = 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) ;}``

## 6网络邻接的计算

``softPower <-  6adjacencies <-  array(0, dim = c(nSets, nGenes, nGenes))for (set in 1:nSets){  adjacencies[set, , ] <-  abs(cor(multiExpr[[set]]\$data, use = "p"))^softPower  }``

## 7拓扑重叠的计算

``TOM <-  array(0, dim = c(nSets, nGenes, nGenes))for (set in 1:nSets){  TOM[set, , ] <-  TOMsimilarity(adjacencies[set, , ])  }``

## 8给TOM做个scale

### 8.1 scale一下

``scaleP = 0.95set.seed(12345)nSamples = as.integer(1/(1-scaleP) * 1000)scaleSample = sample(nGenes*(nGenes-1)/2, size = nSamples)TOMScalingSamples = list()scaleQuant = rep(1, nSets)scalePowers = rep(1, nSets)for (set in 1:nSets){  TOMScalingSamples[[set]] = as.dist(TOM[set, , ])[scaleSample]  scaleQuant[set] = quantile(TOMScalingSamples[[set]],                             probs = scaleP, type = 8)if (set>1){  scalePowers[set] = log(scaleQuant[1])/log(scaleQuant[set])  TOM[set, ,] = TOM[set, ,]^scalePowers[set]}  }``

### 8.2 可视化scale效果

``scaledTOMSamples = list()for (set in 1:nSets){scaledTOMSamples[[set]] = TOMScalingSamples[[set]]^scalePowers[set]}sizeGrWindow(6,6)qqUnscaled = qqplot(TOMScalingSamples[[1]], TOMScalingSamples[[2]], plot.it = T, cex = 0.6,                    xlab = paste("TOM in", setLabels[1]), ylab = paste("TOM in", setLabels[2]),                    main = "Q-Q plot of TOM", pch = 20)qqScaled = qqplot(scaledTOMSamples[[1]], scaledTOMSamples[[2]], plot.it = FALSE)points(qqScaled\$x, qqScaled\$y, col = "red", cex = 0.6, pch = 20);abline(a=0, b=1, col = "blue")legend("topleft", legend = c("Unscaled TOM", "Scaled TOM"), pch = 20, col = c("black", "red"))``

## 9共识拓扑重叠的计算

``consensusTOM = pmin(TOM[1, , ], TOM[2, , ])``

## 10聚类与模块识别

### 10.1 计算模块

``consTree <-  hclust(as.dist(1-consensusTOM), method = "average")minModuleSize <-  30unmergedLabels <-  cutreeDynamic(dendro = consTree, distM = 1-consensusTOM,                                 deepSplit = 2, cutHeight = 0.995,                                 minClusterSize = minModuleSize,                                 pamRespectsDendro = F)unmergedColors <-  labels2colors(unmergedLabels)table(unmergedLabels)``

### 10.2 可视化

``sizeGrWindow(8,6)plotDendroAndColors(consTree, unmergedColors, "Dynamic Tree Cut",                    dendroLabels = F, hang = 0.03,                    addGuide = T, guideHang = 0.05)``

### 10.3 计算模块MEs

``unmergedMEs <-  multiSetMEs(multiExpr, colors = NULL, universalColors = unmergedColors)consMEDiss <-  consensusMEDissimilarity(unmergedMEs)consMETree <-  hclust(as.dist(consMEDiss), method = "average")``

### 10.4 可视化聚类结果

``sizeGrWindow(7,6)par(mfrow = c(1,1))plot(consMETree, main = "Consensus clustering of consensus module eigengenes",     xlab = "", sub = "")abline(h=0.25, col = "red")``

### 10.5 合并

``merge <-  mergeCloseModules(multiExpr, unmergedLabels, cutHeight = 0.25, verbose = 3)``

### 10.6 提取重要数据

``moduleLabels <-  merge\$colorsmoduleColors <-  labels2colors(moduleLabels)consMEs <-  merge\$newMEs``

### 10.7 最终可视化

``sizeGrWindow(9,6)plotDendroAndColors(consTree, cbind(unmergedColors, moduleColors),                    c("Unmerged", "Merged"),                    dendroLabels = F, hang = 0.03,                    addGuide = T, guideHang = 0.05)``

## 11Save一下

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

📍 往期精彩

##### jamesbang
V1

wx🔍: Grassssss 卷起来了