1写在前面

最近遇到了非常不讲理的病人和家属，真是忍不住想要吐槽两句。😤

为了让病人做上手术，求了输血科半天才给备了血，家属也答应去献血。😒

万万没想到，术后都2天了还是没去献血，无论是问病人还是问家属，都说家里没有一个人适合献血的。🤒

用别人的血就天经地义、理所当然，让你去献血就是大夫在害你了，我也是没得话讲，农夫与蛇的现代翻版。😷

言归正传了，继续之前的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

接着我们为每个数据集计算一下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.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]);
} else
text(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) ;
} else
legend("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$multiMEs
moduleLabels <- net$colors

moduleColors <-  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一下

保存一下data吧，毕竟运行了那么久。🥳

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

点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰