MRD-KM生存分析实操:用R语言绘制肿瘤预后生存曲线
2026/9/15 22:30:59 网站建设 项目流程

MRD-KM这个组合,最近几年在血液肿瘤和实体瘤的预后分析里出现频率越来越高。它本质上不是什么高深算法,而是把“微小残留病灶(MRD)检测”和“Kaplan-Meier生存曲线”这两样东西绑在一起,用来回答一个临床问题:MRD状态不同,患者的长期结局到底差多少?这篇内容我会从头到尾拆一遍这个分析怎么做,数据怎么准备、分组怎么定、R代码怎么写、有哪些坑,尽量让刚接触的人也能照着做出来一份能放进论文里的结果。

1. 项目整体设计与思路拆解

1.1 MRD是什么,为什么要看MRD

MRD,minimal residual disease,中文常翻译成微小残留病灶或可测量残留病灶。简单说,就是治疗后用高灵敏度的方法去检测体内还剩下多少肿瘤细胞或者肿瘤来源的核酸片段。传统影像学能看到的是毫米级以上的病灶,而MRD能把检测灵敏度推进到10^-4甚至10^-6的级别,也就是说每10万个正常细胞里如果混了1个肿瘤细胞,也能被揪出来。

那为什么临床上这么看重MRD?因为MRD状态直接反映治疗后的“深度缓解”程度。一个患者化疗后影像学看完全缓解,但MRD是阳性,说明体内还有潜伏的肿瘤细胞,这些人往往复发风险更高、生存期更短。相反,MRD持续阴性的患者,通常是真正意义上的“深度缓解”,预后要好得多。所以MRD已经不只是一个检测指标,在很多血液肿瘤诊疗指南里已经成了疗效评估和后续治疗决策的重要参考。

在科研场景里,我们经常需要把MRD跟生存数据结合起来,用数据说明“MRD阴性组比MRD阳性组活得久”。这时候KM分析就是最常用、也最容易讲清楚的统计工具。

1.2 KM生存分析解决什么问题

Kaplan-Meier生存分析,医疗圈几乎没有人不知道。它做的事情其实很朴素:在不同时间点上,计算还有多少人“活着”或者说“没有发生目标事件”,然后画出一条阶梯状的曲线。这里的“事件”可以是死亡、复发、疾病进展,都可以。KM曲线最直观的价值是让我们能目测两组或多组人群的生存差异,同时通过log-rank检验给出一个P值,判断差异是否具有统计学意义。

KM分析独有的优势在于它能处理删失数据。什么叫删失?比如一个患者入组后随访了18个月,到截止日期还活着,没有发生目标事件,那他的数据就不能简单当作“活了18个月”或“没活到终点”来处理。KM法会把他标记为删失,在曲线里用竖线/加号标识,并把他“未发生事件”的信息保留在风险集中,直到他离开随访队列。这个机制保证了生存估计不会因为部分人还没结束观察就被错误低估。

在做MRD相关生存分析时,KM图基本是标配。因为MRD本身是个分类变量(阴性/阳性),恰好天然适合用KM曲线去画两组甚至三组的生存差异。

1.3 为什么把MRD和KM放在一起看

把MRD和KM结合,本质上是在做一个“暴露因素与临床结局”的关联分析。MRD是暴露(或者说预测因素),临床结局是生存指标。通过KM曲线,我们能看到MRD阴性组的生存曲线明显“高位走”,MRD阳性组的曲线“低位掉”,两条曲线早早分开,log-rank P值越小,说明MRD状态和预后的关联越强。

但这里要强调一句:这种分析得出的是“相关”,不是“因果”。MRD阳性并不一定直接导致死亡或复发,它可能只是反映了肿瘤生物学行为更凶险、或者对治疗不敏感。所以论文里写结果时,要注意措辞,一般说“MRD阳性与更差的无复发生存显著相关”,而不要直接写“MRD阳性导致复发”。

从研究设计的角度看,MRD-KM相关分析常见于这几类场景:

  • 新诊断患者诱导治疗后评估MRD状态,预后分组;
  • 移植前后MRD动态变化与移植预后的关联;
  • 某种新药治疗过程中MRD转阴率与长期生存的关系;
  • 不同MRD检测方法(流式、PCR、NGS)分层后的一致性分析。

无论哪种场景,分析流程大方向一致:确定时间零点、定义MRD分组、整理随访数据、画KM曲线、做log-rank检验、如果有必要再上Cox回归校正混杂因素。

1.4 核心输出物长什么样

一份完整的MRD-KM分析,核心输出物通常包括三个部分:

  • KM生存曲线图(带风险表、标注组别、事件数、P值);
  • Log-rank检验结果(chi-square值、自由度、P值);
  • 如果做了多因素分析,还需要森林图或回归表,展示校正后的HR和置信区间。

很多期刊对KM图有具体要求,比如要标注删失、要显示number at risk、要标明横纵轴含义和单位。这些细节后面实操部分我会一一拆解。

2. 数据准备与分组策略

2.1 MRD检测方法与时间点选择

MRD检测方法会直接影响分组标准,这个必须先在数据层面定义清楚,否则后面统计再漂亮也是空中楼阁。

流式细胞术是血液肿瘤最常用的MRD检测手段,灵敏度通常在10^-4,也就是0.01%的水平。PCR方法比如针对融合基因的定量PCR,灵敏度能到10^-5甚至10^-6,但前提是你得有明确的可追踪分子标志。NGS测序方法近年用得越来越多,可以追踪多个克隆,灵敏度和PCR相当,但成本高、周期长、标准还没完全统一。

方案选型上有一个原则:全组患者最好用同一种检测方法。如果一部分人用流式、一部分人用PCR,那MRD阳性/阴性的判定标准本身就有差异,直接合并分析会引入很大的方法学偏倚。

时间点的问题也容易踩坑。MRD检测的时间点通常会取某个治疗节点,比如诱导治疗后、巩固治疗后、移植前、移植后某个时间点。在KM分析里,这个时间点必须跟“生存时间起点”对齐。最常见的设计是:以某个治疗节点作为时间零点,比如移植日,之后采集MRD状态,然后开始随访生存事件。如果患者在MRD检测之后又过了很久才移植,那时间零点的选择就要格外小心,要不要把检测时间点和移植时间点拉开,需要结合临床问题判断。

2.2 分组逻辑:阴性/阳性、持续阴性/转阳、清零时间

MRD分组方式不是只有“阴性vs阳性”这一种,具体怎么做取决于你想回答什么临床问题。

最基础的分组是单时间点MRD状态:比如取诱导治疗结束后某个时间点的MRD结果,分阴性和阳性两组,然后做生存曲线对比。这个做法简单直接,也是最多论文里出现的用法。问题在于它忽略了MRD的动态变化。

稍微进阶一点的分组是按MRD动态轨迹分类。比如:

  • 持续阴性:多个时间点MRD均为阴性;
  • 持续阳性:多个时间点MRD始终未转阴;
  • 先阴后阳:有过MRD阴性,随访中又转阳(分子复发)。

还有一种分组方式是按“MRD清零时间”分:治疗后多长时间内MRD达到阴性,早清零和晚清零的患者预后可能是不一样的。

具体选哪种,没有标准答案,核心逻辑是分组方式要能代表你要验证的假设。如果关注治疗后早期深度缓解的预后价值,单时间点就够了;如果要说明MRD动态监测的意义,那持续阴性/转阳这种轨迹分类更有说服力。

2.3 需要收集的临床数据字段

做MRD-KM分析,数据收集阶段就要规划好字段,后面统计时才不会手忙脚乱。我建议至少准备以下几类变量:

  • ID:患者唯一标识;
  • 分组变量:MRD阴性/阳性,或者MRD轨迹分类;
  • 生存时间变量:从时间零点到事件发生的时间(单位通常用月);
  • 结局事件变量:0=删失/未发生事件,1=发生了目标事件(死亡、复发等);
  • 随访截止日期和事件发生日期(用来计算生存时长);
  • 协变量:年龄、性别、疾病分型、危险分层、治疗方案、移植方式等,后续可能用于Cox回归校正。

另外,要特别确认结局事件的定义。比如“无复发生存”的事件定义是复发或任何原因死亡,以先发生者为准;“总生存”的事件定义只有死亡。不同的终点定义,对应的事件变量编码逻辑不一样,千万别搞混。

2.4 样本量预估与统计效能

样本量问题在MRD-KM分析中很容易被忽视。很多人拿到数据就开始画曲线,等画出来两组差异不显著,才开始后悔当初没算样本量。

估算样本量的输入参数包括:预期的两组生存率差异(比如MRD阴性组2年RFS预计80%,阳性组预计50%)、显著性水平(通常取0.05)、检验效能(通常取0.8)、随访时间和入组速度。软件方面可以用PASS、nQuery,R里也有powerSurvEpi、survival包能算。

如果样本量已经固定且不大,那就换一种思路:不要硬做两组比较,考虑用连续变量比如MRD水平的对数转换值做Cox回归,或者用限制性立方样条看剂量效应关系。有时候连续变量分析反而更容易发现关联,因为分组会损失信息。我见过不少样本量只有三四十例的研究,强行分两组阳性组只有十来例,log-rank P值基本不可能显著,这时候更务实的做法是定义为探索性分析,报告效应量和置信区间就好。

3. 实操过程:用R完成MRD-KM分析

3.1 数据清洗与生存对象构建

R语言做生存分析最核心的两个包是survival和survminer。survival提供计算框架,survminer负责画图和美观化。环境没配好的话先去装包。

install.packages("survival") install.packages("survminer")

假设你的数据已经整理成Excel或CSV,包含以下列:ID、MRD_group(Neg/Pos)、time_months(生存时间,单位月)、event(0/1)。读入数据的代码如下:

library(survival) library(survminer) dat <- read.csv("mrd_data.csv", stringsAsFactors = FALSE) str(dat) head(dat)

数据清洗阶段有几个关键动作:

  • 检查缺失值:MRD分组缺失、生存时间缺失、事件缺失,这些记录是否剔除要有明确规则。MRD缺失的记录通常直接剔除,因为分组本身没有了;生存时间缺失如果是因为失访,可以保留为删失,但必须记录清楚。
  • 检查生存时间是否合法:所有时间都应该是正数。如果出现0,要核对患者是不是在时间零点当天就发生了事件,如果是,要不要保留,取决于临床定义。一般来说,随访极短且事件发生在零点附近会影响曲线开头,需谨慎处理。
  • 检查事件变量编码:event必须是0或1。数据类型要转成数值型,否则后续Surv函数会报错。
dat$MRD_group <- factor(dat$MRD_group, levels = c("Neg", "Pos")) dat$event <- as.numeric(dat$event) dat$time_months <- as.numeric(dat$time_months) dat <- dat[complete.cases(dat$time_months, dat$event), ]

构建生存对象是KM分析的起点。R里用Surv()函数:

surv_obj <- Surv(time = dat$time_months, event = dat$event)

这里event=1表示事件发生,event=0表示删失。

3.2 KM曲线绘制与log-rank检验

用survival包拟合整体KM模型,再用survdiff做组间比较:

km_fit <- survfit(surv_obj ~ MRD_group, data = dat) summary(km_fit)

summary输出里可以看到各组每个时间点的生存概率、置信区间、删失信息和风险人数。

log-rank检验:

logrank_test <- survdiff(surv_obj ~ MRD_group, data = dat) print(logrank_test)

输出里有chisq、df和p-value。p值小于0.05说明两组生存曲线差异有统计学意义。

画图推荐survminer的ggsurvplot,因为它能自动生成专业级KM图,并且自带风险表:

pdf("km_plot.pdf", width = 8, height = 6) p <- ggsurvplot( km_fit, data = dat, pval = TRUE, pval.method = TRUE, conf.int = TRUE, risk.table = TRUE, risk.table.col = "strata", xlab = "Time (months)", ylab = "Relapse-free survival probability", legend.title = "MRD status", legend.labs = c("MRD negative", "MRD positive"), palette = c("#2E86AB", "#A23B72"), surv.median.line = "hv" ) print(p) dev.off()

pval = TRUE会自动添加log-rank检验P值,pval.method = TRUE会显示“Log-rank”字样。conf.int画出置信区间带。risk.table = TRUE会在曲线下方显示各时间点风险人数,这个几乎是SCI期刊标配。surv.median.line = "hv"会在中位生存时间位置画辅助线。

如果你想把曲线横纵轴刻度调整得更精细,可以用break.time.by参数设置横轴断点间距,比如break.time.by = 6代表每6个月一个刻度。y轴用ylim = c(0, 1)固定范围。

3.3 风险表(number at risk)与Cox回归辅助

风险表是KM图里很容易被忽略但审稿人特别爱看的东西。它告诉读者在每个时间点,队列里还有多少人在观察范围内没有发生事件、也没有被删失。没有风险表的KM曲线就像没有坐标的图,别人无法判断后半段曲线是基于多少人画的。

如果数据集够大,我建议也可以输出中位生存时间及其置信区间:

print(km_fit)

这个输出会显示每组的中位生存时间。中位生存时间的含义是生存概率降到50%的时间点,和风险表配合使用很能说明问题。

如果只是做“MRD阳性组中位无复发生存时间显著短于阴性组”,这还不够。生存分析更高级的阶段是用Cox比例风险回归,在考虑MRD分组的同时校正其他混杂因素:

cox_fit <- coxph(surv_obj ~ MRD_group + age + risk_group, data = dat) summary(cox_fit)

输出里最核心的是exp(coef),也就是风险比HR。比如MRD阳性组的HR=3.2,意味着在调整了年龄和危险分层后,MRD阳性患者的复发/死亡风险是阴性组的3.2倍。同时需要关注exp(coef)的95%置信区间和p值。

需要注意的是Cox回归的前提之一是比例风险假设,即各协变量对风险的影响在不同时间点是成比例的。可以用cox.zph()检验:

zph_test <- cox.zph(cox_fit) print(zph_test)

如果p值小于0.05,说明比例风险假设可能不满足,这时候要么考虑分层分析,要么用时间依协变量模型。

3.4 森林图快速绘制

做完多因素Cox回归,很多场合需要输出一张森林图来展示各组变量的HR和置信区间。R里forestplot包或者survminer自带的ggforest都可以,我一般用ggforest:

library(survminer) ggforest(cox_fit, data = dat)

一行代码就能生成比较规范的森林图,变量名、HR、95%CI、P值都在上面。如果想自定义样式,可以改用forestplot包,但基础需求ggforest足够。

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

4.1 删失数据处理不当

删失数据在KM分析里是核心机制,但也是新手最容易搞错的地方。最常见的错误是把失访当成“存活到研究结束”来处理,也就是把event标记成0(未发生事件),但时间却填了失访之前的最后一次随访时间,这在做法上其实是对的。如果是把失访患者当作“发生了死亡”来处理,比如错误地把时间填到死亡日期,那就是大问题,这等于人为创造了虚假事件,会严重低估生存概率。

另一个常见问题是删失标记不统一。有人用0表示存活、1表示死亡,有人用1表示事件、2表示删失,这些编码本身没问题,但关键在于读入R之后要正确转换。建议在读入数据后先table()检查一下event变量的取值分布,如果发现有些值是NA、9或者中文标签,那就要清洗后再建模。

4.2 零点(time zero)设错导致结果失真

这是MRD-KM分析里最隐蔽的坑。KM分析要求所有人从同一个起点(time zero)开始进入队列,然后观察后续事件发生情况。如果一部分人从入组开始计时,另一部分人从治疗结束开始计时,这就等于把不同起跑线的人放在一起比,结果必然失真。

临床研究中,时间零点的选择一定要和分组变量的评估时间对齐。比如你在移植后第100天测定MRD,那时间零点就应该是移植日,生存时间=事件发生日期-移植日期。如果你把生存时间算成了从确诊开始,那MRD是在确诊后很久才测的,两组在测MRD之前就已经有了一部分“随访期”,这些前置时间会让曲线开头失真。

解决思路是:把所有患者的时间轴头尾对齐到同一个治疗节点。如果某些患者的MRD检测时间点差异很大(比如有人3个月测、有人6个月测),严谨的做法是做landmark分析。也就是选定一个统一的时间点(比如治疗后6个月),只纳入随访到那个时间点且没有发生事件的患者,以那个时间点作为time zero,再按MRD状态分组。这个方法能有效消除时间窗口不齐的问题,在很多前瞻性研究里也是默认的分析策略。

4.3 样本量小、组间基线不齐怎么办

MRD相关研究经常面临样本量小的问题,尤其是单中心研究,某亚组可能只有十几例甚至几例。这时候log-rank检验的检验效能很低,即使实际预后差异很大,P值也可能不显著。

我的建议是“看效应量而不是只盯P值”。报告KM曲线时,标明不同时间点的生存率差值,比如“MRD阴性组2年RFS为78.5%,MRD阳性组为42.1%”,同时给出HR和95%CI,这比单一个P值更有信息量。

如果两组基线特征差异比较大,直接比较KM曲线可能会被年龄、分期等混杂因素干扰。这时候需要做多因素Cox回归进行调整,不能只靠KM图说话。如果样本量太小没法做多因素调整,那就明确写清楚这是探索性结果,结论慎重。

4.4 结果解读的坑:MRD不是静态标签

MRD-KM相关分析最大的一个容易被人诟病的问题是:MRD是检测时间点的状态,但在随访过程中可能发生变化。一个人可能在3个月时MRD阴性,到了9个月转阳。如果只拿3个月的MRD状态去做KM分析,那么“阴性组”里其实混了一部分后来转阳的人,这会稀释阴性组的生存优势。

这在方法学上属于“时间依赖偏倚”或者“不朽时间偏倚”。解决办法之一是前面提到的landmark分析,另一个更复杂的办法是做时间依赖ROC和time-dependent Cox模型。如果你的数据里有多个时间点的MRD结果,强烈建议做一次动态轨迹分析,看看跟单时间点分析相比,结论有没有变化。这个步骤在论文里写出来,审稿人会很认可,因为它表明你考虑到了MRD的动态属性。

4.5 软件与作图细节问题

用R做KM分析最常遇到的报错是Surv函数里time和event长度不一致,通常是数据里有NA被na.action处理掉之后,模型拟合时数据行数对不上。解决办法是在建模前主动用na.omit(dat)清理数据,或者在survfit和survdiff里设置na.action = na.exclude。

画图时的中文字体问题也很烦人。如果数据分组名或图例里用了中文,PDF输出很可能出现方块字。建议作图前先把分组名改成英文,或者用showtext包引入中文字体。我一般习惯直接全英文出图,后续在PPT或论文里再做中文标注。

ggsurvplot里有时候风险表被截断,多半是绘图区域高度不够或者risk.table.height参数没设好。可以适当调整height和width,或者用survminer包的arrange_ggsurvplots把图和风险表分开排布。

如果你最终需要把曲线放进Word里,导出PNG比PDF更方便。PNG分辨率至少300dpi,宽度一般要求10-15cm。R里可以用ggsave()指定输出参数,例如:

ggsave("km_plot.png", plot = p, width = 8, height = 6, dpi = 300)

有时候为了组合多个KM图(比如不同终点分别作图),可以用survminer包的arrange_ggsurvplots函数,或者用gridExtra拼图。这个在要压缩图片数量时特别好用。

4.6 辅助的敏感性分析

好的SCI论文,光有主分析的KM曲线还不够,通常还需要做敏感性分析来验证结果的稳健性。MRD-KM分析里常见的敏感性分析方法包括:

  • 换用不同的终点定义,比如从“无复发生存”换成“无事件生存”,看看结论是否一致;
  • 把MRD分组阈值在不同水平下重跑,比如用0.01%和0.001%两个阈值分别分组,观察HR变化;
  • 排除早期复发或早期死亡的患者重新做一次KM分析;
  • 用倾向性评分匹配法平衡基线特征后再比较MRD两组的生存差异。

这些分析不需要全部做,但至少选一两个,能大幅提升文章说服力。我做这个MRD-km项目的时候就发现,如果只看单时间点分析,MRD两组差异虽然显著,但Cox模型里加入高风险的细胞遗传学分层后,MRD的独立预后价值会被削弱。后面补做了一个landmark分析和一个多时间点轨迹分析,结果更稳定,故事也讲得更圆满。

我个人在实际操作中的体会是,MRD-KM相关分析真正的难点不在代码,而在数据准备和研究设计。画一条KM曲线五分钟就能搞定,但把MRD检测时间点、随访时间和分组逻辑理清楚,往往要花上好几天。尤其是随访数据不完整的中心,还得跟临床医生反复核对患者记录,那才是最耗心力的环节。做科研统计,前面的数据功夫到位了,后面的分析就是水到渠成的事情。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询