Benchmark rebuild, dataset reprocessing, TNBC cohort, and manuscript figures - #81
Open
MohamedOmar2020 wants to merge 279 commits into
Open
Benchmark rebuild, dataset reprocessing, TNBC cohort, and manuscript figures#81MohamedOmar2020 wants to merge 279 commits into
MohamedOmar2020 wants to merge 279 commits into
Conversation
Addresses reviewer concern that 50-gene/20%-signal benchmark is regime-specific. Tests whether dreamlet/NEBULA null-gene FPR inflation attenuates with larger panels and lower signal fractions. New sensitivity grid: Panel sizes: 50, 200, 500, 2000 genes Signal fractions: 1%, 5%, 10%, 20% + pure null per size = 20 scenarios × 200 iterations × 4 methods (two-arm only) Changes: - orchestrator.py: add build_sensitivity_grid() and run_sensitivity_benchmark() functions - run_benchmark.py: add --phase sensitivity entry point - slurm_sensitivity.sh: dedicated SLURM job (72h, 32 CPUs, 256GB) - supp_fig7_signal_fraction_sensitivity.py: new SF7 with 4 panels: A) FPR heatmap (method × panel size × signal fraction) B) FPR line plots per panel size C) Pure-null lambda_GC across panel sizes D) QQ contrast: 50g/20% vs 2000g/1% Smoke tested locally: all 4 methods produce valid results across 50-gene to 2000-gene panels. 2000-gene dreamlet ~65s/iter, NEBULA ~318s/iter; estimated HPC total ~4-8 hours with 30 workers.
- dreamlet: timeout 600s → 1800s for large-panel scenarios - NEBULA: timeout 1200s → 2400s for large-panel scenarios - New SLURM script slurm_sensitivity_2k_rerun.sh: re-runs only the 5 × 2000-gene scenarios with reduced parallelism (10 workers vs 30) to avoid memory pressure and subprocess timeouts
Replaces the seven original benchmark panels with publication-quality versions driven by the signal-fraction sensitivity dataset (11M rows, 4 panel sizes x 5 signal fractions x 4 methods). New panels: H Null-gene FPR curves faceted by panel size (50/200/500/2000) I Null-gene FPR heatmap (method x panel size x signal fraction) J Pure-null lambda_GC across panel sizes K QQ plots at a challenging regime (200 genes, 10% signal) L Pure-null FPR dot-and-whisker with Wilson 95% CIs M Effect-size estimation accuracy (bias + RMSE on signal genes) N Runtime scaling across panel sizes (log-log) Key findings surfaced by the new panels: - All four methods are nominally calibrated under pure null (lambda_GC ~= 1, FPR within 3-7% band) - dreamlet's null-gene FPR inflation scales with signal fraction and is independent of panel size (same inflation at 50 vs 2000 genes) - dreamlet also exhibits ~50% positive bias in effect-size estimation on signal genes, with 3-5x higher RMSE than sctrial, Wilcoxon, NEBULA - sctrial and Wilcoxon (Delta scores) are the only methods that maintain both nominal calibration and unbiased estimation across all tested regimes - sctrial and Wilcoxon are ~1000x faster than the R-based methods Implementation notes: - New helper functions: _compute_null_fpr_table, _compute_signal_power_table, _compute_signal_bias_rmse_table - Consistent styling via _method_style, _add_nominal_band, _style_axis - Data source updated to manuscript/benchmark/sensitivity/sensitivity_combined.csv - Removed the standalone supp_fig7_signal_fraction_sensitivity.py (content folded into SF4)
Panel H (FPR curves): Wilcoxon was hidden underneath sctrial at the nominal 5% line. Added tiny per-method x-offsets (-0.3 to +0.3 signal %-points) so overlapping methods are visible side-by-side. Panel I (FPR heatmap): y-axis signal-fraction labels were only drawn on the first subplot because sharey=True stripped them from the others. Switched to sharey=False so every subplot shows its own y-tick labels. The first subplot still carries the axis-label text. Panel L (pure-null calibration): the old horizontal dot-and-whisker with per-method y-offsets made all four panel sizes cluster visually on top of each other. Redesigned as a log-scale line plot: x = panel size, y = pure-null FPR, one line per method with 95% Wilson CIs. All four panel sizes are now distinct points along each line. Panel M (effect-size accuracy): the bias axis used a symmetric y-range that wasted the bottom half (no method has large negative bias), making bars look "cut in half". Switched to asymmetric [min, max] limits. The in-axis legend overlapped with dreamlet's tall bars; moved it to a figure-level legend above the subplots. sctrial/Wilcoxon/NEBULA bars remain near zero because those methods ARE essentially unbiased — that is the correct, honest visual encoding.
Panels H, J, L, and N previously plotted panel size (50/200/500/2000) or signal fraction (1/5/10/20%) at their literal numeric positions on a linear or log axis, which produced uneven visual spacing between the four discrete levels. Switched all four panels to categorical x-positions (0, 1, 2, 3) with the raw values used only as tick labels, so every subplot now has four evenly-spaced x-ticks regardless of value span. Also set explicit MultipleLocator on y-axes that relied on matplotlib auto-ticking: H: y step 0.1 on [0.0, 0.7] J: y step 0.05 on [0.90, 1.15] L: y step 0.01 on [0.025, 0.085] M (bias row): y step 0.05 M (RMSE row): y step 0.05 N: y log scale (decade ticks) unchanged Panel I (heatmap) and Panel K (QQ with shared axes) already have uniform tick grids so they are unchanged.
added 9 commits
July 28, 2026 22:26
…chmark.py Compartmentalization + no-redundancy: the benchmark constants, manifest-addressed loaders, aggregation helpers, and shared renderers (mixed-FPR, pure-null Type I, BH-FDR, marginal power, runtime, QQ, beta-envelope, lambda_GC, bias/RMSE, family FPR, broken-axis + legend helpers) now live in ONE shared module. figure3.py and supp_fig5 both import from .._benchmark; no figure imports another figure's panel code. Plotting logic unchanged (verbatim move).
…valuability/convergence New per-replicate aggregation helpers (_per_scenario_tpr end-to-end/tested; _per_scenario_quality evaluability/convergence) and two renderers in the shared toolkit for the benchmark supplement: end-to-end BH TPR vs signal fraction, and gene evaluability / convergence-among-attempted-fits per method.
…e by robustness family (core grid) Discovery sensitivity now uses the shared broken-axis renderer (calibrated in the main region, NEBULA on the strip). Evaluability/convergence computed over the core robustness families so dreamlet's filtering appears exactly where it occurs (low cell-yield / empirical heterogeneous yield ~0.85-0.90) while convergence stays 1.0.
…-FPR architecture param
…ss (A-M) New self-contained supp figure drawing every panel from the shared _benchmark toolkit (no cross-figure imports): full mixed-signal FPR (balanced+one-directional), complete pure-null calibration, QQ, beta-envelope, NEBULA hierarchy validation (from the provenance-stamped diagnostic), full marginal curves, FDR-controlled discovery sensitivity, bias/RMSE, family null-FPR and end-to-end TPR, evaluability, convergence, end-to-end-vs-tested. Also moved _panel_bench_pure_null_fpr into the shared toolkit.
…ves become panel H SF5 is now A-G real-data sensitivity + H empirical power (four-row composite). The mixed-signal FPR and pure-null Type I benchmark panels moved to Supp Fig 8, removing the cross-figure redundancy.
Shared _benchmark.py: - Canonical legend/placement order (sctrial, Wilcoxon, limma-voom, dreamlet, NEBULA) for the QQ grid, envelope heatmap, bias/RMSE, family-FPR and power panels - Neutral-gray nominal reference line (was red, clashed with dreamlet's series); drop the shaded band so every calibration panel matches - Split-scale sequential colormap for the null-p envelope heatmap (expand 0-12% and 90-100%, reference at 5%) - architectures= loader param so the one-directional arm can feed the SF8 panel; keep bias/RMSE on the balanced architecture - Adaptive bias tick step; runtime speed-ratio callout now includes NEBULA - Pure-null Type I y-range/ticks harmonized with the vs-participants panel Figure 3: - Panel A leader lines terminate at label edges (white bbox), labels solid black; dropped in the composite cell only (standalone keeps all of them) - Panel B leave-one-out legend moved into a clear margin, ordered top-to-bottom - Forest title uses a colon; centered two-line composite titles for D and E - Remove a stray ax.legend() that reverted panel 3G to draw order Supp Fig 5: - Remove the duplicated benchmark renderers (now shared in _benchmark.py) and their now-unused imports/constants - Analytical/bootstrap key into the empty right margin; stronger, boxed scatter labels - Docstring/print cleanup to the A-H layout Supp Fig 8: - Full A-O benchmark composite; one-directional mixed-FPR panel populated - NEBULA hierarchy panel from the provenance-stamped sigma_u ablation diagnostic scripts/verify_vaccine_forest.py: - Full-precision check confirming the two coincident Vaccine forest effect sizes are a real standardized-effect coincidence, not a duplicated row
MohamedOmar2020
changed the title
![chatgpt-codex-connector[bot]](https://avatars.githubusercontent.com/in/1144995?s=60&v=4){kind=small}
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: aac89051f7
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| panels = nested_panels(probe, rng=np.random.default_rng(seed + 1)) | ||
| panel = [f"gene_{i}" for i in panels[50]] | ||
| effects = ( | ||
| make_signal(panel, frac, "balanced", 0.5, rng=np.random.default_rng(seed + 2)) |
There was a problem hiding this comment.
make_signal accepts an integer count and immediately applies int(n_signal), so passing frac=0.20 here creates zero affected genes. Consequently, the phase labelled signal is another null simulation; the same mistake in gates.py:487-489 makes Gate D's default signal_fraction=0.2 produce empty signal arrays and NaN ablation summaries. Convert each fraction to a realizable count such as round(frac * len(panel)) before both calls.
Useful? React with 👍 / 👎.
Comment on lines
+713
to
+716
| if n_done >= n_iterations: | ||
| print(f" [{si + 1}/{len(scenarios)}] {name} — CACHED, skipping") | ||
| all_results.append(existing) | ||
| continue |
There was a problem hiding this comment.
If a 72-hour job is killed after the base batch is flushed or during an adaptive extension, its CSV already has at least n_iterations but deliberately has no completion record. This branch then skips it forever based only on row count, while aggregate_benchmark.py rejects that same shard as an orphan, so resubmitting cannot recover the run. Treat a shard as cached only when its completion record exists and agrees with the CSV; otherwise rerun or continue it.
Useful? React with 👍 / 👎.
Comment on lines
+20
to
+21
| ENV_YML="$REPO_ROOT/src/sctrial/benchmark/sctrial_bench_environment.yml" | ||
| ENV_NAME="sctrial_benchmark" # must match slurm scripts |
There was a problem hiding this comment.
A fresh setup creates sctrial_benchmark, while the environment YAML names sctrial_bench and every submitted SLURM/deployment command invokes micromamba run -n sctrial. Following the documented setup therefore leaves no environment named sctrial, causing jobs to fail immediately or, worse, use an unrelated stale pre-existing environment. Align the setup, YAML, and all job commands on one name.
Useful? React with 👍 / 👎.
Comment on lines
+44
to
+46
| # Drop cells with no library — a zero offset is undefined | ||
| keep <- colSums(counts) > 0 & meta$lib_size > 0 | ||
| if (sum(keep) < 2) stop("Too few cells with non-zero counts after filtering") |
There was a problem hiding this comment.
At this point counts contains only the tested panel, so colSums(counts) > 0 removes cells whose full-transcriptome library is valid but which happen to have no reads in that panel. This is common for the 50-gene low-expression panels, and the probability of removal can differ by treatment and visit when signal is present, making NEBULA's analyzed cell population outcome-dependent and different across panel sizes. Filter only on the supplied full-transcriptome meta$lib_size > 0 condition.
Useful? React with 👍 / 👎.
Comment on lines
+44
to
45
| 7: ("Supp Fig 7: Heterogeneity and Temporal Dynamics", "supp.supp_fig7_heterogeneity_temporal"), | ||
| } |
There was a problem hiding this comment.
The commit adds supp_fig8_benchmark_calibration.py, but the supplementary registry still ends at 7. As a result, both the advertised all/supplementary generation paths omit the new benchmark figure, and --supp-fig 8 reports it as unknown. Add Figure 8 to SUPP_FIGURES so the manuscript entry point actually generates it.
Useful? React with 👍 / 👎.
MohamedOmar2020
changed the title
added 4 commits
July 29, 2026 10:46
…e, SF8 registry - Ablation study and Gate D passed a signal FRACTION to make_signal, which takes an integer gene count (int(0.2) is 0), silently turning the "signal" runs into null runs. Convert the fraction to a realizable count on the panel at both call sites (scripts/run_benchmark.py, gates.composition_ablation). The production benchmark path was already correct (orchestrator passes scenario counts and verifies them). - Align the HPC conda env name on "sctrial" across setup_hpc_env.sh, the env YAML (name, and drop the stale machine-specific prefix), and the setup doc. Every SLURM job runs `micromamba run -n sctrial`, so the previous "sctrial_benchmark" and "sctrial_bench" names left a fresh setup with no environment the jobs could find. - Register Supplementary Figure 8 in manuscript_figures/run_all.py so --supp-fig 8 and the all/supplementary generation paths actually produce it.
…tract Add a cell_filter argument to nebula_runner.run with two options: - panel_and_lib (default): colSums(counts) > 0 & meta$lib_size > 0, the historical filter retained verbatim so the frozen v1.0.0 run stays byte-for-byte reproducible. - lib_only: meta$lib_size > 0 only. This is the correct filter: the offset is the full-transcriptome library size, and the panel-count term drops cells with no reads in the tested panel, which is outcome- and panel-size-dependent. The default is unchanged, so all existing results are identical. Add scripts/verify_nebula_cell_filter.py: a Phase-0 sensitivity check that runs both filters on the SAME seeded simulated data, through the canonical simulate -> contracts.prepare_inputs -> nebula path (reusing run_benchmark._load_frozen_config), and reports the null-FPR, mean|beta| and retained-cell deltas per scenario. It decides whether the frozen NEBULA numbers need regenerating with lib_only.
…source-hash guard The guard in _load_frozen_config asserts the source tree matches the frozen benchmark, which is correct for a frozen production run but wrong for this diagnostic, which deliberately varies the NEBULA cell filter. Read the same config fields (blob[config] minus scenario-owned) directly instead.
Reuses the manuscript_figures._benchmark loaders and aggregation helpers (the same code behind Figure 3 and Supp Fig 8) to compute every endpoint per method from the frozen combined CSVs: calibration (pure-null Type I, mixed-signal FPR balanced and one-directional, % outside 95% CI, realized FDR), power and discovery (marginal detection, end-to-end BH TPR), estimation (bias/RMSE), robustness families (FPR/TPR/evaluability/convergence), and runtime.
MohamedOmar2020
force-pushed
the
omar/qc-ambient-fixes
branch
from
July 29, 2026 21:38
89ca435 to
11d1a4b
Compare
added 13 commits
July 29, 2026 22:50
The permutation-validation panels ran a single-threaded 999-iteration loop that took ~4 h to rebuild cold on the largest cohort. Each permutation is independent, so fan them out across the node's cores with a fork Pool: workers inherit the large read-only AnnData via copy-on-write instead of re-pickling it, and every permutation reseeds independently (base_seed + i) so the loop is both parallel and reproducible. Falls back to a serial loop when one worker is requested or FIGURE_PERM_SERIAL=1 is set. Because parallel reseeding changes the exact null draws (statistically equivalent), the SF and TNBC permutation cache keys are bumped (figure6_sf_perm_v1 to v2, figure5_tnbc_perm_v2 to v3).
The parallel permutation render OOM-killed the node on the Sade-Feldman cohort (16k cells x 55,330 genes): forking the full transcriptome across 32 workers balloons memory through copy-on-write refcount churn, even though the matrix is never read. The permutation features are participant-level signature scores held in .obs, and did_table resolves obs-column features without indexing the gene matrix, so the fit only needs .obs. Slim the AnnData to a single token gene before forking; each worker copy drops from gigabytes to a few MB and the numeric result is bit-identical (same obs, same seeds). TNBC (20k genes) completed at 32 workers but is slimmed too. Recycle workers with maxtasksperchild to bound any residual growth, and stop mutating the caller's AnnData (the slim is a copy). Also correct three completion-print labels left from the figure renumbering (Figure 6 -> 5, SuppFig3 -> 4, SuppFig6 -> 7; output filenames were already right) and tidy pre-existing import lint in the two supp modules.
… SF7 legend Figure 4 panel B (top-genes waterfall): the significance shading was silently defeated. Per-bar colors encoded n.s. via an alpha suffix, but ax.barh was called with a scalar alpha=0.9 that overrides the per-color alpha, so every top-effect gene rendered as if p<0.05. This made the waterfall contradict the volcano (panel A), where the same large-beta, large-SE genes are correctly non-significant. Encode significance as per-bar RGBA (0.9 sig / 0.35 n.s.) and drop the scalar alpha so n.s. genes render faded. Figure 2 panel B (cell vs participant beta scatter): a signature marker sat in the lower-left cluster while its label was flung to the corner with no leader line, reading as a dropped point. Set explicit symmetric limits so no marker can clip, and place labels with adjustText + thin leader lines. Supp Fig 7 panels A/D (participant x feature heatmaps): the arm/response legend, anchored just above the axes, overprinted the centered title in the standalone panels. The row-annotation helper now reports its legend row count and the standalone title is lifted clear of it. Also apply pre-existing ruff whitespace/import fixes in the two touched main modules.
… forest legend SF6 panel A volcano: strengthen adjustText (force_static replacing the deprecated force_points, add force_explode, wider expand/max_move) so dense gene labels stop overprinting (CCL8 over DKFZP434A062). SF3 panel B (beta scatter): replace hand-tuned index-parity offsets with adjustText + leader lines and explicit symmetric limits, resolving the Memory/Regulatory T Cell label collision. SF3 panel I (per-participant heatmap): keep the arm-colour sidebar inside the axes (widened left xlim) so it no longer overprints the participant-ID y-tick labels, and drop the response legend below the rotated x-tick labels. SF3 composite Row 1: widen the C|D|E gutter so panel C's title and p=0.05 annotation stop bleeding into panel D. Fig2 panel D forest: anchor the arm legend below the x-axis title (ncol=2) instead of overprinting it.
The 2x3 interaction grid was packed with inner_hspace=0.45, so the top-row Pre/Post x-tick labels overprinted the bottom-row subplot titles in the combined artboard. Widen to 0.85.
SF1 panel G (post-QC threshold compliance): in the small composite cell the six datasets' rotated count labels overprinted each other and the top-anchored legend occluded the tallest bars. Add a compact mode (used by the composite) that omits the per-bar value labels and drops the legend below the axes; the standalone panel keeps both. SF2 panel A (cell-type UMAP row): the 2-column per-UMAP legends were wider than their narrow column and overprinted the neighbouring dataset's legend. Truncate long cell-type names, cap the entries, and use a single narrow column in the composite; the standalone panels are unchanged.
…djustText Figure 4: the TNBC panels A-D carried copy-pasted melanoma legends (Responder/Non-responder, 'Melanoma DiD') that a post-hoc string-replacement helper rewrote to the TNBC arm names. Emit the correct labels at the source using the existing arm-name constants and delete that helper, so there is one source of truth. Panel B (top genes) and the pathway panels are unchanged: they already select the top up/down genes by effect size and the top pathways by NES, nothing hardcoded. Replaced en/em-dash panel titles with colons and dropped an unused variable. SF2: panel A cell-type UMAP legends are two columns again, with shorter labels and tighter spacing so the two columns still fit inside each UMAP cell. SF1: panel C arm legend moved below the axes (the many-arm upper-right legend overlapped the tallest points once shrunk into the composite); panel B gains bar headroom so value labels and the arm legend clear the tallest bar. SF4: panel D gene labels now repel from points via force_static (the old force_points is silently ignored by adjustText 1.3.0, so labels never separated from the dense near-origin cluster), with explode, leader arrows, in-axes constraint, and a legibility bbox; colon in place of a dash in the subtitle.
SF1: the composite reimplements panel B (participants per arm x visit) inline rather than calling the standalone panel function, so the standalone headroom fix did not reach it. Raise the inline y-axis ceilings (Melanoma 17.5 to 23, TNBC 8.8 to 11) so the upper-right arm legend clears the tallest bar. SF4: set the baseline-comparability axes frame with a 15% margin and draw the identity line before labelling, so adjustText has a stable frame and empty margin to spread the signatures that share the (0,0) point (CD19/IL2, FOXP3/IL2 were overprinting); unify the composite/standalone branches, raise the explode force, and switch to set_aspect(adjustable='box') so the aspect no longer reshapes the limits after the labels are placed.
The composite has a legend-restyling pass (_inside) that re-created panel C's 12-entry arm legend inside the plot, overriding the standalone placement and overlapping the boxplots at every in-axes location. Anchor it below the axes in six columns, matching the standalone panel; the hspace=0.32 row gap leaves room.
The below-axes six-column arm legend ran into panel D's title. Give panel C its own dedicated block: 7 columns, smaller font, so the ~13-label union of every dataset's arms fits in two rows within the row gap without touching panel D.
Below-the-axes still grazed panel D's title. The boxplots top out near log10 ~4.6 while the axis reaches 5.7, so raise the ceiling to 6.3 and put the compact two-row arm legend in that top headroom, clear of both the data and panel D.
…filter The GSEA pathway panels (Fig4 C and N, the response-stratified R-vs-NR GSEA) and SF6 restricted pathways to an immune+metabolic keyword include-list before taking the top-N by NES, so they showed the top immune/metabolic pathways rather than the global top. Remove that filter at all four Fig4 sites and delete the now-unused keyword list, regex, and _is_immune_or_metabolic helper; drop the unused import in SF6 (it never called the filter). Pathway panels now select the global top pathways by NES. Gene panels were already unfiltered.
Rebuild and expand the curated GENE_SIGNATURES gene sets against MSigDB (Hallmark/Reactome/KEGG) and the dataset source papers: - Rebuild Apoptosis (Hallmark/p53), T Cell Exhaustion (terminal, TOX-anchored) and NK Cell Activity (effector; drop KIR dropouts and cross-lineage adaptors) - Split the old Immune Exhaustion into terminal Exhaustion plus acute Checkpoint Receptors - Split Interferon Response into Type I and Type II - Move MHC-II (HLA-DR) from T Cell Activation to Antigen Presentation - Add Humoral/Plasma, Monocyte/Macrophage, Tissue-Resident Memory, Tfh/Tph and Hypoxia modules - Drop EOMES from Memory (cross-dataset sign flip), IL10/TGFB1 from Treg, and FTL/FTH1 from Oxidative Stress (myeloid loading) Update display names, the clinical signature subset, and the pre-specified endpoint references in Figure 3 and SuppFig 5. Fix the SuppFig5 pickle cache key to include the code version so signature changes invalidate the cache instead of serving stale results. Move the SuppFig3 forest legend to clear space.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Merges the benchmark, dataset, and figure work into dev. It combines commits from several feature branches (contributors: Mohamed Omar, Priyanka Vasanthakumari, Itzel). Merges cleanly, since dev is fully contained in this branch.
Simulation benchmark
NEBULA hierarchy diagnostic
Dataset reprocessing (4-team dataset review)
TNBC / breast-cancer cohort
Manuscript figures
Testing, docs, infra
Verification
Notes