ReaderMX.
pythonic生物人

复现顶刊Kaplan-Meier curves ！

ppythonic生物人约 6 分钟

本次使用s41591-023-02790-x中的代码、数据来复现Nature medcine中的Kaplan-Meier curves。

Kaplan-Meier曲线常用于生存分析（Survival analysis）评估疗效。例如，下图中，在一个肝内胆管癌免疫治疗队列研究中，使用了某种免疫治疗药物后，TNB-high组患者的OS明显长于TNB-other组患者，可初步得出TNB-high组患者更受益于该药物。

doi: 10.1136/jitc-2019-000367
doi: 10.1136/jitc-2019-000367

下面使用R制作Kaplan-Meier curves，

✅数据预处理

# 设置图形输出参数：宽度8英寸，高度6英寸，分辨率200 DPI
options(repr.plot.width = 8, repr.plot.height = 6, repr.plot.res = 200)

# 加载生存分析和可视化所需的R包
library(survival)   # 用于生存分析，如Cox模型和生存曲线拟合
library(survminer)  # 用于绘制专业的生存曲线（Kaplan-Meier曲线）

#----------------------------------------------------------------------------------
#  Step 1: 加载数据
#----------------------------------------------------------------------------------
# (1) 读取191例PDAC患者（RJ队列1）的临床信息表
rj1.cohort <- read.xlsx("r_data/PDAC2_data_results/data/Extended Data Table 2.xlsx", startRow = 2)

# (2) 读取32个功能模块在癌组织（PDAC）和癌旁组织（TAC）中的ssGSEA评分矩阵，并进行转置（行为模块，列为样本）
module.ssgsea.pro.all <- read.xlsx("r_data/PDAC2_data_results/data/Extended Data Table 4.xlsx", sheet = 2, startRow = 2, rowNames = T)
module.ssgsea.pro.all <- t(module.ssgsea.pro.all) # 转置矩阵，使行名为模块名（如ME11），列名为样本ID

#----------------------------------------------------------------------------------
#  Step 2: 基于ME11模块进行总生存期（OS）分析
#----------------------------------------------------------------------------------
# 准备绘图数据：将ME11模块的评分添加到临床数据中，并根据其中位数将患者分为"High"/"Low"两组
plot.cohort <- rj1.cohort %>% 
  mutate(
    ME11_pro = as.numeric(module.ssgsea.pro.all["ME11", Proteomic_ID]), # 提取每个样本的ME11评分
    ME11_pro_l = ifelse(ME11_pro >= median(ME11_pro), "High", "Low"),    # 按中位数分组
    ME11_pro_l = factor(ME11_pro_l, levels = c("Low", "High")),         # 将分组设为因子，确定比较顺序
# 以下两行处理化疗状态因子，创建两种水平的因子用于可能的不同分析对比
    Censoring_chemo_rev = Censoring_chemo,
    Censoring_chemo = factor(Censoring_chemo, levels = c("Chem", "NoChem")),
    Censoring_chemo_rev = factor(Censoring_chemo_rev, levels = c("NoChem", "Chem"))
  )

head(plot.cohort)
Image

✅ 计算p值、CI等统计指标

  • Hazard Ratio：HR，风险比；
  • Confidence Interval：CI，置信区间；
  • Cox P-value；
  • Kaplan-Meier P-value；
# 执行Cox比例风险回归分析，检验ME11分组（High vs Low）与总生存期（OS）的关系
info <- summary(coxph(Surv(OS_month, OS_status) ~ ME11_pro_l, data = plot.cohort))

# 生成统计结果注释文本，包括风险比(HR)、置信区间(CI)和P值
anno.text <- ""
for (i in1:nrow(info$conf.int)) {
  anno.text <- paste0(anno.text, "\n", paste0(rownames(info$conf.int)[i], " HR=", round(info$conf.int[i, 1], 3), " CI=", round(info$conf.int[i, 3], 3), "-", round(info$conf.int[i, 4], 3), " P=", signif(info$coefficients[i, 5], 4) ))
}
# 添加Log-rank检验（Kaplan-Meier法）的P值
anno.text <- paste0(anno.text, "nKaplan-Meier P=", signif(survdiff(Surv(OS_month, OS_status) ~ ME11_pro_l, data = plot.cohort)$pvalue, 4) )
anno.text <- str_replace_all(anno.text, "ME11_pro_l", "") # 清理文本中的变量名

anno.text
Image

✅拟合Kaplan-Meier生存曲线

# 使用survfit函数拟合Kaplan-Meier生存曲线
fit <- survfit(Surv(OS_month, OS_status) ~ ME11_pro_l, data = plot.cohort)
fit
Image

✅绘制Kaplan-Meier生存曲线

# 使用ggsurvplot可视化fit
p1 <- ggsurvplot(fit, 
                 data = plot.cohort,
                 xlab = 'Time (Months)',       # X轴标签
                 pval = TRUE,                  # 在图上显示Log-rank检验的P值
                 risk.table = TRUE,           # 在曲线下方添加风险表，显示各时间点的风险患者数
                 risk.table.height = 0.28,    # 风险表的高度占比
                 conf.int.alpha = 0.05,      # 置信区间的透明度
                 conf.int = TRUE,             # 绘制生存曲线的置信区间
                 palette = c("#00599F","#d80700"), # 定义颜色，蓝色代表"Low"组，红色代表"High"组
                 axes.offset = TRUE,          # 坐标轴偏移
                 break.time.by = 12,          # X轴刻度间隔为12个月
                 xlim = c(0, 48),             # X轴范围限制在0到48个月
                 font.title = 12,
                 title= paste0("OS ME11_pro_l \n", anno.text)) # 图形标题，包含统计结果
p1 # 显示图形
Image

✅完整Kaplan-Meier生存曲线代码

# 设置图形输出参数：宽度8英寸，高度6英寸，分辨率200 DPI
options(repr.plot.width = 8, repr.plot.height = 6, repr.plot.res = 200)

# 加载生存分析和可视化所需的R包
library(survival)   # 用于生存分析，如Cox模型和生存曲线拟合
library(survminer)  # 用于绘制专业的生存曲线（Kaplan-Meier曲线）

#----------------------------------------------------------------------------------
#  Step 1: 加载数据
#----------------------------------------------------------------------------------
# (1) 读取191例PDAC患者（RJ队列1）的临床信息表
rj1.cohort <- read.xlsx("r_data/PDAC2_data_results/data/Extended Data Table 2.xlsx", startRow = 2)

# (2) 读取32个功能模块在癌组织（PDAC）和癌旁组织（TAC）中的ssGSEA评分矩阵，并进行转置（行为模块，列为样本）
module.ssgsea.pro.all <- read.xlsx("r_data/PDAC2_data_results/data/Extended Data Table 4.xlsx", sheet = 2, startRow = 2, rowNames = T)
module.ssgsea.pro.all <- t(module.ssgsea.pro.all) # 转置矩阵，使行名为模块名（如ME11），列名为样本ID

#----------------------------------------------------------------------------------
#  Step 2: 基于ME11模块进行总生存期（OS）分析
#----------------------------------------------------------------------------------
# 准备绘图数据：将ME11模块的评分添加到临床数据中，并根据其中位数将患者分为"High"/"Low"两组
plot.cohort <- rj1.cohort %>% 
  mutate(
    ME11_pro = as.numeric(module.ssgsea.pro.all["ME11", Proteomic_ID]), # 提取每个样本的ME11评分
    ME11_pro_l = ifelse(ME11_pro >= median(ME11_pro), "High", "Low"),    # 按中位数分组
    ME11_pro_l = factor(ME11_pro_l, levels = c("Low", "High")),         # 将分组设为因子，确定比较顺序
# 以下两行处理化疗状态因子，创建两种水平的因子用于可能的不同分析对比
    Censoring_chemo_rev = Censoring_chemo,
    Censoring_chemo = factor(Censoring_chemo, levels = c("Chem", "NoChem")),
    Censoring_chemo_rev = factor(Censoring_chemo_rev, levels = c("NoChem", "Chem"))
  )

# 执行Cox比例风险回归分析，检验ME11分组（High vs Low）与总生存期（OS）的关系
info <- summary(coxph(Surv(OS_month, OS_status) ~ ME11_pro_l, data = plot.cohort))

# 生成统计结果注释文本，包括风险比(HR)、置信区间(CI)和P值
anno.text <- ""
for (i in1:nrow(info$conf.int)) {
  anno.text <- paste0(anno.text, "n", paste0(rownames(info$conf.int)[i], " HR=", round(info$conf.int[i, 1], 3), " CI=", round(info$conf.int[i, 3], 3), "-", round(info$conf.int[i, 4], 3), " P=", signif(info$coefficients[i, 5], 4) ))
}
# 添加Log-rank检验（Kaplan-Meier法）的P值
anno.text <- paste0(anno.text, "\nKaplan-Meier P=", signif(survdiff(Surv(OS_month, OS_status) ~ ME11_pro_l, data = plot.cohort)$pvalue, 4) )
anno.text <- str_replace_all(anno.text, "ME11_pro_l", "") # 清理文本中的变量名

# 使用survfit函数拟合Kaplan-Meier生存曲线
fit <- survfit(Surv(OS_month, OS_status) ~ ME11_pro_l, data = plot.cohort)

# 使用ggsurvplot绘制生存曲线图，并设置多种可视化参数
p1 <- ggsurvplot(fit, 
                 data = plot.cohort,
                 xlab = 'Time (Months)',       # X轴标签
                 pval = TRUE,                  # 在图上显示Log-rank检验的P值
                 risk.table = TRUE,           # 在曲线下方添加风险表，显示各时间点的风险患者数
                 risk.table.height = 0.28,    # 风险表的高度占比
                 conf.int.alpha = 0.05,      # 置信区间的透明度
                 conf.int = TRUE,             # 绘制生存曲线的置信区间
                 palette = c("#00599F","#d80700"), # 定义颜色，蓝色代表"Low"组，红色代表"High"组
                 axes.offset = TRUE,          # 坐标轴偏移
                 break.time.by = 12,          # X轴刻度间隔为12个月
                 xlim = c(0, 48),             # X轴范围限制在0到48个月
                 font.title = 12,
                 title= paste0("OS ME11_pro_l n", anno.text)) # 图形标题，包含统计结果
p1 # 显示图形

本文结束！

图片
复现效果图-abcd图❤️
复现效果图-b图❤️
图片
❤️复现效果图-ABCD图❤️
图片
❤️复现效果图-bcd图❤️
图片
图片
❤️复现效果图-b图❤️
图片
图片
图片
图片
图片
图片
图片
图片
图片
图片
图片

图片

图片
图片
图片

（加入学习，收费，备注：299）

图片

KEEP THE THREAD GOING

更多来源文章