TCGA 数据挖掘实战 —— WGCNA 模块与临床表型关联分析

1. WGCNA与TCGA数据挖掘的核心价值

在癌症研究领域，TCGA数据库就像一座金矿，存储着来自上万例肿瘤样本的多组学数据。而WGCNA（加权基因共表达网络分析）则是挖掘这座金矿的强力工具，它能从海量基因表达数据中识别出具有生物学意义的共表达模块。我处理过上百个TCGA数据集，发现这种组合特别适合寻找肿瘤发生发展中的关键分子机制。

与传统差异表达分析不同，WGCNA关注的是基因之间的协同变化模式。举个例子，就像在社交网络中识别兴趣小组——不是看单个人的发言内容，而是观察哪些人经常互动、形成稳定社群。这种方法能捕捉到那些表达量变化不大但功能协同的重要基因，这些基因往往被常规分析方法遗漏。

实战中最关键的三个优势：

• 系统视角：将数万个基因简化为几十个功能模块
• 抗噪能力：通过拓扑重叠矩阵(TOM)过滤随机噪声
• 表型关联：直接链接分子特征与临床指标

2. 数据预处理的关键细节

2.1 表达矩阵标准化处理

拿到TCGA的RNA-seq数据后，我通常从FPKM或TPM值开始。这里有个容易踩的坑：不要直接用raw counts做WGCNA！曾经有个合作者坚持用counts数据，结果模块划分完全混乱。正确的做法是：

# 从TCGA获取TPM数据示例 library(TCGAbiolinks) query <- GDCquery( project = "TCGA-BRCA", data.category = "Transcriptome Profiling", data.type = "Gene Expression Quantification", workflow.type = "STAR - Counts" ) GDCdownload(query) data <- GDCprepare(query) tpm <- assay(data, "tpm_unstrand")

2.2 基因过滤的智能策略

过滤低表达基因时，我推荐使用**中位数绝对偏差(MAD)**而不是简单按平均值。因为某些重要调控基因可能在某些亚型中高表达，而在其他亚型中低表达。我的经验法则是：

# 取MAD前5000的基因 mad_values <- apply(tpm, 1, mad) filtered_genes <- names(sort(mad_values, decreasing=TRUE))[1:5000] exp_filtered <- tpm[filtered_genes,]

2.3 样本质量控制的实战技巧

样本聚类时经常会发现离群样本，这时需要谨慎处理。我开发过一个双重判断流程：

1. 先用hclust观察整体聚类结构
2. 再用PCA确认离群样本是否在主要成分上偏离

# 样本聚类诊断 sample_tree <- hclust(dist(t(exp_filtered)), method="average") plot(sample_tree, cex=0.6)

遇到批次效应时，我更喜欢用ComBat而不是直接删除样本，毕竟TCGA的珍贵样本删一个少一个。

3. 共表达网络构建的进阶技巧

3.1 软阈值选择的经验法则

选择软阈值功率(power)时，新手常犯的错误是盲目追求scale-free拓扑。实际上当样本量<50时，R²达到0.8可能很困难。我的调整策略是：

powers <- c(1:10, seq(12,20,2)) sft <- pickSoftThreshold(t(exp_filtered), powerVector=powers) # 当R²不理想时的备选方案 if(sft$powerEstimate < 5){ power <- ifelse(ncol(exp_filtered)<30, 8, 6) } else { power <- sft$powerEstimate }

3.2 模块识别中的参数调优

blockwiseModules函数有多个关键参数需要特别关注：

• mergeCutHeight：控制模块合并的阈值，我通常设为0.15-0.25
• minModuleSize：根据数据集大小调整，小样本设为30-50
• pamRespectsDendro：设为FALSE可以获得更自然的模块划分

net <- blockwiseModules( t(exp_filtered), power = power, TOMType = "unsigned", minModuleSize = 30, mergeCutHeight = 0.2, numericLabels = TRUE, verbose = 3 )

3.3 模块可视化技巧

用plotDendroAndColors展示结果时，我习惯添加两个改进：

1. 调整hang参数使树状图更紧凑
2. 使用自定义颜色板突出关键模块

library(WGCNA) moduleColors <- labels2colors(net$colors) plotDendroAndColors( net$dendrograms[[1]], moduleColors, "Module colors", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05 )

4. 模块-表型关联的深度分析

4.1 临床数据整合的实用方法

TCGA临床数据需要特别注意三点：

1. 样本ID匹配（barcode前12位）
2. 分类变量的标准化处理
3. 缺失值处理策略

# 获取ER/PR/HER2状态示例 clin_query <- GDCquery( project = "TCGA-BRCA", data.category = "Clinical", data.type = "Clinical Supplement", file.type = "patient" ) clin_data <- GDCprepare(clin_query)$clinical_patient_brca

4.2 关联分析的统计策略

除了默认的Pearson相关，我在不同场景下会选用：

• 二分变量：点二列相关(point-biserial)
• 生存数据：Cox回归
• 多分类变量：方差分析(ANOVA)

# 构建设计矩阵示例 design <- model.matrix(~0 + clin_data$er_status_by_ihc) modTraitCor <- bicor(MEs_col, design) modTraitP <- corPvalueStudent(modTraitCor, nSamples)

4.3 关键基因筛选的四象限法

我开发了一个可视化筛选方法：

1. x轴：基因与模块的相关性(MM)
2. y轴：基因与表型的相关性(GS)
3. 划分四个象限，优先选择右上象限基因

module <- "brown" genes <- names(net$colors)[moduleColors==module] MM <- abs(geneModuleMembership[genes, paste0("ME",module)]) GS <- abs(geneSignificance[genes, "Trait"]) plot(MM, GS, col=module, pch=16, xlab="Module Membership", ylab="Gene Significance") abline(h=median(GS), v=median(MM), lty=2)

5. 生物学解释与验证策略

5.1 模块功能注释的三步法

1. GO/KEGG富集分析：使用clusterProfiler
2. 蛋白互作网络验证：导入STRING数据库
3. 文献挖掘：用pubmed.mineR查找支持证据

library(clusterProfiler) ego <- enrichGO( gene = moduleGenes, OrgDb = org.Hs.eg.db, keyType = "SYMBOL", ont = "BP" ) dotplot(ego, showCategory=15)

5.2 生存分析的注意事项

做KM生存曲线时要注意：

• 用中位数表达量分组更稳健
• 多因素Cox回归控制临床混杂因素
• 考虑时间依赖性ROC分析

library(survival) fit <- survfit( Surv(OS.time, OS) ~ ifelse(MM>median(MM),"High","Low"), data = survival_data ) ggsurvplot(fit, risk.table=TRUE)

5.3 实验验证的设计建议

从计算到实验，我推荐这些验证策略：

1. qPCR验证：选择top 5-10个关键基因
2. 免疫组化：在独立队列验证蛋白表达
3. 细胞实验：敲除/过表达关键基因

在最后一次乳腺癌项目中，我们通过这种方法发现了TUBB6作为新的预后标志物，相关成果已发表在肿瘤学期刊上。整个过程从数据分析到实验验证大约需要6-8个月时间，但这样的闭环研究才能真正产生临床价值。