TCGA甲基化数据挖掘这个系列写到第三期，前两期把下载和清洗讲完以后，收到的反馈里出现频率最高的一句话是：拿到甲基化矩阵之后，下一步到底该怎么验证一个分子表型？这一期我把完整的验证闭环拆开讲一遍，包括从差异甲基化位点筛选到分子表型评分、生存分析、独立队列复验，附上可以直接改路径就跑的R代码。适合已经开始用TCGA做表观遗传研究、但还不清楚怎么把甲基化信号和临床表型关联起来的同学。

这期代码的核心是解决一个常见尴尬：TCGA的甲基化数据下载下来是几十万行探针的β值矩阵，但想要回答“某个分子表型是否真的存在甲基化差异”时，很多人会被卡在中间环节——不知道该用什么模型算差异位点，也不知道算完之后怎么验证这个结果是稳健的，而不是某一组样本的偶然扰动。我把这个流程按我自己平时跑数据的习惯整理成四段：准备、筛选、验证、可视化，每段都有对应代码和参数说明。

1. 项目思路与整体流程设计

1.1 为什么要用甲基化数据验证分子表型

分子表型这个词听起来抽象，落到TCGA数据里其实就是你关心的样本分组，比如突变型对野生型、高表达组对低表达组、不同免疫亚型之间的比较。DNA甲基化作为表观遗传调控的重要机制，启动子区高甲基化通常伴随基因沉默，基因体区甲基化又和组织特异性表达相关，所以当一个分子表型成立时，大概率能在甲基化层面看到系统性的差异信号。

这就是验证的逻辑基础：如果A组和B组确实是分子层面不同的两类样本，它们之间的差异甲基化位点(DMP)和差异甲基化区域(DMR)应当显著多于随机分组。反过来，如果某个表型在甲基化层面完全无信号，那基本可以怀疑这个分组的生物学意义。用TCGA做这个验证还有个额外优势——大部分样本有完整的临床随访数据，可以直接把甲基化特征和预后做生存关联，这是很多自测样本不具备的。

1.2 常用数据源与选型对比

做TCGA甲基化挖掘，数据从哪拿、用哪种格式，直接决定后续代码怎么写。我整理了一个对比表，列出了三种主流获取方式，你可以根据自己的网络条件和后续分析需求选。

我自己常用的是UCSC Xena，因为它的甲基化矩阵是合并好的，样本注释也比较完整，拿下来就能跑，不用处理GDC那套复杂的文件目录结构。如果你要做离线归档或者后续要重新做探针滤波，GDC的原始层级数据是更好的起点，代价是需要处理更长的预处理链路。

1.3 分析流程总览

整套流程我拆成六个环节，每个环节都是下一步的输入依赖，执行顺序不建议乱跳：

1. 数据获取与探针过滤，输出干净的甲基化矩阵
2. 标准化与批次效应校正，消除非生物学差异
3. 差异甲基化分析，输出DMP/DMR表
4. 分子表型评分构建，把位点信号压缩成分数
5. 临床关联与生存分析，验证表型与预后的关系
6. 独立队列复验，排除偶然性和过度拟合

这六个环节里，第1、3、5步是卡住最多人的地方。后面我会把每步的关键参数和排查方式展开讲，重点解释“为什么这样选参数”而不是只丢一段代码。

2. 甲基化数据准备与预处理要点

2.1 三个下载渠道的取舍与注意事项

如果你决定用Xena的数据，下载ILMN 450K甲基化矩阵时有一个细节容易踩坑：Xena上提供的是基于不同注释版本的处理结果，有的列名带cg前缀，有的带ch前缀，后者的分析价值和研究范围都和前者不同。按我的经验，常规基因启动子甲基化分析选择cg探针即可，ch探针在Xena里通常指非CpG位点，分析工具支持有限，还容易和cg探针混在一起导致后续注释出错。

用GDC客户端的同学要注意版本对齐：TCGA的甲基化数据有多个数据版本（比如hg19和hg38对应的注释文件不同），下载时最好记录你用的是哪个版本的manifest，因为450K探针的坐标注释在不同基因组版本下会有偏移。我见过有人用hg19的探针注释去和hg38的表达矩阵做相关性分析，结果自然是错位百出，排查了很久才发现是参考基因组版本不一致。

Tcgabiolinks这个包的优势是它封装了GDC查询、下载、对象转换的完整流程，代码看起来会少很多。但它会创建一系列R对象，如果前面的探针过滤步骤没做，最后血糖矩阵可能会有百万行级别的探针，内存占用会很大，建议在分析前主动过滤。

2.2 β值、M值与探针过滤原理

450K芯片输出的原始信号经过分位数归一化后，通常得到β值，计算公式是：

β = Methylated / (Methylated + Unmethylated + offset)

β值反映的是甲基化水平百分比，0到1之间，解释直观，适合做可视化。但它有两个问题：第一，方差在接近0和接近1的两端被压缩，不符合常规线性模型的方差齐性假设；第二，不同探针的总体方差差异很大，直接做线性模型会放大噪声。

转换后的M值可以解决这个问题，它是β值经过logit变换得到的：

M = log2(β / (1 - β))

M值在全区间范围内近似正态分布，因此差异分析一般用M值跑模型，拿到差异位点后再换算回β值用于报告和可视化。建议在分析时同时保留β矩阵和M矩阵，两者有各自的用途，不要混用。

探针过滤的步骤以ChAMP包默认逻辑为基础，但我会加上自己的几条筛选条件。需要过滤的探针类型包括：检测p值缺失或大于阈值的探针、位于性染色体上的探针（避免性别差异干扰）、SNP位点附近的探针、以及多次比对到参考基因组的多义探针。过滤之后的探针数量通常会在42万左右，具体数量和原数据质量有关。

2.3 批次效应与标准化处理

TCGA的数据来自多个研究中心，不同批次之间的甲基化信号会存在系统性偏移。虽然450K在实验设计上已经做了很多技术标准化，但批次效应依旧存在，处理方式分为两种：一种是用ComBat做经验贝叶斯校正，另一种是用RCP或BMIQ做探针类型校正。

从我实际跑过的结果看，ComBat对跨平台合并数据效果更明显，但如果你的分析只用一个癌种的TCGA队列，批次的干扰通常不至于推翻差异结论，更多时候是让P值分布略微膨胀。这时候可以用一个更轻量的策略：在线性模型中加入样本来源批次作为协变量，而不是对全矩阵做全局校正。全局校正的风险是可能把真实的生物学差异抹掉一部分，尤其当批次和分组不完全正交时。

标准化这一步做完后，有一个值得做的质控：画一下样本间相关系数热图，确认同组样本聚在一起。如果发现明显的离群样本，先检查是不是探针过滤没生效，再考虑样本本身是否混入了重复或注释错误的样本。我遇到过一个个案，某样本在聚类热图中单独落在所有人之外，查了后发现它的芯片类型其实是上一代27K的数据，拿进450K矩阵里一起分析，自然怎么调参数都救不回来。

3. 差异甲基化分析核心代码与参数细节

3.1 基于limma的DMP分析完整代码

差异甲基化位点分析，我用的主力工具是limma包，底层逻辑和转录组差异表达类似。核心代码分为设计矩阵构建、线性拟合、经验贝叶斯检验三步：

# 假设M_mtx是探针x样本的M值矩阵，group_vec是分组因子 # 先做一个简单的两分组比较 library(limma) design <- model.matrix(~ 0 + group_vec) colnames(design) <- levels(factor(group_vec)) # 定义比较：A组减去B组 contrast.matrix <- makeContrasts(A - B, levels = design) fit <- lmFit(M_mtx, design) fit2 <- contrasts.fit(fit, contrast.matrix) fit3 <- eBayes(fit2) DMP_results <- topTable(fit3, coef = 1, number = Inf, sort.by = "none") DMP_results$beta_diff <- aveBeta_mtx[match(rownames(DMP_results), rownames(aveBeta_mtx)), "A"] - aveBeta_mtx[match(rownames(DMP_results), rownames(aveBeta_mtx)), "B"]

注意这里有个容易漏掉的点：拟合用的是M值矩阵，但筛选阈值时通常要看β值差异(Δβ)。因为M值的差异和β值差异不是线性关系，M值变化大不代表β值变化能通过临床阈值（常见Δβ绝对值>0.1或0.2）。所以我会额外加一列β值差异，筛选标准就定为“调整后P值小于0.05且Δβ绝对值大于0.1”。

3.2 DMR分析与注释

单个位点分析之后，如果信号比较零散，或者想找启动子区域的连续差异信号，建议做差异甲基化区域分析。我常用的是ChAMP包里的DMR模式，或者DMRcate包，两者选一个就行。

我做DMR分析时习惯这样设置参数：区域最小探针数3个，区域长度上限1000bp，每个探针的加权p值阈值取0.05。这类参数没有绝对标准，目的是把分散的位点压缩成有生物学意义的区域，阈值太严格会丢失真实信号，太宽松会堆积大量无意义区域。经验上先用默认参数看一遍区域数量，如果不超过200个，就可以直接继续往下走；如果数量爆炸，再逐步收紧探针数量阈值。

注释这一步主要是把位点/区域映射到基因、CpG岛位置、以及DNAse敏感区等特征上。可以借用annotatr或者ChAMP自带的注释函数。启动子区域的定义常见是转录起始位点上游2000bp到下游200bp之间，做注释时要把这个窗口写清楚，避免把远端的增强子区域和启动子混淆。

3.3 关键参数选择与计算逻辑

差异甲基化分析的参数选择，最核心的是两个：统计阈值和Δβ阈值。前者控制错误发现率，后者控制生物学效应大小。我遇到过不少初学同学只在代码里改了P值阈值就不管Δβ，结果筛出几千个位点，人工复核时发现绝大多数Δβ只有0.02——这种程度的甲基化差异在绝大多数生物学场景里都不能可靠区分样本，写文章也很容易被审稿人质疑。

另一个参数是探针筛选时要不要按方差排序。limma本身对低方差的探针会给出较小的权重，但如果你发现结果被一批恒定高甲基化的探针主导，可以考虑先对M矩阵做一次无监督的方差筛选，保留前25%高变异的探针再跑差异分析。这样做会让后续聚类和表型评分更稳健，缺点是可能丢失某些在某个亚组里才有差异信号的位点。

4. 分子表型验证的完整实现

4.1 构建甲基化评分与亚型聚类

拿到差异位点之后，下一个任务是把几十甚至几百个位点压缩成一个可解释的分子表型分数。两种常用方式：

第一种是直接基于DMP的方向打分：选定一组代表位点，计算每个样本在这些位点上甲基化水平的平均值，如果分组的表型是二分类，这个平均值就可以直接用于组间比较，给出检验结果。

第二种是聚类打分：用差异位点的β值做无监督聚类（层次聚类或一致性聚类），把样本分成K类，然后和你的目标表型做交叉表卡方检验。这种方式的优势是能发现预设分组外的自然亚型，缺点是需要选K，且可能得到与已知表型不一致的样本归类。

我自己的习惯是两种都做，先看聚类结果大概有几群，再根据先验知识决定是直接用目标表型还是聚类得到的表型。代码上注意一点：聚类之前把β矩阵做一次z-score标准化操作，不然甲基化水平整体偏高的样本会主导聚类距离。

4.2 生存分析与Cox回归验证

验证分子表型有没有临床价值，最直观的方法是做生存分析。TCGA的临床数据一般包含生存时间和生存状态，合并时要注意样本ID的匹配，尤其Xena下载的表达矩阵和临床表单的样本编码可能存在格式差异（有短横线和没有短横线）。

KM曲线用survival和survminer两个R包就能完成：

library(survival) library(survminer) # m_score是每个样本的甲基化评分，可取中位数分成高/低两组 surv_obj <- Surv(time = clin$OS_time, event = clin$OS_status) km_fit <- survfit(surv_obj ~ group) ggsurvplot(km_fit, pval = TRUE, risk.table = TRUE)

这里有一个很常见的坑：TCGA里有些样本的生存时间是NA，有些是0天。直接丢进Surv函数会报警告。处理方式是过滤生存时间缺失且状态缺失的样本，而不是把0直接当作缺失删掉。0天生存的样本通常代表术后早期死亡事件，不适合简单剔除，除非你有明确的入组排除标准。

Cox回归则用于多因素校正，把年龄、性别、分期作为协变量，验证甲基化评分是否独立于这些临床因素影响预后：

cox_fit <- coxph(surv_obj ~ m_score + age + stage + gender, data = clin) summary(cox_fit)

这里的核心检验指标是HR（风险比）及其95%置信区间，以及Wald检验的P值。如果单因素KM显著，多因素Cox回归不显著，说明甲基化评分可能只是临床分期的伴随信号，不算独立的预后因子，结论要谨慎下。

4.3 关联分析与独立队列复验

分子表型验证不能只看生存，还要看看它和已知生物学特征是否一致。比较通用的做法是：分析差异甲基化位点对应基因的表达水平，确认启动子甲基化与表达呈负相关。

# 假设beta值是启动子探针，expr是对应基因的mRNA表达 cor_res <- cor.test(beta_value, expr_value, method = "spearman")

对于探针与基因的对应关系，注意一个探针可能注释多个基因，一个基因也有多个探针，建议按“探针到基因最近转录起始位点距离最短”这个原则做唯一对应，避免重复计算。

独立队列复验是整套验证最有说服力的一环。理想情况是在另一个公开甲基化数据集（比如GEO里的某癌症队列）中重复同一种差异甲基化方向和表型关联。但独立队列常常没有和TCGA完全相同的肿瘤类型或甲基化芯片平台，这时候可以降低标准：不要求DMP完全一致，而是在GO/KEGG通路层面验证DMP富集的生物学通路与主队列一致，或者在独立队列中验证你构建的甲基化评分的预后方向是否一致（HR>1还是<1）。

5. 可视化与结果汇报

5.1 核心可视化图表出图思路

甲基化数据挖掘的报告展示一般需要四类图。第一类是火山图，展示DMP的P值和Δβ分布；第二类是热图，展示差异位点在样本间的表达模式；第三类是KM生存曲线，展示高低分组的预后差异；第四类是森林图，展示Cox回归中各因素的风险比。代码上火山图和森林图我常用ggplot2手工画，因为可以灵活调整配色和标签位置：

library(ggplot2) # 火山图：x轴是Δβ, y轴是-log10(adjusted P) ggplot(DMP_results, aes(x = beta_diff, y = -log10(adj.P.Val))) + geom_point(color = ifelse(abs(DMP_results$beta_diff) > 0.1 & DMP_results$adj.P.Val < 0.05, "red", "grey")) + theme_minimal() + labs(x = "Methylation difference (β)", y = "-log10(adjusted P)")

热图建议只展示筛选后的显著位点，不要一股脑把几万个探针画上去，否则聚类分支会被噪声淹没，也没有人能从图里获取有效信息。位点数量在20到100之间时热图的可读性最好。色板方面，甲基化热图常用蓝色到黄色或绿色到紫色，Z-score处理后用红蓝渐变更常见，我一般用pheatmap包。

5.2 结果解读与报告撰写的三条心得

结果解读这里容易犯的错误是只看P值不看效应量。统计显著和生物学显著不是一回事，一个位点调整后P值到小数点后十位但Δβ只有0.03，这对临床标志物的价值有限。报告时要同时标注P值和Δβ分布区间，让读者能判断这个位点的甲基化差异是否具备实际区分能力。

第二条心得是关注方向一致性。差异甲基化位点富集到某通路时，要单独检查该通路里位点的甲基化方向是整体上调还是整体下调。如果方向混乱，说明筛选可能混入了噪声或探针注释问题。通路方向的单一性比富集显著性更值得在报告里强调。

第三条心得是保存中间分析对象。R的RDS格式建议每个步骤都存一份，后面想换阈值重跑时不用从头下载数据。我自己会按01_raw、02_filtered、03_dmp、04_score这样的命名顺序存档，每一步的结果都能追溯到输入和参数，排查问题会高效得多。

6. 常见问题与排查技巧实录

6.1 高频报错与解决方案

我把这段时间被问得最多的问题整理成一张速查表，都是实际运行中会出现的情况和对应的处理思路：

6.2 数据一致性与复现性排查

甲基化数据挖掘的复现性很大程度取决于版本锁定。R版本、包版本、注释文件版本、参考基因组版本，任何一个变动都可能导致DMP表微调。如果要把分析流程写成可复现脚本，建议在开头的代码里用sessionInfo()记录环境，再用renv固定包的版本。

除此之外，还有一个容易被忽略的细节：450K探针在不同版本注释中的基因映射可能有变化。比如某探针在新版注释中映射到基因A，旧版却映射到基因B的转录起始位点上游。这种差异在通路富集分析时会被放大。稳妥的方案是以当前主流注释版本为准，不要混用多个版本的注释。

6.3 我在实际操作中总结的三个避坑技巧

第一，无论用什么流程拿到TCGA数据，第一步永远是检查样本ID和分组标签的对应关系。TCGA的barcode含义丰富（前四位是癌种、第七到第八位是样本类型），建议先把01A/01B这种肿瘤样本和11A这种正常样本挑出来，不要混在一起分析，不然正常样本混入会带来强烈但不真实的差异信号。

第二，差异甲基化位点的筛选阈值不要一开始就用很严格的FDR。先用宽松的阈值（比如P<0.05但是不做多重校正）跑一遍看看整体分布，确认没有技术性偏大效应后，再用FDR校正收窄结果。这样可以快速发现A/B组间是否有一整块芯片信号异常导致的假差异。

第三，把检索结果和已发表文献做一次交叉比对。TCGA的每个癌种都有大量已发表的甲基化标志物研究，跑完分析后去PubMed搜一下同类肿瘤中的热门基因，如果完全找不到已知基因，先怀疑数据或方法，而不是急着下结论说发现了新标志物。这个交叉验证习惯能帮你省下大量无效的后续实验验证时间。

这一系列流程走下来，整体收益很明显：它不仅验证了分子表型在甲基化层面的真实性，还把筛选出来的甲基化信号转化成了一个可以用于预后分层的量化指标。对我来说，最值得的是把独立队列复验纳入标准动作，这会逼你在分析的最开始就定义好组间差异的生物学预期，而不是到后面堆了一堆表格再试图解释。如果后续你在实际跑这些代码时遇到报错或者结果异常，建议先做减法——回到探针过滤和样本匹配这两步，大部分定位不准确的问题都藏在那里。