Skip to main content

Kaplan–Meier Survival Analysis

Generate Kaplan–Meier curve plots for Overall survival (OS) and progression-free survival (PFS)

Demo Data

The workflow follows the same logic as Figure_1.r in the reference repository (brucezhou9527/hello-world, branch Cell_Reports_Medicine_2025):

  1. A per-patient senescence score (in the original: ssGSEA of GOBP_CELLULAR_SENESCENCE) is dichotomised into High / Low.
  2. survfit(Surv(time, event) ~ status) fits the KM curves.
  3. survminer::ggsurvplot() draws them (palette #E93E3F for High, #A6CEE3 for Low, x-axis 0–24 months in steps of 6, y-axis in percent).

Because the real clinical and expression data are not public here, We simulate demo data whose group sizes, curve shapes and log-rank p-values reproduce the figure.

library(survival) # Surv(), survfit(), survdiff()
library(survminer) # ggsurvplot(), surv_pvalue()
library(ggplot2)
library(patchwork) # combine the two panels
library(tidyverse) # dplyr, tibble, readr

col_high <- "#E93E3F" # red
col_low <- "#A6CEE3" # light blue

# target survival curves read from the figure
target <- list(
    OS = list(
        High = list(
            t = c(0, 2, 4, 6, 8, 10, 12, 14, 18, 24),
            S = c(1, .97, .83, .69, .53, .44, .36, .30, .27, .25)
        ),
        Low = list(
            t = c(0, 2, 4, 6, 9, 12, 14, 18, 21, 24),
            S = c(1, .98, .88, .80, .66, .58, .53, .46, .37, .36)
        )
    ),
    PFS = list(
        # steep early drop around 2-3 months (first imaging)
        High = list(
            t = c(0, 1.5, 1.9, 2.2, 2.5, 3, 5, 9, 15, 24),
            S = c(1, .97, .90, .66, .46, .33, .23, .17, .12, .11)
        ),
        Low = list(
            t = c(0, 1.5, 2, 2.3, 2.6, 3.2, 4.5, 7, 12, 18, 24),
            S = c(1, .97, .90, .70, .54, .50, .36, .30, .23, .19, .18)
        )
    )
)

# inverse of a piece-wise linear survival curve
inv_surv <- function(u, tg, tail_rate = 1 / 12) {
    f <- approxfun(rev(tg$S), rev(tg$t), ties = "ordered") # S -> t
    smin <- min(tg$S)
    tmax <- max(tg$t)
    # u below the plateau: event would happen after 24 months (will be censored)
    ifelse(u > smin, f(pmin(u, 1)), tmax + rexp(length(u), tail_rate))
}

simulate_cohort <- function(seed, n_high = 120, n_low = 88) {
    set.seed(seed)
    n <- n_high + n_low
    score <- rnorm(n, mean = 0.35, sd = 0.08) # ssGSEA-like score
    ord <- order(score)
    grp <- rep("High", n)
    grp[ord[seq_len(n_low)]] <- "Low"

    u <- runif(n) # shared by OS & PFS
    cens <- pmin(6 + rexp(n, 1 / 22), 24) # random censoring, max 24 mo

    os_lat <- pfs_lat <- numeric(n)
    for (g in c("High", "Low")) {
        i <- grp == g
        os_lat[i] <- inv_surv(u[i], target$OS[[g]])
        pfs_lat[i] <- inv_surv(u[i], target$PFS[[g]])
    }

    data.frame(
        patient_id = sprintf("P%03d", seq_len(n)),
        CANCER_TYPE = "mUC",
        GOBP_CELLULAR_SENESCENCE = round(score, 6),
        OS_MONTHS = round(pmin(os_lat, cens), 2),
        OS_CNSR_REV = as.integer(os_lat <= cens), # 1 = death observed
        PFS_MONTHS = round(pmin(pfs_lat, cens), 2),
        PFS_CNSR_REV = as.integer(pfs_lat <= cens) # 1 = progression/death
    )
}

demo <- simulate_cohort(seed = 1967) |> as_tibble()

# write.csv(demo, "demo_survival_data.csv", row.names = FALSE)
head(demo)
## # A tibble: 6 × 7
##   patient_id CANCER_TYPE GOBP_CELLULAR_SENESC…¹ OS_MONTHS OS_CNSR_REV PFS_MONTHS
##   <chr>      <chr>                        <dbl>     <dbl>       <int>      <dbl>
## 1 P001       mUC                          0.394      3.3            1       1.93
## 2 P002       mUC                          0.362      5.6            1       2.13
## 3 P003       mUC                          0.275     15.2            0       3.91
## 4 P004       mUC                          0.316     24              0      16.2 
## 5 P005       mUC                          0.359      1.37           1       1.03
## 6 P006       mUC                          0.445     14.9            1       3.74
## # ℹ abbreviated name: ¹​GOBP_CELLULAR_SENESCENCE
## # ℹ 1 more variable: PFS_CNSR_REV <int>

Dichotomise the score

Same logic as Figure_1.r (status <- ifelse(score > cutoff, "High", "Low")), but with a fixed cut-point that gives 120 / 88 patients. To use a data-driven optimal cut-point instead, replace cutoff with the result of survminer::surv_cutpoint() (see the commented lines).

km_data <- subset(demo, CANCER_TYPE == "mUC")

sorted <- sort(km_data$GOBP_CELLULAR_SENESCENCE)

# 88 lowest scores -> "Low"
cutoff <- mean(sorted[88:89])

# as in Figure_1.r
# cutoff <- median(km_data$GOBP_CELLULAR_SENESCENCE)
# cp <- surv_cutpoint(km_data, "OS_MONTHS", "OS_CNSR_REV", "GOBP_CELLULAR_SENESCENCE")
# cutoff <- cp$cutpoint$cutpoint                                # optimal cut-point

km_data$status <- factor(
    ifelse(km_data$GOBP_CELLULAR_SENESCENCE > cutoff, "High", "Low"),
    levels = c("High", "Low")
)
table(km_data$status)
## 
## High  Low 
##  120   88

Kaplan–Meier Fits

fit_os <- survfit(Surv(OS_MONTHS, OS_CNSR_REV) ~ status, data = km_data)
fit_pfs <- survfit(Surv(PFS_MONTHS, PFS_CNSR_REV) ~ status, data = km_data)

# log-rank tests
p_os <- surv_pvalue(fit_os, data = km_data)$pval
p_pfs <- surv_pvalue(fit_pfs, data = km_data)$pval
c(OS = p_os, PFS = p_pfs)
##          OS         PFS 
## 0.004938708 0.021325871

Survival Plots

ggsurvplot() supplies the curves, the censor ticks and the percent axis. Everything else that makes the figure look like the original is added with ggplot2:

  • a plain theme_classic() (only left and bottom axis lines),
  • the p-value as a text annotation (not the default pval = TRUE position),
  • a hand-drawn legend (hollow coloured squares + text, with the title Cellular Senescence Score), placed in data coordinates; clip = "off" lets it sit on the axis edge.
km_panel <- function(fit, data, ylab, pval, p_xy, leg) {
    # leg: list(tx, ty, ey, hx, lx)  title x/y, entry y, High x, Low x
    g <- ggsurvplot(
        fit,
        data = data,
        surv.scale = "percent",
        break.time.by = 6,
        xlim = c(0, 24),
        linewidth = 0.9,
        censor.shape = "+",
        censor.size = 2.6,
        linetype = c(1, 1),
        palette = c(col_high, col_low),
        xlab = "Months after ICB",
        ylab = ylab,
        pval = FALSE,
        conf.int = FALSE,
        legend = "none",
        ggtheme = theme_classic(base_size = 12)
    )$plot

    n_h <- sum(data$status == "High")
    n_l <- sum(data$status == "Low")
    bw <- 1.1
    bh <- 0.065 # legend box width / height

    g +
        coord_cartesian(xlim = c(0, 24), ylim = c(0, 1), clip = "off") +
        scale_x_continuous(
            breaks = seq(0, 24, 6),
            expand = expansion(mult = c(0.02, 0.02))
        ) +
        annotate(
            "text",
            x = p_xy[1],
            y = p_xy[2],
            hjust = 0,
            size = 4.6,
            label = sprintf("p = %s", pval)
        ) +
        # legend title
        annotate(
            "text",
            x = leg$tx,
            y = leg$ty,
            hjust = 0,
            size = 4,
            label = "Cellular Senescence Score"
        ) +
        # High
        annotate(
            "rect",
            xmin = leg$hx,
            xmax = leg$hx + bw,
            ymin = leg$ey - bh / 2,
            ymax = leg$ey + bh / 2,
            fill = NA,
            colour = col_high,
            linewidth = 0.8
        ) +
        annotate(
            "text",
            x = leg$hx + bw + 0.3,
            y = leg$ey,
            hjust = 0,
            size = 3.6,
            label = sprintf("High (n=%d)", n_h)
        ) +
        # Low
        annotate(
            "rect",
            xmin = leg$lx,
            xmax = leg$lx + bw,
            ymin = leg$ey - bh / 2,
            ymax = leg$ey + bh / 2,
            fill = NA,
            colour = col_low,
            linewidth = 0.8
        ) +
        annotate(
            "text",
            x = leg$lx + bw + 0.3,
            y = leg$ey,
            hjust = 0,
            size = 3.6,
            label = sprintf("Low (n=%d)", n_l)
        ) +
        theme(
            axis.text = element_text(colour = "black"),
            axis.title = element_text(colour = "black"),
            axis.line = element_line(colour = "black")
        )
}

fmt_p <- function(p) formatC(p, digits = if (p < 0.01) 3 else 3, format = "g")

g_os <- km_panel(
    fit_os,
    km_data,
    "Overall survival probability",
    fmt_p(round(p_os, 3)),
    p_xy = c(9.7, 0.79),
    leg = list(tx = 1.4, ty = 0.19, ey = 0.05, hx = -0.3, lx = 10.5)
)

g_pfs <- km_panel(
    fit_pfs,
    km_data,
    "Progression-free survival",
    fmt_p(round(p_pfs, 3)),
    p_xy = c(4.8, 0.85),
    leg = list(tx = 6.4, ty = 0.57, ey = 0.44, hx = 4.9, lx = 15.6)
)

g_os | g_pfs

Kaplan–Meier curves of overall survival (left) and progression-free survival (right) by Cellular Senescence Score (demo data).

# fig <- (g_os | g_pfs)
# ggsave("KM_OS_PFS.png", fig, width = 11, height = 4.2, dpi = 300)
# ggsave("KM_OS_PFS.pdf", fig, width = 11, height = 4.2)
# tiff("KM_OS_PFS.tiff", units = "in", width = 11, height = 4.2, res = 300, compression = "lzw"); print(fig); dev.off()

Session Info

sessionInfo()
## R version 4.5.3 (2026-03-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Asia/Shanghai
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lubridate_1.9.5 forcats_1.0.1   stringr_1.6.0   dplyr_1.2.1    
##  [5] purrr_1.2.2     readr_2.2.0     tidyr_1.3.2     tibble_3.3.1   
##  [9] tidyverse_2.0.0 patchwork_1.3.2 survminer_0.5.2 ggpubr_1.0.0   
## [13] ggplot2_4.0.3   survival_3.8-12
## 
## loaded via a namespace (and not attached):
##  [1] utf8_1.2.6         generics_0.1.4     rstatix_1.1.0      stringi_1.8.9     
##  [5] lattice_0.23-1     hms_1.1.4          digest_0.6.39      magrittr_2.0.5    
##  [9] timechange_0.4.0   evaluate_1.0.5     grid_4.5.3         RColorBrewer_1.1-3
## [13] fastmap_1.2.0      jsonlite_2.0.0     Matrix_1.7-5       backports_1.5.1   
## [17] Formula_1.2-6      gridExtra_2.3.1    scales_1.4.0       abind_1.4-8       
## [21] cli_3.6.6          rlang_1.3.0        splines_4.5.3      withr_3.0.3       
## [25] yaml_2.3.12        otel_0.2.0         tools_4.5.3        parallel_4.5.3    
## [29] tzdb_0.5.0         ggsignif_0.6.4     broom_1.0.13       vctrs_0.7.3       
## [33] R6_2.6.1           lifecycle_1.0.5    car_3.1-5          htmlwidgets_1.6.4 
## [37] pkgconfig_2.0.3    pillar_1.11.1      gtable_0.3.6       glue_1.8.1        
## [41] xfun_0.61          tidyselect_1.2.1   knitr_1.52         farver_2.1.2      
## [45] htmltools_0.5.9    labeling_0.4.3     rmarkdown_2.32     carData_3.0-6     
## [49] compiler_4.5.3     S7_0.2.2
Back to top