WGCNA分析（一）导入表达数据和寻找最佳β值

安装WGCNA

BioManager::install("impute")
BioManager::install("preprocessCore")
install.packages("WGCNA")

安装过程可能会不太顺利，大家可以自行查找解决方案。

导入表达数据

library(WGCNA)
exp0 <- read.delim(test1, row.names=1)

有很多基因是不表达或者表达了很低，所以我们过滤掉在所有样品中表达之和小于5的基因

exp <- exp0[rowSums(exp0) > 5, ]

我们的表达矩阵通常每一行是一个基因，每一列是一个样本；但WGCNA要求每一行是一个样本，每一列是一个基因，所以我们需要进行转置

Texp0 = t(exp0)

寻找最佳 β 值

# 开启多线程模式
  enableWGCNAThreads(nThreads = 20)
  
# 通过对 power 的多次迭代，确定最佳 power
  powers <- c(1:20)
x = pickSoftThreshold(Texp0, powerVector = powers, verbose = 5,
                          networkType = "signed")
# 计算出来的最佳 β 存放在
 x$powerEstimate
  
# 画图
## 设置画布
  sizeGrWindow(9, 5)
  par(mfrow = c(1,2))
  cex1 = 0.9

#  Scale independence
  plot(x$fitIndices[,1], -sign(x$fitIndices[,3])*x$fitIndices[,2],
       xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
       main = paste("Scale independence"));
 
# 增加一条阈值线
  abline(h=0.8,col="red")

# Mean connectivity 
  plot(x$fitIndices[,1], x$fitIndices[,5],
       xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
       main = paste("Mean connectivity"))