SIDISH — Cell-Death / High-Risk Cell Identification

Overview

What it does. SIDISH (Semi-supervised Iterative Deep Learning for Identifying Single-cell High-Risk Populations) integrates three data types — scRNA-seq, bulk RNA-seq and patient survival — into one iterative deep-learning loop:

• Phase 1 trains a VAE on the single-cell expression matrix.
• Phase 2 trains a Deep Cox survival model on bulk RNA-seq, using the VAE latent space to score cells with the patient-level hazard model.
• Between iterations, cell weights and patient weights are updated (SHAP-based feature selection + a Weibull/fitted hazard distribution), and cells above the chosen hazard percentile are relabelled High-Risk (h) vs Background (b).

Its downstream module performs in-silico perturbation: knock a gene (or a gene pair) out, re-run the hazard model, and count how many High-Risk cells flip back to Background. That gives a ranked list of candidate targets.

When to use it.
- You have scRNA-seq (or spatial) plus a bulk RNA-seq cohort with survival data (duration + event), and you want cell populations whose transcriptional state tracks the worst clinical outcome.
- You want a therapeutic-target shortlist from perturbation scoring, or a risk-associated cell-state marker-gene list to test on the bulk cohort (Kaplan-Meier).

When NOT to use it.
- No survival data / no bulk cohort → the whole Phase 2 signal is missing; use a death-signature or apoptosis-scoring approach instead.
- You literally need dead/dying cells (Annexin V, Caspase, cell-cycle arrest, pct MT): SIDISH has no cell-death score in its objective. It labels cells that are survival-linked high-risk; interpreting them as "dying/damaged cells" is a biological interpretation that must be validated (see Pitfalls).

Upstream.
- Repo: https://github.com/mcgilldinglab/SIDISH (verified HTTP 200)
- Docs/API: https://sidish.readthedocs.io/en/latest/api.html (verified HTTP 200)
- PyPI: SIDISH==1.0.0 (https://pypi.org/pypi/SIDISH/json verified 200, only release 1.0.0)
- Paper-style description and tutorials (4 notebooks) live in the repo tutorials/.

This skill's script. scripts/sidish_cell_death.py is a single command-line wrapper that reproduces the reference notebook chain (script/SIDISH.ipynb in the user's copy, equivalent to upstream tutorials/tutorial_0..3) with every path as a CLI flag and no hardcoded directories.

Inputs

Verified against the reference dataset shipped with the local copy of the tool (LUAD / TCGA-LUAD):

Single-cell .obs columns the pipeline reads or writes: cells, n_genes, scissors, scAB, degas, SIDISH_value, risk_value, SIDISH (labels h/b). preprocess(processed=True) optionally renames a cell-type column to celltype_major; the reference data has no cell-type column, so celltype_name is a no-op there. var columns: gene_ids, n_cells — preprocess() prefers var['gene_ids'] over var_names when intersecting with bulk genes.

Data you must bring for a new project: (1) an .h5ad whose var_names (or var['gene_ids']) are gene symbols, (2) a bulk matrix of the same gene space, (3) sample-level survival with 1/0 status whose sample IDs match the bulk columns (the notebook normalises TCGA.97.7938.01A → TCGA-97-7938-01A by replacing . with -).

Steps

Script: scripts/sidish_cell_death.py (all flags, no hardcoded paths). Run with the SIDISH environment's python on a GPU machine.

1. Check the data contract first (cheap, CPU, no torch needed) — this is the only stage that works in a CPU-only container:

python3 scripts/sidish_cell_death.py --stage dry-run \
    --input-h5ad data/adata.h5ad --bulk-csv data/bulk.csv \
    --survival-csv data/sample_info.csv --ppi-dir data/PPI \
    --output-dir out

Exit 0 = compatible; exit 2 = contract violation (it reports gene/sample overlap, event-value validity, PPI completeness). Writes out/Stat/dry_run_report.json.

2. Prepare (notebook Step 0) — SIDISH.SIDISH.preprocess(processed=True) intersects sc/bulk genes, merges survival (duration, event), writes the two files the trainer reads:

python3 scripts/sidish_cell_death.py --stage prepare --input-h5ad data/adata.h5ad \
    --bulk-csv data/bulk.csv --survival-csv data/sample_info.csv --output-dir out
# -> out/Rdata/processed_adata.h5ad (4102 x 1208)
# -> out/Stat/processed_bulk.csv (572 x [duration,event,1208 genes])

Use --no-processed to make SIDISH do its own QC/normalisation/HVG selection, and --batch-correction + --batch-key <obs col> for Harmony (sc) / ComBat (bulk).

3. Train (notebook Step 1) — VAE (Phase 1) + Deep Cox (Phase 2), iterations outer rounds, checkpoints into --model-dir:

python3 scripts/sidish_cell_death.py --stage train --device cuda \
    --processed-h5ad out/Rdata/processed_adata.h5ad \
    --processed-bulk-csv out/Stat/processed_bulk.csv \
    --output-dir out --model-dir out/Rdata/model
# library API used: sidish(adata, bulk, "cuda", seed=0)
#   .init_Phase1(225, 20, 32, [512,128], 256, "Adam", 1e-4, 1e-4, 0)
#   .init_Phase2(500, 128, 1e-4, 0, 0.2, 256)
#   .train(5, 0.95, 30, model_dir, 0, True, "fitted")
# -> out/Rdata/model/{vae_transfer, deepCox, W_matrix_{i}.csv, adata_SIDISH.h5ad}

4. Reload + visualise + markers + KM (notebook Step 2):

python3 scripts/sidish_cell_death.py --stage infer --device cuda \
    --processed-h5ad out/Rdata/model/adata_SIDISH.h5ad \
    --processed-bulk-csv out/Stat/processed_bulk.csv \
    --output-dir out --model-dir out/Rdata/model --distribution-fit fitted
# reload(model_dir); get_percentille(0.95)
# plot_CellType_UMAP(resolution=0.8, celltype=False); plot_HighRisk_UMAP(resolution=0.08, celltype=False)
# get_MarkerGenes(logfc_threshold=1); plot_KM(penalizer=10)

The UMAP calls write adata_SIDISH_embedding.h5ad with obsm['latent'], obsm['X_umap'], obs['leiden'], obs['SIDISH_'], obs['SIDISH__'].

5. In-silico perturbation (notebook Step 3) — needs the PPI network; --ppi-dir is staged via a symlink so the library's '../PPI/...' lookups resolve (see Pitfalls):

python3 scripts/sidish_cell_death.py --stage perturb --device cuda \
    --processed-bulk-csv out/Stat/processed_bulk.csv --ppi-dir data/PPI \
    --output-dir out --model-dir out/Rdata/model \
    --perturb-genes AGRN CEBPB RHOB CDK2
# get_percentille(0.95); run_Perturbation(n_jobs=4)
# plot_top_perturbed_genes(percent_change/delta_change, top_n=30)
# run_double_Perturbation(top 20 genes, threshold=0.8)

6. --stage all chains prepare → train → infer → perturb in one process (same SIDISH object, so p_flip/percent_change are available to the double-perturbation step). --smoke shrinks every numeric setting (epochs=1, iterations=1, small layers, ≤64-cell batches) for a best-effort CPU sanity run — never verified to converge (see Pitfalls).

Parameters

Notebook values are what script/SIDISH.ipynb actually passes; the script's defaults equal them.

Outputs

--output-dir is a root, not a flat folder: the script creates three sub-directories on start-up (idempotently) and writes every artifact into the one that matches its type — <out>/Img/ figures, <out>/Stat/ tables & text, <out>/Rdata/ objects & binaries. A run ends by printing the contents of all three.

--model-dir (default <out>/Rdata/model) and --processed-h5ad / --processed-bulk-csv (defaults <out>/Rdata/processed_adata.h5ad, <out>/Stat/processed_bulk.csv) are the only knobs that move files out of the layout; an explicit value always wins, and a --model-dir outside the root keeps its Rdata/-style checkpoint files wherever you point it.

Reading the results: risk_value is the Deep Cox output per cell (higher = worse predicted survival); SIDISH == 'h' ⇔ risk_value >= percentile_cells where percentile_cells is the --percentile quantile of the fitted hazard distribution (default 95th → ~5 % High-Risk cells; the local tutorial model gives 228/4102 = 5.6 % h). In perturbation tables, positive percent_change = knocking that gene pushes High-Risk cells back to Background (a candidate vulnerability); negative = the knockout increases risk. p_flip is a one-sided sign test on the flips, p_score tests the continuous shift.

Pitfalls

• GPU is effectively mandatory. Upstream README (local copy) states 环境：（需要GPU）. Phase 1 runs 225 epochs and the full loop repeats the VAE + 500-epoch Deep Cox model --iterations times, plus one forward pass per cell per gene in the perturbation step. Treat "hours on a GPU box" as the realistic budget; there is no published CPU timing.
• This container cannot run SIDISH. No torch, no SIDISH package, no GPU here (verified: import torch → ModuleNotFoundError, import SIDISH → ModuleNotFoundError, /opt/conda/bin/python3 3.12.14). Only --stage dry-run is runnable locally; every other stage exits 3 with an install hint. Do not promise training results from this container. --smoke exists for a GPU-less sanity run but its convergence is UNVERIFIED.
• --model-dir must end with /. The library builds checkpoint paths by string formatting ("{}vae_transfer".format(self.path) with no os.path.join); the script appends the separator for you, but hand-written calls must do the same.
• SIDISH reads '../PPI/*.txt' relative to the process CWD. ppi_network_handler.py does os.path.dirname('../data/') → '..', so the network is looked up in the parent of the working directory. The notebook comment ("脚本和data一定要放在平行文件夹") refers to this. The script copies the pattern by symlinking <output-dir>/.sidish_work/PPI -> --ppi-dir and chdir-ing into .sidish_work/run/ for perturbation calls.
• PPI files are not in the repo. The GitHub tree (checked) has no PPI/ directory and PyPI only ships the package. 9606.protein.links.v11.5.txt must come from STRING (https://stringdb-downloads.org/download/protein.links.v11.5/9606.protein.links.v11.5.txt.gz, verified 200) and the info file from the matching protein.info.v11.5 URL (verified 200). HIPPIE: http://cbdm-01.zdv.uni-mainz.de/~mschaefer/hippie/hippie_current.txt (301 → 200, 138 MB).
• The local hippie_current.txt is a broken download. It is 9,557 bytes and starts with <!DOCTYPE html … — an HTML error page, not the 138 MB HIPPIE table. load_network() splits on tabs and silently skips lines with <5 fields, so nothing crashes: you just lose the HIPPIE contribution and keep STRING. Verify the file's first bytes before trusting a perturbation run.
• load_network() loads the whole 549 MB STRING file into RAM (f.readlines()[1:]) and builds a Python dict of interactions; and plot_perturbation_UMAP_default/differential call load_network() once per gene inside the loop (line ~982 / ~1069). On a 4-gene panel that is 4× the I/O. Budget tens of GB of RAM for the STRING step, and prefer running it once per process.
• Gene symbols are required. load_network() maps STRING proteins to gene symbols and keeps only pairs whose symbol is in adata.var.index. An .h5ad indexed by Ensembl IDs gives an empty network → every "knockout" silently degrades to a plain single-gene zeroing, and the perturbation scores become meaningless. Check var_names[:5] first.
• get_percentille() always fits a Weibull, regardless of --distribution-fit: Utils.get_threshold calls weibull_min.fit/ppf directly, while train() may have selected Weibull/Gamma/Exponential by AIC/BIC/KS/R². So the reload-side threshold can differ in family from the training-side labelling — keep --distribution-fit consistent between train and infer, and expect the h/b counts to shift slightly.
• get_MarkerGenes() mutates the AnnData in place (normalize_total + log1p on self.adata). Call it after the embeddings are saved (the script does), and don't reuse the object expecting raw counts — in --stage all the script deliberately re-reads adata_SIDISH.h5ad before the perturbation step so the hazard percentile is computed on the un-normalised matrix.
• run_Perturbation() re-reads <model-dir>/adata_SIDISH.h5ad and overwrites self.adata, so it ignores anything you did in memory; the script re-reads the file and recomputes get_percentille() before it (as the notebook does). plot_perturbation_UMAP_* and run_double_Perturbation also re-read that file internally.
• --n-jobs is currently ignored: run_Perturbation(n_jobs=…) then hardcodes run_parallel_processing(self.adata, n_jobs=4) upstream; the PPI threshold inside that call is likewise hardcoded to 0.7 while the UMAP panels use 0.8. Reported here because it silently changes what you think you configured.
• plot_double_Perturbation_Heatmap() is broken upstream — it evaluates self.percent_change[gene][id] (builtin id) → TypeError. The reference notebook replaces it with plain matplotlib, and the script reproduces that fallback (Img/double_perturbation_heatmap.*). Note the double-perturbation keys are strings "GENE1+GENE2".
• plot_CellType_UMAP/plot_HighRisk_UMAP raise ValueError if celltype=True and a resolution is passed, and each call recomputes latent → neighbours → UMAP → Leiden (two calls = twice the work). plot_CellType_UMAP colours by the leiden column, so it only works after clustering; --skip-umap avoids the cost.
• Notebook inconsistency: init_Phase1 uses batch_size=256 in the training cell but 512 in the reload cells. Only architecture-defining arguments (--latent-size, --layer-dims) must match across train/reload; a mismatch there makes load_state_dict fail.
• The tutorial .h5ad already carries SIDISH, SIDISH_value, risk_value columns (168 h / 3934 b in data/adata.h5ad) — leftovers from an earlier run. Training overwrites them; don't read them as ground truth. --label-key only selects which obs column the script exports/compares.
• Event status is validated strictly: map_event_column accepts only dead/alive/1/0 and raises ValueError otherwise (e.g. censored, 2, NaN). The --stage dry-run check flags this before you allocate a GPU.
• bulk.csv orientation: genes are rows and samples are columns (55,402 × 572 here); the notebook's .T is what makes preprocess()'s "columns are genes" assumption true. A transposed input file silently yields 0 overlapping genes.
• Memory model: adata.X is a dense float32 array and init_Phase1 allocates W_matrix = np.ones(adata.X.shape). Large cohorts (10⁵ cells × 2 000 genes) mean ~800 MB per dense copy — subset cells rather than hoping.

Environment

This container (verified):

Upstream installation (per README, verified URLs):

conda create --name sidish_env python=3.12 && conda activate sidish_env
# install the PyTorch build matching your GPU first (https://pytorch.org/get-started/locally/)
git clone https://github.com/mcgilldinglab/SIDISH.git && cd SIDISH
pip install -e .            # or: pip install .      /   pip install SIDISH==1.0.0
pip install -r requirements.txt

requirements.txt (upstream) pins: anndata 0.11.1, scanpy 1.10.4, numpy 1.26.4, pandas 2.2.3, scipy 1.14.1, scikit-learn 1.5.2, matplotlib 3.9.2, seaborn 0.13.2, shap 0.46.0, statsmodels 0.14.4, lifelines 0.30.0, pyro-ppl 1.9.1, leidenalg 0.10.2, python-igraph 0.11.8, imbalanced-learn 0.12.4, h5py 3.12.1, numba 0.60.0, torchvision 0.21.0, tqdm 4.67.1, umap-learn 0.5.7, bioinfokit 2.1.3.

Environment the reference notebook was produced in (a dedicated conda env, per the project README, which states 环境：（需要GPU）): python 3.12 with torch 2.6.0, scanpy 1.10.4, anndata 0.11.1, lifelines 0.30.0, shap 0.46.0, pyro-ppl 1.9.1, torch-geometric 2.7.0, matplotlib 3.9.2 (versions read from that environment's site-packages). Run the script with that interpreter (or an equivalent GPU env) for the training stages.

UNVERIFIED: no end-to-end SIDISH training was executed while writing this skill, so absolute runtimes, memory peaks and whether --smoke converges are unmeasured.