| name | bulkrna-survival |
| description | Load when stratifying patients by gene expression and testing for survival differences (Kaplan-Meier + Cox) in bulk RNA-seq. Skip if no time-to-event clinical data exists, or for non-bulk cohorts (single-cell / spatial survival is not supported). |
| version | 0.3.0 |
| author | OmicsClaw |
| license | MIT |
| tags | ["bulkrna","survival","Kaplan-Meier","Cox","hazard-ratio","clinical"] |
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
| Input | Format | Required |
|---|
| Expression matrix | .csv (gene × sample) | yes (or --demo) |
| Clinical data | .csv (sample, time, event cols) via --clinical | yes (or --demo) |
| Genes to test | --genes TP53,BRCA1 (comma-separated) | optional, defaults to all in expression matrix |
| Output | Path | Notes |
|---|
| Per-gene survival stats | tables/survival_results.csv | columns: gene, cutoff, n_high, n_low, hazard_ratio, log_rank_chi2, log_rank_pval, median_survival_high, median_survival_low |
| KM curves | figures/km_<gene>.png | one per gene tested |
| Forest plot | figures/forest_plot.png | HR + CI across genes |
| Report | report.md + result.json | summary contains n_genes and per-gene results array (bulkrna_survival.py:153-160) |
Flow
- Load expression matrix + clinical data; align by sample ID.
- 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%.
- Try R
survival package first; fall back to Python lifelines (:626 warns "R survival not available (...); using Python fallback.").
- 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.
Key CLI
python omicsclaw.py run bulkrna-survival --demo
python omicsclaw.py run bulkrna-survival \
--input expression.csv --clinical clinical.csv \
--genes TP53,BRCA1,EGFR --output results/
python omicsclaw.py run bulkrna-survival \
--input expression.csv --clinical clinical.csv \
--genes TP53 --cutoff-method optimal --output results/
See also
references/parameters.md — every CLI flag and tuning hintreferences/methodology.md — KM + log-rank + Cox theory, R vs Python backend differences, optimal-cutoff caveatsreferences/output_contract.md — exact output directory layout- Adjacent skills:
bulkrna-de (parallel: differential expression — survival adds the time-to-event dimension), bulkrna-coexpression (parallel: module-level survival via eigengene if traits include time-to-event)
