Skip to contents

Kaplan-Meier survival plots for bulk expression data

SurvPlot.Rd

Generates publication-ready Kaplan-Meier survival plots for one or more genes or gene-set signatures against bulk RNA-seq or microarray expression data. The function has two input modes:

  1. TCGA mode — supply tcga_project and the function downloads, processes, and caches the expression matrix and clinical metadata automatically via TCGAbiolinks. Subsequent calls with the same project load from the local cache instantly.

  2. Custom matrix mode — supply your own expr_mat (genes × samples) and meta (samples × variables) data frames. This is the path for non-TCGA cohorts such as the Cavalli medulloblastoma dataset, PBTA, or any in-house bulk RNA-seq study.

Samples are classified as High / Low (or Q1–Q4) using the chosen split_method, and KM curves are drawn with log-rank p-values. Results can be faceted across two clinical variables simultaneously (e.g. sex as rows, tumour subtype as columns) producing one grid page per gene/signature. A sidecar CSV records every cutoff, group size, and p-value for full reproducibility.

Usage

SurvPlot(
  expr_mat = NULL,
  meta = NULL,
  tcga_project = NULL,
  tcga_data_type = c("rnaseq", "microarray"),
  tcga_backend = c("recount3", "TCGAbiolinks"),
  tcga_cache_dir = NULL,
  genes = NULL,
  gene_lists = NULL,
  time_col = "OS",
  status_col = "status",
  split_method = c("median", "quartile", "value", "optimal"),
  split_value = NULL,
  split_per_group = TRUE,
  facet_row = NULL,
  facet_col = NULL,
  output_file = NULL,
  pdf_width = NULL,
  pdf_height = NULL,
  palette = c("Red", "Blue"),
  quartile_palette = c("#2166AC", "#92C5DE", "#F4A582", "#D6604D"),
  show_pval = TRUE,
  show_n = TRUE,
  write_sidecar = TRUE,
  min_group_n = 5L,
  multivariate_covariates = NULL
)

Arguments

expr_mat

Numeric matrix, genes × samples. Rownames must be gene symbols; colnames must be sample identifiers matching rownames(meta). Values should be on a log scale (e.g. log2 TPM + 1 or log2 FPKM + 1). Supply either this + meta, or tcga_project.

meta

Data frame of sample-level clinical metadata. Rownames must match colnames(expr_mat). Must contain at minimum the columns named by time_col and status_col.

tcga_project

Character scalar. TCGA (or TARGET) project identifier. See the TCGA projects section above for the full list of available IDs. Requires TCGAbiolinks to be installed: BiocManager::install("TCGAbiolinks"). Data are downloaded once and cached as an .rds file; subsequent calls load from the cache.

tcga_data_type

Data modality to download from GDC. One of:

"rnaseq" (default)

STAR-aligned RNA-seq counts (workflow.type = "STAR - Counts"). TPM values are extracted and log2-transformed (log2(TPM + 1)). Falls back to FPKM if TPM is unavailable.

"microarray"

Legacy microarray data via the GDC legacy archive. Use this for older TCGA accessions that pre-date RNA-seq, or for external cohorts processed on Affymetrix / Illumina arrays (e.g. the Cavalli medulloblastoma dataset).

tcga_backend

Download engine for TCGA / TARGET data. One of:

"recount3" (default)

Downloads one pre-consolidated genes × samples count matrix per project — a single file, no per-patient downloads. Counts are converted to log2(TPM + 1). Requires recount3: BiocManager::install("recount3"). Automatically switches to "TCGAbiolinks" when tcga_data_type = "microarray".

"TCGAbiolinks"

Downloads one file per patient from the GDC API (e.g. 391 files for TCGA-GBM). Much slower and prone to partial-download crashes. Use only when recount3 does not carry the project, or for microarray data.

tcga_cache_dir

Path to a directory where downloaded data are cached as .rds files. Defaults to the OS user cache directory for scSidekick (tools::R_user_dir("scSidekick", "cache")). Set a shared network path to allow multiple users to reuse the same cache. Delete the .rds file to force a fresh download.

genes

Character vector of gene symbols to analyze in single-gene mode. Each gene produces one set of KM panels. Genes absent from expr_mat trigger a warning and are skipped.

gene_lists

Named list of character vectors for signature-score mode. Each element is a gene set; the name becomes the panel title. Scores are the mean log2 expression of genes present in the matrix. Example: list(YAP_up = c("PRC1","HPCA"), YAP_dn = c("ZIC1","CHD7")).

time_col

Name of the survival time column in meta (or in the standardised TCGA metadata). For TCGA downloads this is "OS" (overall survival in days). For disease-specific survival use the column produced by your preprocessing pipeline.

status_col

Name of the event indicator column. Must be numeric: 1 = event occurred (death / recurrence), 0 = censored. For TCGA downloads this is "status" (standardised from vital_status by SurvPlot automatically).

split_method

Strategy for classifying samples into expression groups. One of "median" (default), "quartile", "value", or "optimal". See the Splitting strategies section for full details of each option.

split_value

Numeric cutoff used when split_method = "value". Samples with expression / score \(\geq\) split_value are labeled High; all others are Low. Commonly set to 0 for gene-set enrichment scores.

split_per_group

Logical (default TRUE). When TRUE, the threshold (median / quantiles) is computed independently within each facet group. This ensures a 50/50 split within every subtype or sex strata, preventing a globally-low subgroup from being assigned entirely to "Low". Set to FALSE to use a single global threshold across all samples (equivalent to not splitting per group).

facet_row

Name of a column in meta whose levels define the rows of the output grid. Common choices: "gender", "race", "age_group". Set to NULL (default) for no row faceting.

facet_col

Name of a column in meta whose levels define the columns of the output grid. Common choices: "Subgroup", "paper_IDH.status", "tumor_stage". Set to NULL (default) for no column faceting.

output_file

Output path. Two modes depending on the path:

Ends in .pdf

All genes are written as pages into a single multi-page PDF.

Does not end in .pdf

Treated as a directory. One PDF per gene/signature is written (<gene>_survival.pdf), plus one shared survival_sidecar.csv. The directory is created if it does not exist.

Set to NULL (default) to suppress file output and return plots invisibly for further manipulation.

pdf_width, pdf_height

Numeric. PDF dimensions in inches. When NULL (default) dimensions are auto-calculated from the number of genes and facet levels: approximately 4.5 in per column and 4.5 in per row.

palette

Character vector of length \(\geq 2\). Colors for the High and Low groups respectively. Default: c("Red", "Blue"). Any R color name or hex code is accepted.

quartile_palette

Character vector of length \(\geq 4\). Colors for Q1 through Q4 when split_method = "quartile". Default is a blue-to-red diverging scale: c("#2166AC", "#92C5DE", "#F4A582", "#D6604D").

show_pval

Logical (default TRUE). Display the log-rank p-value on each KM panel.

show_n

Logical (default TRUE). Append the group sample size to each legend label, e.g. "High (n = 47)".

write_sidecar

Logical (default TRUE). Write a companion CSV file recording the cutoff, group sizes, split method, and log-rank p-value for every panel. Ignored when output_file = NULL. See the Sidecar CSV section for column definitions.

min_group_n

Integer (default 5). Minimum number of samples required in each group after splitting. Panels where any group falls below this threshold are skipped and recorded as NA in the sidecar.

multivariate_covariates

Character vector of meta column names to include as covariates in a multivariate Cox proportional hazards model. When supplied, a forest plot (survminer::ggforest) is printed as an additional page after the KM panels for each gene / signature. The Cox model uses continuous expression values (not High/Low labels) as the feature term, which is more powerful than treating it as a binary predictor. All samples with complete data across the expression feature and covariates are used (no faceting for Cox). Example: c("gender", "age_at_index", "paper_IDH.status"). Set to NULL (default) to skip multivariate analysis.

Value

Invisibly returns a named list with two elements:

plots

Named list, one entry per gene/signature. Each entry is itself a list of ggsurvplot objects (one per facet combination), or NULL for skipped panels.

sidecar

Data frame with one row per KM panel documenting the feature name, facet levels, cutoff, group sizes, and log-rank p-value.

TCGA projects

Pass any of the following strings to tcga_project. All IDs follow the format "TCGA-<CODE>".

Brain / CNS

"TCGA-GBM"

Glioblastoma multiforme (n ≈ 173)

"TCGA-LGG"

Brain Lower Grade Glioma (n ≈ 516)

Breast

"TCGA-BRCA"

Breast invasive carcinoma (n ≈ 1 098)

Lung

"TCGA-LUAD"

Lung adenocarcinoma (n ≈ 585)

"TCGA-LUSC"

Lung squamous cell carcinoma (n ≈ 504)

Colorectal

"TCGA-COAD"

Colon adenocarcinoma (n ≈ 521)

"TCGA-READ"

Rectum adenocarcinoma (n ≈ 177)

Liver / Biliary

"TCGA-LIHC"

Liver hepatocellular carcinoma (n ≈ 377)

"TCGA-CHOL"

Cholangiocarcinoma (n ≈ 51)

Kidney

"TCGA-KIRC"

Kidney renal clear cell carcinoma (n ≈ 611)

"TCGA-KIRP"

Kidney renal papillary cell carcinoma (n ≈ 321)

"TCGA-KICH"

Kidney chromophobe (n ≈ 113)

Gynaecological

"TCGA-OV"

Ovarian serous cystadenocarcinoma (n ≈ 608)

"TCGA-UCEC"

Uterine corpus endometrial carcinoma (n ≈ 587)

"TCGA-UCS"

Uterine carcinosarcoma (n ≈ 57)

"TCGA-CESC"

Cervical squamous cell carcinoma (n ≈ 309)

Skin

"TCGA-SKCM"

Skin cutaneous melanoma (n ≈ 473)

"TCGA-UVM"

Uveal melanoma (n ≈ 80)

Gastrointestinal / Pancreatic

"TCGA-STAD"

Stomach adenocarcinoma (n ≈ 443)

"TCGA-ESCA"

Esophageal carcinoma (n ≈ 185)

"TCGA-PAAD"

Pancreatic adenocarcinoma (n ≈ 185)

Thoracic

"TCGA-MESO"

Mesothelioma (n ≈ 87)

"TCGA-THYM"

Thymoma (n ≈ 124)

Head and Neck / Thyroid

"TCGA-HNSC"

Head and neck squamous cell carcinoma (n ≈ 528)

"TCGA-THCA"

Thyroid carcinoma (n ≈ 507)

Haematological

"TCGA-LAML"

Acute myeloid leukaemia (n ≈ 200)

"TCGA-DLBC"

Diffuse large B-cell lymphoma (n ≈ 48)

Sarcoma / Mesenchymal

"TCGA-SARC"

Sarcoma (n ≈ 265)

Urological

"TCGA-BLCA"

Bladder urothelial carcinoma (n ≈ 433)

"TCGA-PRAD"

Prostate adenocarcinoma (n ≈ 551)

"TCGA-TGCT"

Testicular germ cell tumours (n ≈ 156)

Endocrine

"TCGA-ACC"

Adrenocortical carcinoma (n ≈ 92)

"TCGA-PCPG"

Pheochromocytoma and paraganglioma (n ≈ 187)

Paediatric (TARGET programme)

"TARGET-AML"

Acute myeloid leukaemia

"TARGET-NBL"

Neuroblastoma

"TARGET-OS"

Osteosarcoma

"TARGET-WT"

Wilms tumour

"TARGET-RT"

Rhabdoid tumour

"TARGET-CCSK"

Clear cell sarcoma of the kidney

Standard TCGA clinical columns

After download, SurvPlot standardises OS (overall survival in days, combining days_to_death and days_to_last_follow_up) and status (1 = deceased, 0 = living/censored). Common columns available for facet_row / facet_col include:

gender

Reported sex ("male" / "female")

race

Self-reported race

age_at_index

Age at diagnosis (years)

tumor_stage / pathologic_stage

Pathological stage

paper_*

TCGA marker-paper columns, e.g. paper_IDH.status, paper_Histology, paper_CIMP.status, paper_Pan.Glioma.RNA.Expression.Cluster

Available columns vary by cancer type. To inspect what is available after the first download run: colnames(result$sidecar) or examine meta from the returned list.

Splitting strategies

"median" (default)

Splits samples at the median expression / score value into High and Low groups. When split_per_group = TRUE (recommended), the median is computed independently within each facet cell, so that a subgroup with uniformly low expression is still split 50/50 — matching the approach used in most published survival analyses.

"quartile"

Divides samples into four equally-sized groups (Q1 = lowest 25\ four KM curves per panel colored by quartile_palette. Useful for detecting non-linear dose–response relationships.

"value"

Splits at a user-supplied fixed numeric threshold (split_value). Applies the same cutoff across all groups and facets. Set split_value = 0 for gene-set enrichment scores computed by gene_lists mode (positive score = enriched).

"optimal"

Uses survminer::surv_cutpoint() to find the threshold that maximises the log-rank statistic. Computationally intensive and prone to overfitting on small cohorts; not recommended for discovery. Falls back to "median" if the optimisation fails. Requires survminer \(\geq\) 0.4.9.

Gene-set signature scoring

When gene_lists is supplied instead of genes, the function computes a per-sample score for each gene set as the mean log2 expression across all genes in the set that are present in the matrix. This is the same mean-expression scorer used internally by PlotPathwayButterfly; no additional packages are required. Missing genes are silently dropped and a warning is issued if more than 50\

For enrichment-style interpretation, use split_method = "value" with split_value = 0 (score \(\geq\) 0 = enriched) or split_method = "median" (top half = enriched).

Output layout

No faceting

All genes / signatures appear on a single row in one page (or one PDF file per gene in directory mode).

facet_row only

Panels are stacked vertically, one row per level of the row variable.

facet_row + facet_col

One page per gene/signature. Panels are arranged in an (n\_row\_levels × n\_col\_levels) grid. Cells with fewer than 2 × min_group_n samples are left blank to preserve alignment. PDF dimensions are auto-scaled to the grid size unless overridden with pdf_width / pdf_height.

Output file path

Single PDF

Set output_file to a path ending in .pdf (e.g. "~/results/survival.pdf"). All genes appear as separate pages in one file.

Directory

Set output_file to a folder path (existing or new; must not end in .pdf). One PDF per gene/signature is written, named <gene>_survival.pdf. A single survival_sidecar.csv covering all genes is written to the same folder.

Sidecar CSV

Written alongside every PDF output, the sidecar records one row per KM panel with the following columns:

feature

Gene name or signature name

facet_row, facet_col

Level values ("all" when no faceting)

split_method

Method actually used (may differ from request if fallback occurred)

cutoff

Threshold applied; for "quartile" this is a semicolon-separated list of the five quantile boundaries

n_total

Total samples in the subset

n_High, n_Low (or n_Q1–n_Q4)

Samples per group

pval

Log-rank p-value

See also

PlotPathwayButterfly for single-cell pathway state scoring using the same mean-expression scorer.

TCGAbiolinks for the GDC download engine used in TCGA mode.

survminer for the underlying ggsurvplot plotting engine.

Examples

if (FALSE) { # \dontrun{
# ── Example 1: TCGA-GBM, two genes, sex × IDH status grid ────────────────
SurvPlot(
  tcga_project = "TCGA-GBM",
  genes        = c("YAP1", "S100A4"),
  time_col     = "OS",
  status_col   = "status",
  facet_row    = "gender",
  facet_col    = "paper_IDH.status",
  output_file  = "~/Desktop/GBM_survival/"   # directory → one PDF per gene
)

# ── Example 2: TCGA-GBM, quartile split, all in one PDF ──────────────────
SurvPlot(
  tcga_project = "TCGA-GBM",
  genes        = c("YAP1", "CD276", "CD47"),
  split_method = "quartile",
  output_file  = "~/Desktop/GBM_quartiles.pdf"
)

# ── Example 3: Gene-set signature scores, value split at 0 ───────────────
SurvPlot(
  tcga_project = "TCGA-LGG",
  gene_lists   = list(
    YAP_Activated = c("PRC1", "HPCA", "SRSF7", "TIMM17A", "PIM1"),
    YAP_Repressed = c("ZIC1", "ZIC2", "CHD7", "KIF5C", "DPYSL2")
  ),
  split_method  = "value",
  split_value   = 0,       # >= 0 = enriched
  facet_row     = "gender",
  facet_col     = "paper_Pan.Glioma.RNA.Expression.Cluster",
  output_file   = "~/Desktop/LGG_signatures/"
)

# ── Example 4: Custom matrix (e.g. Cavalli medulloblastoma) ──────────────
Cav.EXP  <- t(read.table("Cavalli_expression.txt", header = TRUE,
                           row.names = 1))
Cav.meta <- read.table("Cavalli_pheno.txt", header = TRUE, row.names = 1)
# Cavalli metadata uses "survival" and "status" already
res <- SurvPlot(
  expr_mat     = Cav.EXP,
  meta         = Cav.meta,
  genes        = c("YAP1", "CD276", "CD47", "CCL5"),
  time_col     = "survival",
  status_col   = "status",
  facet_row    = "Gender",
  facet_col    = "Subgroup",
  split_method = "median",
  output_file  = "~/Desktop/Cavalli_survival/"
)
# Inspect the sidecar
head(res$sidecar)

# ── Example 5: Paediatric cancer via TARGET ───────────────────────────────
SurvPlot(
  tcga_project  = "TARGET-NBL",
  genes         = c("MYCN", "ALK"),
  time_col      = "OS",
  status_col    = "status",
  split_method  = "optimal",
  output_file   = "~/Desktop/NBL_survival.pdf"
)

# ── Example 6: No faceting, return plots for manual arrangement ───────────
res <- SurvPlot(
  tcga_project = "TCGA-SKCM",
  genes        = "S100A4",
  split_method = "median"
)
# Access the ggsurvplot object for the first panel
print(res$plots[["S100A4"]][[1]])
} # }