B-ALL concordance analysis: in-silico perturbation vs. CRISPR essentiality

Goal

Test whether an in-silico model's predicted "impact" of perturbing ATAC-seq peaks is concordant with real CRISPR knockout essentiality, in B-cell acute lymphoblastic leukemia (B-ALL). There is no direct peak-to-gene assignment in the data, so the bridge between the two is transcription-factor identity: chromatin-factor binding annotations tell us which peaks each TF/factor is predicted to bind, and most of those factor names are themselves genes that were CRISPR-screened.

Repository layout

data/                              # raw inputs (see Data files below)
  BALL_CRISPRGeneEffect.tsv
  BALL_Expression.tsv
  BALL_Model.tsv
  factor_overlap_matrix.tsv
  impact_score_matrix.tsv
  subtype_meta.tsv
scripts/
  concordance_workflow.py          # main pipeline: builds factor-level scores, runs concordance test
  confounder_model.py              # multivariate model: essentiality ~ impact + confounders
  run_app.py                       # start FastAPI concordance UI (port 8765)
outputs/
  factor_concordance_results.tsv   # factor-level table from concordance_workflow.py
  concordance_scatter.png
  extended_factor_features.tsv     # factor-level table from confounder_model.py (more features)
  confounder_model_summary.txt     # full OLS / VIF / CV output, plain text
  partial_regression_plot.png
  feature_importance.png
  cache/                           # disk cache for on-demand bundles (created by run_app.py)
src/isgs/
  static/                          # web UI (index.html, app.css, app.js)
  ...                              # installable Python package (logic shared by scripts)
pyproject.toml

Quick start

Batch analysis (TSV + plots)

cd /path/to/ISGS
python3 -m venv .venv && source .venv/bin/activate
pip install -e .
python scripts/concordance_workflow.py
python scripts/confounder_model.py

Interactive app (FastAPI + disk cache)

The web UI runs as a local server that loads peak-filter bundles on demand and caches them to disk (~20–30s the first time per (impact, variance, subtype) tuple; instant on repeat).

source .venv/bin/activate
pip install -e .
python scripts/run_app.py
# open http://127.0.0.1:8765

Entry point: isgs-app (same server as scripts/run_app.py).

Frontend assets live in src/isgs/static/ and are served directly by the app.

Cache location: outputs/cache/ — analysis bundles, meta.json, peak summaries, and UI state (app_prefs.json, favorites.json). Cache entries are invalidated when any data/*.tsv file is newer than the cache file (UI prefs are not invalidated by data changes).

Both batch scripts accept --data-dir (default ./data) and --out-dir (default ./outputs). Entry points: isgs-concordance, isgs-confounders, isgs-app.

Data files

File	Shape	Rows	Columns	Notes
BALL_CRISPRGeneEffect.tsv	13 × 18,531	cell line (ModelID)	gene symbol	DepMap Chronos/CERES gene effect score. More negative = more essential. Has structural NaNs: 11 of 13 lines share one set of 744 unscreened genes, the other 2 lines share a different 372/373-gene set (different sgRNA library versions) — see "Gotchas" below.
BALL_Expression.tsv	26 × 19,215	cell line (ModelID)	gene symbol	Expression values (looks like log2(TPM+1)). 26 lines, a superset of the 13 CRISPR-screened lines. Has an extra IsDefaultEntryForModel column that scripts drop on load.
BALL_Model.tsv	33 × 48	cell line (ModelID)	metadata	Per-line metadata. PatientSubtypeFeatures has free-text fusion calls (e.g. KMT2A-AFF1, TCF3::PBX1) — this is the only way to approximate a subtype label for the CRISPR/expression lines, and its vocabulary does not match subtype_meta.tsv. Currently hand-curated into a dict (CELL_LINE_SUBTYPE in concordance_workflow.py) for the 13 CRISPR lines only.
factor_overlap_matrix.tsv	94,478 × 205	ATAC peak (peak_id)	chromatin factor name	Binary: 1 if that factor's binding site overlaps the peak. Row order is identical to impact_score_matrix.tsv (verified), but scripts still join on peak_id defensively rather than assume positional alignment.
impact_score_matrix.tsv	94,478 × 155	ATAC peak (peak_id)	chr, start, end, then 151 patient GEO IDs	In-silico perturbation impact score per peak per patient. Raw TSV values are small and strictly positive (~0.0004–0.15). The app multiplies by IMPACT_SCORE_SCALE (1000) at load time for display (~0.4–150). Treated as an unsigned magnitude of predicted importance. Worth confirming the sign/direction convention with whoever produced this model before trusting downstream interpretation.
subtype_meta.tsv	153 × 2	GEO ID	Subtype	Patient subtype labels for the impact-score columns. 151 of 153 patients here are present in impact_score_matrix.tsv. Subtype sizes are very uneven (e.g. Unknown=25, BCR-ABL1=20 down to iAMP21=1) — anything below ~10 patients shouldn't be trusted for subtype-level aggregation.

Gene/factor name reconciliation

Of the 205 factor names in factor_overlap_matrix.tsv, 202 match a gene symbol in the CRISPR data 1:1, and 201 also exist in the expression data. The three that don't: MED (a complex placeholder, not a single gene), PBX1-2-3 (a paralog group), and TCF3-PBX1 (the fusion product itself, not a wild-type gene). All three are dropped automatically (via .intersection() on column names) — nothing hardcoded, this just falls out of the join.

Scripts

concordance_workflow.py

1. Loads all six files, aligns factor_overlap_matrix and impact_score_matrix on peak_id.
2. Builds a Factor × Patient "predicted impact" matrix via one matrix multiply: factor_overlap.T @ impact_scores, normalized by peak count per factor. This is the computational core — avoids looping over 205 factors × 94,478 peaks.
3. Collapses to one global value per factor (mean across all patients) and per-subtype values (mean within subtype_meta.tsv groups).
4. Computes CRISPR essentiality per gene as -mean(gene_effect) across the 13 lines (so higher = more essential, matching the impact score's direction), and a mean-expression filter to drop factors not expressed in B-ALL (MIN_MEAN_EXPRESSION = 1.0).
5. Runs the primary concordance test: Spearman correlation + a 10,000-iteration label-permutation null (more robust than the parametric p-value at this n). Also reports median-aggregation and no-expression-filter robustness variants.
6. Runs a subtype-specificity check: for subtypes with ≥2 matched CRISPR lines (currently only KMT2A, n=4, and TCF3-PBX1, n=2, after excluding the Unknown/B-Other catch-all bins), compares matched vs. mismatched subtype pairings.
7. Saves factor_concordance_results.tsv and concordance_scatter.png.

Key tunables at the top of the file: MIN_PEAKS_PER_FACTOR (10), MIN_MEAN_EXPRESSION (1.0), N_PERMUTATIONS (10,000), CELL_LINE_SUBTYPE (the hand-curated mapping).

confounder_model.py

Takes the question further: is the impact-score/essentiality relationship confounded by other factor characteristics? Builds one extra engineered feature beyond what's in factor_concordance_results.tsv:

• mean_co_occupancy — average number of other factors co-bound at a factor's peaks (whether it mostly binds generically "busy" multi-TF hub regions vs. more selective sites).

(extended_factor_features.tsv also includes impact_sd_across_patients for reference, but it is not used in the OLS model.)

Fits an OLS regression (crispr_essentiality_global ~ predicted_impact_global + log10(n_bound_peaks) + mean_expression_BALL + mean_co_occupancy, all predictors z-scored), checks VIF for collinearity, compares nested models (impact-only vs. full), and cross-validates a linear model and a random forest (5-fold) to check out-of-sample predictive power. Saves extended_factor_features.tsv, confounder_model_summary.txt, partial_regression_plot.png, feature_importance.png.

Findings so far (as of last run)

• Global concordance is weak: Spearman rho = 0.092 (permutation p = 0.20, n=190 expressed factors). Holds up similarly under median aggregation / no expression filter.
• Restricting to the 20 most essential factors (essentiality > 1.0) strengthens it: Spearman rho = 0.347 (p=0.13), Pearson r = 0.415 (p=0.069) — but this is a non-random, outcome-selected subset, so it's descriptive, not a fair hypothesis test.
• The strongest single predictor of CRISPR essentiality in the multivariate model is mean_expression_BALL (β=0.221, p<0.001) — a generic, chromatin-unrelated effect. The impact score's adjusted coefficient (β=0.075, p=0.052) barely changes from its unadjusted value (β=0.079, p=0.031), so it isn't confounded by the other covariates tested, but it's weak and right at the significance boundary.
• 5-fold cross-validated R² is negative for both a linear model (-0.04) and a random forest (-0.08) — the in-sample fit (R²=0.24) does not generalize. This is the most important caveat on everything above: at n≈190 factors / 13 cell lines, don't treat the in-sample significance as confirmatory.
• Subtype-specificity check is inconclusive: KMT2A shows the expected matched > mismatched pattern, TCF3-PBX1 doesn't, but both groups have only 2-4 cell lines.
• Top concordant hits with good face validity (high impact AND high essentiality): MYC, CDK7, CDK9, CTCF, HCFC1, UBTF, MYB — all well-known pan-essential transcriptional/chromatin machinery genes. EBF1 (a core B-ALL lineage TF) also shows up in the high-essentiality subset.

Known data quirks / "gotchas" to remember

• CRISPR NaN structure: BALL_CRISPRGeneEffect.tsv is not fully dense. When subsetting to a small group of cell lines (e.g. for a subtype-specific essentiality mean), scipy.stats.spearmanr/pearsonr will silently return nan if even one paired value is missing — always drop NaNs pairwise first (see safe_spearman() in concordance_workflow.py) rather than calling the scipy functions directly on a small subset.
• subtype_meta.tsv vocabulary ≠ BALL_Model.tsv vocabulary. The former has clean canonical subtype names; the latter only has raw fusion-call free text. The current cell-line-to-subtype mapping is manually curated and only covers the 13 CRISPR lines — re-derive it if you add more lines or want a more principled mapping.
• Peak-to-gene linking is TF-mediated, not coordinate-based. impact_score_matrix.tsv has chr/start/end columns that are currently unused. If you get access to a gene annotation/TSS file, a nearest-gene assignment could be a complementary (and probably cleaner) linking strategy worth comparing against the current TF-binding-based approach.
• Sign convention of the impact score is assumed, not confirmed. All observed values are small and positive; if the underlying model is actually signed (e.g., direction of expression change rather than magnitude), the aggregation logic (mean of impact scores per factor) should be revisited.

Environment

pandas
numpy
scipy
matplotlib
statsmodels
scikit-learn
seaborn
plotly
fastapi
uvicorn[standard]

Optional for API tests: httpx (used by FastAPI's TestClient).

No GPU or special hardware needed; the heaviest step is the factor_overlap.T @ impact_scores matrix multiply (94,478 × 205 by 94,478 × 151), which runs in a few seconds with BLAS-backed numpy. impact_score_matrix.tsv is the large file (~270MB) — loading it is the main I/O cost (a few seconds with pandas.read_csv).

Suggested next steps

• Get clarity on the impact score's sign/direction convention from its source.
• If a gene-annotation file becomes available, add a coordinate-based peak→gene assignment as an alternative to the TF-binding bridge, and compare results.
• Expand/validate the CELL_LINE_SUBTYPE mapping (currently hand-curated from free text) — possibly cross-check against OncotreeSubtype or external annotations of these specific cell lines.
• Given the cross-validation results, consider whether a larger panel of CRISPR-screened B-ALL lines (if available from a broader DepMap pull) would give enough power for a more trustworthy multivariate fit.
• Consider restricting the concordance test a priori to known B-ALL lineage-defining TFs (PAX5, EBF1, IKZF1, TCF3, RUNX1, etc.) as a more targeted, biologically motivated hypothesis test rather than testing all 190+ factors.