
WGCNA分析:基因共表达网络、表型关联及聚类模块分析
最近在搞基因数据分析,发现WGCNA这玩意儿真是科研狗的救命稻草。今天就唠唠怎么用R语言玩转
这个工具,顺便把代码掰碎了喂给大伙儿。
先来个数据预处理热热身。拿到基因表达矩阵先别急着跑分析,这行代码能帮你快速排查坑爹的缺
失值:
```r
library(WGCNA)
options(stringsAsFactors = FALSE)
datExpr <- read.csv("你的表达矩阵.csv")
gsg = goodSamplesGenes(datExpr, verbose = 3)
if (!gsg$allOK) {
datExpr = datExpr[gsg$goodSamples, gsg$goodGenes]
}
```
这段代码的骚操作在于自动识别异常样本和基因。那个goodSamplesGenes函数就像个安检机,把表
达量全为零或者方差太小的基因直接踢出群聊。记得检查输出里的"Flagging genes and samples with
too many missing values..."这句话,它能告诉你数据被清洗得多干净。
接下来是重头戏——构建共表达网络。这里有个玄学参数叫软阈值,选不好整个分析直接翻车。试试
这个暴力测试法:
```r
powers = c(1:20)
sft = pickSoftThreshold(datExpr, powerVector = powers, networkType = "signed")
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)", ylab="Scale Free Topology Model Fit")
```
画出来的曲线要是像过山车似的,选那个R接近0.9的power值就稳了。一般哺乳动物数据用12左右
,植物可能得飙到18,具体看你的数据硬不硬核。
模块挖掘环节绝对是见证奇迹的时刻。跑完下面这段代码,你的基因就会自动抱团:
```r
net = blockwiseModules(datExpr, power = 12,