bulkrna-survival

bulkrna-survival

When to use

Run on a bulk RNA-seq cohort with paired clinical survival data (time-to-event + censoring) when you want to ask "does high vs low expression of gene X predict survival?". Default workflow: per-gene median-cutoff stratification, log-rank p-value, Kaplan-Meier curve, and Cox proportional-hazards hazard ratio.

Inputs & Outputs

Inputs

• File types: .csv

Outputs

• tables/clinical.csv
• tables/expr.csv
• tables/km_data.csv
• tables/survival_results.csv
• figures/forest_plot.png
• report.md
• result.json

Flow

1. Load expression matrix + clinical data; align by sample ID.
2. For each gene in --genes (or all):
   • Skip with warning at bulkrna_survival.py:630 if gene not in expression matrix.
   • Stratify samples by --cutoff-method (default median; alt optimal finds the maxstat cut).
   • Run log-rank test on the stratified groups.
   • Compute a simple events/time hazard ratio. Warn at :326 ("Heavy censoring (X%). KM tail estimates may be unreliable.") when the censoring rate exceeds 80%.
3. Try R survival package first; fall back to Python lifelines (:626 warns "R survival not available (...); using Python fallback.").
4. Render KM curves + forest plot; emit tables/survival_results.csv.

Gotchas

• Genes not in the expression matrix are silently skipped. bulkrna_survival.py:630 logs a warning per missing gene and continues. After the run, count the rows in tables/survival_results.csv (or inspect result.json["results"]) and compare against the --genes list — a typo'd or wrong-namespace gene produces no obvious error.
• --cutoff-method optimal p-values are NOT corrected for multiple testing. The optimal cutoff scans all possible cuts and picks the maximally separating one, which inflates Type I error. Reported log-rank p-values are raw — apply Bonferroni / BH correction externally if you scan many genes.
• The hazard ratio is a simple events/person-time ratio, not a Cox MLE. The script computes (events_high / time_high) / (events_low / time_low) (bulkrna_survival.py:328-333), not a Cox proportional-hazards regression coefficient. This estimator is biased when proportional-hazards holds with unequal exposure — for publication-grade HRs, re-fit a proper Cox model in R or lifelines against the same stratification.
• R-vs-Python backend silently switches. :626 warns and falls back to a NumPy log-rank implementation when R survival isn't importable; the per-gene HR estimator is the same simple events/time ratio in both cases, but the chosen backend isn't recorded in the summary dict — only in the warning log. Verify R availability before relying on the result for downstream papers.
• Heavy censoring distorts KM tail estimates. :326 fires when ≥80% of patients are censored; the printed median survival numbers are dominated by extrapolation past the last event time. Treat median_survival_* as "≥ X" rather than a point estimate when the corresponding gene's censoring rate is high.