加权基因共表达WGCNA¶

一句话概述：WGCNA（加权基因共表达网络分析）将数千个基因按表达模式的相似性分成"模块"，再找出与表型相关的模块和hub基因，是转录组数据挖掘的高级武器。

核心知识点速览¶

概念	白话解释
WGCNA	把"一起涨一起跌"的基因归为同一个模块，再看模块和表型的关系
软阈值（β）	把相关系数"放大差异"的指数，让强相关更强、弱相关消失
无标度网络	少数节点连接很多（hub），大多数节点连接很少，像社交网络
模块（Module）	一群表达模式高度相似的基因，通常对应一个生物功能
模块特征基因（ME）	模块内所有基因表达的"代表值"（第一主成分）
模块-性状相关性	ME和临床表型的相关系数，找出哪个模块和你关心的疾病有关
模块成员（MM）	一个基因与其所属模块ME的相关性，MM高=该模块的核心
基因显著性（GS）	一个基因与临床性状的相关性
Hub基因	MM高且GS高的基因=既是模块核心又与疾病相关
TOM	拓扑重叠矩阵，考虑共享邻居的相似性度量

一、WGCNA原理（白话版）¶

想象你在管理一个100人的公司：

1. 计算相关性：看哪些员工"同进同退"（表达相关性）
2. 选软阈值：确定"多相关算熟人"的标准（β参数）
3. 建网络：画出所有人的"交友关系图"
4. 分部门：把关系紧密的人分到同一部门（模块）
5. 找部门负责人：每个部门选一个代表（ME）
6. 评估部门贡献：看哪个部门和公司业绩最相关
7. 找核心员工：在关键部门中找最核心的人（hub基因）

二、完整实操流程¶

2.1 安装和加载¶

# 安装WGCNA
install.packages("WGCNA")        # 从CRAN安装
# 或 BiocManager::install("WGCNA")  # 从Bioconductor安装

# 加载必要的包
library(WGCNA)        # WGCNA核心包
library(DESeq2)       # 数据标准化（可选）
library(ggplot2)      # 画图

# WGCNA推荐设置：允许多线程
allowWGCNAThreads(nThreads = 4)  # 设置线程数，加速计算

2.2 数据准备¶

# 读取表达矩阵（行=样本，列=基因）
# 注意：WGCNA要求行是样本，列是基因（和常规表达矩阵转置！）
expr_data <- read.csv("expression_matrix.csv", row.names = 1)  # 读取数据
expr_data <- t(expr_data)  # 转置：行变样本，列变基因

# 数据过滤：只保留变异较大的基因（减少计算量）
gene_var <- apply(expr_data, 2, var)           # 计算每个基因的方差
top_genes <- names(sort(gene_var, decreasing = TRUE))[1:5000]  # 取方差最大的5000个
expr_data <- expr_data[, top_genes]            # 过滤

# 检查数据质量
gsg <- goodSamplesGenes(expr_data, verbose = 3)  # 检查样本和基因质量
gsg$allOK  # 如果TRUE表示全部通过

# 如果有不合格的样本或基因，移除它们
if (!gsg$allOK) {
  expr_data <- expr_data[gsg$goodSamples, gsg$goodGenes]  # 只保留合格的
}

# 样本聚类检查离群值
sampleTree <- hclust(dist(expr_data), method = "average")  # 层次聚类
plot(sampleTree, main = "Sample Clustering Dendrogram")      # 画聚类树

# 如果有明显离群样本，手动设阈值移除
# abline(h = 15000, col = "red")  # 画切割线
# clust <- cutreeStatic(sampleTree, cutHeight = 15000)
# keepSamples <- (clust == 1)
# expr_data <- expr_data[keepSamples, ]

2.3 选择软阈值（β参数）¶

# 测试一系列候选软阈值
powers <- c(1:20)  # 测试1到20

# 计算每个候选值的网络拓扑
sft <- pickSoftThreshold(
  expr_data,                # 表达矩阵
  powerVector = powers,     # 候选软阈值列表
  networkType = "unsigned", # 网络类型：unsigned/signed/signed hybrid
  verbose = 5               # 输出详细信息
)

# 画软阈值选择图
par(mfrow = c(1, 2))  # 两张图并排

# 左图：Scale-free topology fit index vs power
plot(sft$fitIndices[, 1],                       # X轴=power
     -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],  # Y轴=R²×符号
     xlab = "Soft Threshold (power)",
     ylab = "Scale Free Topology Model Fit (R²)",
     main = "Scale Independence")
abline(h = 0.85, col = "red")  # R² > 0.85 的参考线
text(sft$fitIndices[, 1],      # 标注数字
     -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
     labels = powers, col = "red")

# 右图：Mean connectivity vs power
plot(sft$fitIndices[, 1],             # X轴=power
     sft$fitIndices[, 5],             # Y轴=平均连接度
     xlab = "Soft Threshold (power)",
     ylab = "Mean Connectivity",
     main = "Mean Connectivity")
text(sft$fitIndices[, 1], sft$fitIndices[, 5], labels = powers, col = "red")

# 选择软阈值的原则：
# 1. R² > 0.85 的最小power值
# 2. 如果找不到满足条件的，推荐用 sft$powerEstimate
# 3. unsigned网络通常 β=6-12，signed网络通常 β=12-20
soft_power <- sft$powerEstimate  # 自动推荐的软阈值
cat("推荐的软阈值:", soft_power, "\n")

2.4 构建网络和检测模块¶

# 一步法构建网络和检测模块（推荐新手使用）
net <- blockwiseModules(
  expr_data,                    # 表达矩阵
  power = soft_power,           # 选定的软阈值
  networkType = "unsigned",     # 网络类型
  TOMType = "unsigned",         # TOM类型
  minModuleSize = 30,           # 最小模块大小（基因数）
  reassignThreshold = 0,        # 重新分配阈值
  mergeCutHeight = 0.25,        # 合并相似模块的阈值（0.25=相关>0.75的合并）
  numericLabels = TRUE,         # 用数字标签
  pamRespectsDendro = FALSE,    # PAM是否遵循聚类树
  saveTOMs = TRUE,              # 保存TOM矩阵
  saveTOMFileBase = "TOM",      # TOM文件前缀
  verbose = 3,                  # 输出详细信息
  maxBlockSize = 6000           # 最大块大小（与内存有关）
)

# 查看结果
table(net$colors)  # 每个模块的基因数

# 画模块聚类树
mergedColors <- labels2colors(net$colors)  # 数字标签转颜色
plotDendroAndColors(
  net$dendrograms[[1]],        # 聚类树
  mergedColors[net$blockGenes[[1]]],  # 对应颜色
  "Module colors",              # 图例标题
  dendroLabels = FALSE,         # 不显示基因名（太多）
  hang = 0.03,                  # 悬挂比例
  addGuide = TRUE,              # 添加引导线
  guideHang = 0.05              # 引导线位置
)

2.5 模块与临床性状关联分析¶

# 读取临床性状数据
traits <- read.csv("clinical_traits.csv", row.names = 1)  # 临床数据
# 确保样本顺序和表达矩阵一致
traits <- traits[rownames(expr_data), ]

# 计算模块特征基因（ME）
MEs <- moduleEigengenes(expr_data, mergedColors)$eigengenes  # 每个模块的代表值
MEs <- orderMEs(MEs)  # 按相似性排序

# 计算模块-性状相关性
moduleTraitCor <- cor(MEs, traits, use = "p")           # Pearson相关系数
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(expr_data))  # p值

# 画热图
textMatrix <- paste(
  signif(moduleTraitCor, 2), "\n(",     # 相关系数
  signif(moduleTraitPvalue, 1), ")", sep = ""  # p值
)
dim(textMatrix) <- dim(moduleTraitCor)

par(mar = c(6, 8.5, 3, 3))  # 调整边距
labeledHeatmap(
  Matrix = moduleTraitCor,          # 相关系数矩阵
  xLabels = names(traits),          # X轴=性状名
  yLabels = names(MEs),             # Y轴=模块名
  ySymbols = names(MEs),
  colorLabels = FALSE,
  colors = blueWhiteRed(50),        # 蓝白红配色
  textMatrix = textMatrix,          # 显示数值
  setStdMargins = FALSE,
  cex.text = 0.5,                   # 文字大小
  zlim = c(-1, 1),                  # 颜色范围
  main = "Module-trait relationships"  # 标题
)

2.6 找Hub基因¶

# 选择感兴趣的模块（比如和疾病最相关的模块）
module_of_interest <- "blue"  # 假设blue模块最相关

# 计算模块成员（Module Membership, MM）
# MM = 基因与模块ME的相关性
MM <- cor(expr_data, MEs[, paste0("ME", module_of_interest)], use = "p")

# 计算基因显著性（Gene Significance, GS）
# GS = 基因与临床性状的相关性
trait_of_interest <- traits$disease_status  # 选择感兴趣的性状
GS <- cor(expr_data, trait_of_interest, use = "p")

# 画 MM vs GS 散点图
plot(abs(MM[mergedColors == module_of_interest]),  # X轴=|MM|
     abs(GS[mergedColors == module_of_interest]),  # Y轴=|GS|
     xlab = "Module Membership",
     ylab = "Gene Significance",
     main = paste("MM vs GS in", module_of_interest, "module"),
     pch = 19, col = module_of_interest)

# 筛选Hub基因：MM > 0.8 且 GS > 0.2
module_genes <- names(mergedColors)[mergedColors == module_of_interest]
hub_filter <- abs(MM[module_genes, ]) > 0.8 & abs(GS[module_genes, ]) > 0.2
hub_genes <- module_genes[hub_filter]
cat("Hub基因数:", length(hub_genes), "\n")
print(hub_genes)  # 打印hub基因列表

三、参数选择指南¶

参数	推荐值	说明
样本数	≥15，推荐≥25	样本太少WGCNA不可靠
基因数	3000-10000	取方差最大的基因
soft power	R²>0.85的最小值	unsigned一般6-12
networkType	unsigned/signed	signed更严格，推荐signed hybrid
minModuleSize	30	模块最少包含30个基因
mergeCutHeight	0.25	相关>0.75的模块合并

常见报错与解决¶

报错信息	原因	解决方案
`Scale free topology R² < 0.8`	数据可能不适合WGCNA	检查离群样本，增加过滤，或手动选β
`Memory allocation error`	基因太多内存不够	减少基因数或使用blockwiseModules
`goodSamplesGenes: removing genes`	有全零或缺失值的基因	让gsg自动过滤即可
`模块太少（只有grey）`	软阈值选错或数据质量差	重新选软阈值，检查数据标准化
`No module-trait correlation`	样本太少或批次效应	增加样本，校正批次效应
`Error in TOM calculation`	矩阵包含NA	检查expr_data中的NA值

速查表¶

# WGCNA完整流程
数据准备(转置+过滤) → 选软阈值(pickSoftThreshold)
→ 建网络(blockwiseModules) → 模块-性状关联(cor+heatmap)
→ 找hub基因(MM>0.8 & GS>0.2)

# 关键函数
goodSamplesGenes()      — 数据质量检查
pickSoftThreshold()     — 选软阈值
blockwiseModules()      — 一步法建网络+分模块
moduleEigengenes()      — 计算模块特征基因
cor() + corPvalueStudent() — 模块-性状相关
labeledHeatmap()        — 画模块-性状热图

# 网络类型选择
unsigned: 不区分正负相关，模块内基因可能正相关也可能负相关
signed:   只考虑正相关，更保守
signed hybrid: 折中方案，推荐

# Hub基因标准
模块成员 MM > 0.8（和模块高度一致）
基因显著性 GS > 0.2（和表型相关）

面试高频问题¶

Q1：WGCNA的"加权"是什么意思？¶

答：传统网络用阈值把相关性二分为"有连接/无连接"（如r>0.5就连上），这很粗糙。WGCNA的"加权"是把相关系数取β次方（如r^6），这样r=0.9变成0.53但r=0.5变成0.016，等于用连续权重代替了0/1的硬分界。β参数通过拟合无标度网络拓扑来确定，使网络符合生物学的"少数hub基因控制多数普通基因"的特征。

Q2：软阈值怎么选？R²的含义是什么？¶

答：R²衡量的是网络在该β值下有多符合无标度拓扑模型。选β的标准：①R²>0.85（达到无标度特征）的最小β值；②同时保持足够的连通性（mean connectivity不能太低）。如果所有β值R²都达不到0.85，可能说明数据有问题（批次效应、样本太少、或生物体系本身不满足无标度假设）。

Q3：WGCNA至少需要多少样本？¶

答：官方建议至少15个样本，推荐25+个。样本太少会导致相关性估计不准确，模块检测不稳定。对于临床数据，两组各10个（共20个）是最低要求。如果只有6-8个样本，WGCNA结果需要非常谨慎解读。

Q4：模块特征基因（ME）是什么？¶

答：ME是模块内所有基因表达矩阵的第一主成分（PC1），代表了该模块的"总体表达趋势"。可以理解为给一个班级选"班长"——这个人的成绩变化趋势最能代表全班平均水平。后续分析中，用ME代表整个模块和临床性状做相关分析，大大简化了计算。

Q5：如何判断Hub基因的可靠性？¶

答：①统计层面——Hub基因要同时满足高MM（>0.8）和高GS（>0.2）；②交叉验证——在独立数据集中验证模块的可重复性；③生物学验证——Hub基因是否有文献支持、是否在PPI网络中也处于核心位置；④功能验证——理想情况下应该有实验验证（如敲低/过表达实验）。发文章时通常结合PPI网络（STRING+CytoHubba）和WGCNA两种方法取交集。