770. 基因共表达网络WGCNA进阶¶

一句话概述：WGCNA进阶分析包括签名网络(signed network)、模块保存性检验、多数据集比较——就像从"基础关系网"升级到"深度社交分析"，不仅知道谁跟谁一起，还能验证这些关系在不同环境下是否稳定。

核心知识点速查表¶

概念	白话解释	关键工具
Signed network	区分正相关和负相关	推荐使用
软阈值(β)	控制网络的scale-free特性	pickSoftThreshold
模块保存性	模块在其他数据集中是否稳定	modulePreservation
模块-表型关联	模块与临床性状的相关性	Module-trait关系
Hub基因	模块内连接度最高的基因	kME值
bicor	双权重中位相关(稳健)	替代Pearson

一、进阶概念¶

1.1 Signed vs Unsigned Network¶

Unsigned network（无符号网络）：
  相关性取绝对值 → |cor(A,B)|
  问题：正相关和负相关的基因被放在同一个模块
  → 生物学上不合理（激活和抑制的基因不应在一起）

Signed network（签名网络，推荐）：
  相关性保留符号 → (1 + cor(A,B)) / 2
  → 只有正相关的基因才可能在同一个模块
  → 负相关的基因在不同模块
  → 生物学解释更合理

Signed hybrid（签名混合，折中）：
  正相关 → cor(A,B)^β
  负相关 → 0
  → 忽略负相关但不混入正相关模块

1.2 相关性方法选择¶

Pearson相关（默认）：
  假设线性关系 + 正态分布
  对异常值敏感

Biweight midcorrelation (bicor，推荐)：
  稳健中位相关，对异常值不敏感
  WGCNA开发者推荐
  函数：WGCNA::bicor()

二、Signed Network构建¶

# ===== WGCNA Signed Network完整流程 =====
library(WGCNA)  # 加载WGCNA

# ===== Step 1: 数据准备 =====
options(stringsAsFactors = FALSE)  # 字符串不自动转因子
allowWGCNAThreads(8)               # 启用多线程

# 读取表达数据（基因为列，样本为行）
datExpr <- read.csv("expression_matrix.csv", row.names=1)
datExpr <- t(datExpr)  # 转置为样本×基因

# 检查缺失值和异常样本
gsg <- goodSamplesGenes(datExpr, verbose=3)  # 检查数据质量
if (!gsg$allOK) {
  datExpr <- datExpr[gsg$goodSamples, gsg$goodGenes]  # 移除问题数据
}

# 样本聚类检查异常值
sampleTree <- hclust(dist(datExpr), method="average")
plot(sampleTree, main="Sample Clustering to Detect Outliers")

# ===== Step 2: 选择软阈值（Signed Network）=====
powers <- c(1:20)  # 候选软阈值

# 关键：networkType = "signed"
sft <- pickSoftThreshold(
  datExpr,
  powerVector = powers,
  networkType = "signed",       # 签名网络（重要！）
  corFnc = "bicor",             # 使用bicor相关（推荐）
  verbose = 5
)

# 绘制拟合图
par(mfrow=c(1,2))
# 左图：Scale-free topology fit index
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)", ylab="Scale Free Topology Model Fit",
     main="Scale independence", type="n")
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers, col="red")
abline(h=0.85, col="red")  # R²>0.85的阈值线

# 右图：Mean connectivity
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)", ylab="Mean Connectivity",
     main="Mean connectivity", type="n")
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, col="red")

# 选择β值（R²>0.85的最小power，signed通常在10-14之间）
softPower <- sft$powerEstimate
cat("Selected soft threshold:", softPower, "\n")

# ===== Step 3: 一步构建网络和模块 =====
net <- blockwiseModules(
  datExpr,
  power = softPower,                # 软阈值
  networkType = "signed",           # 签名网络
  corType = "bicor",                # bicor相关
  TOMType = "signed",               # 签名TOM
  minModuleSize = 30,               # 最小模块大小
  reassignThreshold = 0,            # 重分配阈值
  mergeCutHeight = 0.25,            # 模块合并阈值
  numericLabels = TRUE,             # 数字标签
  pamRespectsDendro = FALSE,        # PAM设置
  saveTOMs = TRUE,                  # 保存TOM矩阵
  saveTOMFileBase = "signed_TOM",   # TOM文件名前缀
  verbose = 3,
  maxBlockSize = 20000              # 最大块大小
)

# 查看模块
table(net$colors)  # 每个模块的基因数
# 0    1    2    3    4    5    ...
# 500  1200 800  600  400  350  ...
# 模块0 = 灰色 = 未分配的基因

# 绘制模块树状图
mergedColors <- labels2colors(net$colors)  # 数字转颜色
plotDendroAndColors(
  net$dendrograms[[1]],
  mergedColors[net$blockGenes[[1]]],
  "Module colors",
  dendroLabels = FALSE,
  main = "Gene Dendrogram and Module Colors (Signed Network)"
)

三、模块保存性检验¶

# ===== Module Preservation Analysis =====
# 检验一个数据集中发现的模块在另一个数据集中是否保存

# 数据集1: 参考数据(已建好网络)
# 数据集2: 测试数据(验证模块是否存在)

# 准备多数据集格式
multiExpr <- list(
  Reference = list(data = datExpr_ref),   # 参考数据集
  Test = list(data = datExpr_test)        # 测试数据集
)

# 参考数据集的模块颜色
multiColor <- list(
  Reference = moduleColors_ref            # 参考数据的模块标签
)

# 运行模块保存性检验
mp <- modulePreservation(
  multiData = multiExpr,               # 多数据集
  multiColor = multiColor,             # 模块颜色
  referenceNetworks = 1,               # 参考网络索引(=1)
  nPermutations = 200,                 # 置换次数（越多越准，但越慢）
  randomSeed = 42,                     # 随机种子
  networkType = "signed",              # 网络类型
  verbose = 3
)

# 提取保存性统计量
stats <- mp$preservation$Z$ref.Test     # Z统计量
# Zsummary > 10: 强保存
# 2 < Zsummary < 10: 中等保存
# Zsummary < 2: 不保存

# 可视化
moduleSize <- stats[, "moduleSize"]
Zsummary <- stats[, "Zsummary.pres"]

plot(moduleSize, Zsummary,
     xlab="Module Size", ylab="Zsummary Preservation",
     main="Module Preservation", pch=19, col="blue")
abline(h=10, col="red", lty=2)   # 强保存线
abline(h=2, col="darkgreen", lty=2)  # 弱保存线
text(moduleSize, Zsummary, labels=rownames(stats), pos=3, cex=0.7)
legend("bottomright", c("Z>10: Strong", "Z>2: Moderate", "Z<2: Weak"),
       col=c("red","darkgreen","gray"), lty=2)

四、模块-表型关联¶

# ===== Module-Trait关系分析 =====
# 读取临床/表型数据
traits <- read.csv("clinical_traits.csv", row.names=1)
# 列：age, gender, tumor_stage, survival_status, treatment_response

# 计算模块特征基因(ME)
MEs <- moduleEigengenes(datExpr, mergedColors)$eigengenes
MEs <- orderMEs(MEs)  # 排序

# 计算模块-表型相关性
moduleTraitCor <- cor(MEs, traits, use="p")          # Pearson相关
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(datExpr))  # p值

# 热图可视化
textMatrix <- paste(
  signif(moduleTraitCor, 2), "\n(",
  signif(moduleTraitPvalue, 1), ")", sep=""
)
dim(textMatrix) <- dim(moduleTraitCor)

labeledHeatmap(
  Matrix = moduleTraitCor,
  xLabels = colnames(traits),
  yLabels = names(MEs),
  ySymbols = names(MEs),
  colorLabels = FALSE,
  colors = blueWhiteRed(50),
  textMatrix = textMatrix,
  setStdMargins = FALSE,
  cex.text = 0.5,
  main = "Module-Trait Relationships (Signed Network)"
)

# ===== Hub基因鉴定 =====
# kME = 基因与模块特征基因的相关性
# kME越高 → 基因对模块越重要（hub基因）

kME <- signedKME(datExpr, MEs, corFnc="bicor")  # 计算kME

# 找蓝色模块的hub基因
blue_genes <- colnames(datExpr)[mergedColors == "blue"]
blue_kME <- kME[blue_genes, "kMEblue"]
hub_genes <- names(sort(blue_kME, decreasing=TRUE))[1:20]  # Top 20 hub
cat("Blue module hub genes:\n", hub_genes, "\n")

五、常见报错与解决¶

报错信息	原因	解决方案
`pickSoftThreshold: R²<0.8`	数据不满足scale-free	增大power范围(至20)或检查数据
`blockwiseModules: memory`	基因太多	减小maxBlockSize或减少基因数
`Only grey module`	基因间相关性太低	降低minModuleSize或增加样本数
`Zsummary negative`	模块不保存	正常结果，说明模块无跨数据集意义
`bicor warning`	数据中有常数列	过滤方差为0的基因
`signed vs unsigned不同`	方法差异	signed更合理（推荐）

六、面试高频问题¶

Q1: 为什么推荐signed network？¶

A: Unsigned把正相关和负相关的基因放一起（只看绝对值），生物学上不合理。Signed只把正相关的基因聚成模块，负相关的在不同模块。例如：激活型基因和抑制型基因不应在同一模块。WGCNA开发者明确推荐signed network。

Q2: 模块保存性检验有什么用？¶

A: 验证在数据集A中发现的模块是否在数据集B中也存在。Zsummary>10说明模块稳健，有跨数据集的生物学意义；Zsummary<2说明模块可能是数据集特异的噪声。这在发现疾病biomarker模块时非常重要——只有保存的模块才可能有普遍的生物学意义。

Q3: 如何选择hub基因？¶

A: Hub基因选择标准：①kME值高(>0.8)：与模块特征基因强相关；②模块内连接度高：与模块内其他基因广泛相关；③基因显著性(GS)高：与感兴趣的表型强相关；④文献支持：已知与该生物学过程相关。最佳候选是同时满足高kME和高GS的基因。

七、速查表¶

# ===== WGCNA进阶速查 =====

# Signed网络构建
sft <- pickSoftThreshold(datExpr, networkType="signed", corFnc="bicor")
net <- blockwiseModules(datExpr, power=sft$powerEstimate,
    networkType="signed", corType="bicor", TOMType="signed",
    minModuleSize=30, mergeCutHeight=0.25)

# 模块保存性
mp <- modulePreservation(multiData, multiColor,
    referenceNetworks=1, nPermutations=200)
# Zsummary > 10: 强保存
# 2-10: 中等  | < 2: 不保存

# 模块-表型关联
moduleTraitCor <- cor(MEs, traits, use="p")

# Hub基因
kME <- signedKME(datExpr, MEs, corFnc="bicor")
hub <- names(sort(kME$kMEblue, decreasing=TRUE))[1:20]

# 关键参数选择：
# networkType = "signed" (推荐)
# corType = "bicor" (稳健)
# minModuleSize = 30 (经验值)
# mergeCutHeight = 0.25 (默认)
# power: signed通常10-14, unsigned通常6-12