From 566b6f46ff440928fea1e6c5306f7b137d6b5d11 Mon Sep 17 00:00:00 2001 From: Maggie Shin-Young Chen Date: Wed, 30 Sep 2026 15:32:18 -0700 Subject: [PATCH] CRT: option to sample NTC null groups with replacement (as SCEPTRE does) Add replace=False to make_ntc_groups_matched_by_freq / make_ntc_groups_ensemble, exposed as --ntc_replace, plus --max_ntc_groups to cap groups per ensemble. Without replacement, NTC guides are consumed as groups are built, so the null runs out after about n_ntc / --number_guide groups, and frequency-bin matching shrinks that further (e.g. ~200 NTC guides with 15 guides per target gave only 3-6 groups per ensemble). With replace=True each group is drawn without replacement but guides return to the pool, and identical groups are rejected, mirroring SCEPTRE's negative-control group construction (sample_combinations_v2, Katsevich-Lab/sceptre). With no cap, one group is built per real target. Default behavior is unchanged: replace=False reproduces the original groups exactly (checked against the unmodified function over 20 configurations). Co-Authored-By: Claude Opus 5.5 (1M context) --- .../references/03-calibration.md | 2 + .../references/parameter-catalog.md | 4 +- CHANGELOG.md | 6 + .../B_Calibration/Slurm_version/CRT/CRT.py | 5 +- .../B_Calibration/Slurm_version/CRT/CRT.sh | 5 + .../Slurm_version/CRT/CRT_parallel.sh | 5 + .../B_Calibration/Slurm_version/CRT/README.md | 4 +- .../B_Calibration/src/CRT/ntc_groups.py | 53 +++++++- .../Stage2_Evaluation/test_crt_ntc_replace.py | 120 ++++++++++++++++++ 9 files changed, 194 insertions(+), 10 deletions(-) create mode 100644 tests/Script/Stage2_Evaluation/test_crt_ntc_replace.py diff --git a/.claude/skills/perturbNMF-runner/references/03-calibration.md b/.claude/skills/perturbNMF-runner/references/03-calibration.md index c4673df..6ac5175 100644 --- a/.claude/skills/perturbNMF-runner/references/03-calibration.md +++ b/.claude/skills/perturbNMF-runner/references/03-calibration.md @@ -50,6 +50,8 @@ Key parameters (always ask): | `--resampling` | `bernoulli` | Null treated-set draw: `bernoulli` (size varies) or `fixed_count` (treated count fixed within each categorical-covariate stratum; exact for categorical-only covariates, approximate Pareto sampling with continuous ones). `fixed_count` is calibrated for few-cell targets; output names don't change, so use `--save_dir` to keep runs apart | | `--outcome` | `clr` | CRT outcome: `clr` (CLR of usage floored at 1e-8; effect `log2FC`) or `usage` (per-cell usage share, no log/floor; effect `usage_share_diff`). Recommend `--resampling fixed_count --outcome usage`; output names don't change, so use `--save_dir` to keep runs apart | | `--matched_ntc_null` | off | Also write `{K}_CRT_matched_null__{condition}.txt`: one NTC pseudo-target per real target with exactly its cell count, scored with the same CRT (~2x runtime) | +| `--ntc_replace` | off | Build NTC null groups with replacement across groups (as SCEPTRE does); use when `--number_guide` is large relative to the NTC pool | +| `--max_ntc_groups` | None | Cap on NTC null groups per ensemble (with `--ntc_replace`, None = one per real target) | | `--save_dir` | auto | Custom output directory (default: `//Evaluation/_/`) | | `--skip_existing` | off | Skip the CRT recompute for a (K, sel_thresh, condition) when **both** its real and fake `.txt` exist, and regenerate the QQ `.png` from the cached raw p-values. Resume a preempted job **or** re-plot without recomputing. | diff --git a/.claude/skills/perturbNMF-runner/references/parameter-catalog.md b/.claude/skills/perturbNMF-runner/references/parameter-catalog.md index c798ba4..bc87933 100644 --- a/.claude/skills/perturbNMF-runner/references/parameter-catalog.md +++ b/.claude/skills/perturbNMF-runner/references/parameter-catalog.md @@ -445,10 +445,12 @@ CRT reads guide info (`obsm['guide_assignment']`, `uns['guide_names']`, `uns['gu | `--resampling` | str | `bernoulli` | `bernoulli` or `fixed_count` (choices enforced); `fixed_count` holds the treated count fixed within each categorical-covariate stratum in every permutation | | `--outcome` | str | `clr` | `clr` or `usage` (choices enforced); `usage` = per-cell usage share (row-normalized, no log/floor), effect column `usage_share_diff` instead of `log2FC` | | `--matched_ntc_null` | flag | off | Also score a cell-count-matched NTC null → `{K}_CRT_matched_null__{condition}.txt` | +| `--ntc_replace` | flag | off | Build NTC null groups with replacement across groups (never twice in one group; identical groups rejected), as SCEPTRE does | +| `--max_ntc_groups` | int | None | Cap on NTC null groups per ensemble (with `--ntc_replace`, None = one per real target) | | `--save_dir` | str | None | Custom save directory (default: `//Evaluation/_/`) | | `--skip_existing` | flag | off | Skip the CRT recompute for a (K, sel_thresh, condition) when **both** its real and fake `.txt` exist, and regenerate the QQ `.png` from the cached raw p-values (resume a preempted job or re-plot without recomputing) | -**Null calibration:** CRT builds NTC (non-targeting control) pseudo-gene groups of size `--number_guide`, frequency-matched to real genes, and compares them to real targets on a real-vs-NTC QQ plot (one `.png` per K/sel_thresh/condition). `--number_guide` therefore controls the null group size (no longer hardcoded to 6). Ensemble count / bin count / seeds are set inside the vendored `src/CRT/` package, not exposed as flags. +**Null calibration:** CRT builds NTC (non-targeting control) pseudo-gene groups of size `--number_guide`, frequency-matched to real genes, and compares them to real targets on a real-vs-NTC QQ plot (one `.png` per K/sel_thresh/condition). `--number_guide` therefore controls the null group size (no longer hardcoded to 6). Ensemble count / bin count / seeds are set inside the vendored `src/CRT/` package, not exposed as flags; `--ntc_replace` and `--max_ntc_groups` control how the groups are sampled. **Output files** (covariate token = `--covariates` then `log_`-prefixed `--log_covariates` joined by `_`, or `no_covariates` when none): - `{K}_CRT__{condition}.txt` — real results: `target_name, program_name, log2FC` (with `--outcome usage`: `usage_share_diff, control_usage_share, usage_share_relative_diff`), `p-value` (skew), `adj_pval`, `p-value_raw`, `adj_pval_raw`. diff --git a/CHANGELOG.md b/CHANGELOG.md index 7c92664..2a3c366 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,12 @@ Entries marked ⚠️ change pipeline output — re-run affected analyses. `--resampling fixed_count`. - CRT `--matched_ntc_null`: cell-count-matched NTC pseudo-targets as a calibration diagnostic for few-cell targets (`{K}_CRT_matched_null_*.txt`). +- CRT `--ntc_replace` (`replace=True` in `make_ntc_groups_matched_by_freq` / + `make_ntc_groups_ensemble`): NTC null groups are drawn with replacement across + groups (never twice in one group; identical groups rejected), following SCEPTRE's + negative-control group construction. Without it the NTC pool runs out after about + n_ntc / `--number_guide` groups. `--max_ntc_groups` caps groups per ensemble. + Default stays without replacement. - CRT NTC null now also carries the skew-normal p-value (`p-value`, `adj_pval` in `{K}_CRT_fake_*.txt`), plus a `_skew.png` QQ plot and a `_skew` NTC-significance summary on that scale. diff --git a/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/CRT.py b/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/CRT.py index 8f32391..3294ef6 100644 --- a/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/CRT.py +++ b/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/CRT.py @@ -140,7 +140,8 @@ def run_CRT(adata, k, sel_thresh, output_folder, args): n_ensemble=10, seed0=7, group_size=args.number_guide, - max_groups=None, + max_groups=args.max_ntc_groups, + replace=args.ntc_replace, ) # CRT p-values for each NTC group in each ensemble (null distribution): the @@ -565,6 +566,8 @@ def main(): parser.add_argument('--resampling', type=str, choices=list(RESAMPLING_METHODS), default='bernoulli', help='How CRT null treated sets are drawn. bernoulli (default): each cell independently ~ Bernoulli(propensity), so the null set size varies around the observed count. fixed_count: every null set keeps the observed treated count within each stratum of the categorical covariates (exact within-stratum permutation when all covariates are categorical; propensity-weighted Pareto sampling within strata, an approximation, when continuous covariates are present). fixed_count is calibrated for few-cell targets where bernoulli is anticonservative.') parser.add_argument('--outcome', type=str, choices=list(OUTCOMES), default='clr', help='CRT outcome Y per cell and program. clr (default): centered log-ratio of usage floored at 1e-8; effect reported as log2FC. usage: per-cell usage share (usage row-normalized to proportions, no log, no floor); effect reported as usage_share_diff, the covariate-adjusted difference in mean usage share (not a fold change), plus control_usage_share (covariate-adjusted baseline share) and usage_share_relative_diff (diff / baseline). usage with --resampling fixed_count is the recommended configuration (calibrated on a matched NTC null, no floor hyperparameter). Output file names do not change; use --save_dir to keep runs apart.') parser.add_argument('--matched_ntc_null', action='store_true', help='Also build a cell-count-matched NTC null (one NTC pseudo-target per real target with exactly its cell count, scored with the same CRT) and write {K}_CRT_matched_null_{covar_tag}_{condition}.txt. Roughly doubles CRT runtime. The frequency-matched guide groups rarely reach the few-cell regime, so this is the calibration diagnostic for small targets.') + parser.add_argument('--ntc_replace', action='store_true', help='Build the NTC null groups with replacement across groups (a guide can appear in several groups, never twice in one; identical groups are rejected), as SCEPTRE does for its negative-control gRNA groups. Default off: guides are consumed, which caps the null at about n_ntc_guides / --number_guide groups per ensemble (fewer after frequency-bin matching), too few when --number_guide is large. With --ntc_replace and no --max_ntc_groups, each ensemble builds one group per real target, so the NTC null costs about 10x (n_ensemble) the real-target CRT.') + parser.add_argument('--max_ntc_groups', type=int, default=None, help='Maximum NTC null groups per ensemble (default: None = as many as the NTC pool allows; with --ntc_replace, one per real target). Use to bound runtime with --ntc_replace.') parser.add_argument('--FDR_method', type=str, choices=['BH', 'StoreyQ'], default='BH', help='FDR correction method: BH (Benjamini-Hochberg) or StoreyQ (Storey Q-value) (default: BH)') parser.add_argument('--save_dir', type=str, default=None, help='Base directory under which {K}_{thresh} subdirs are created. Default: //Evaluation/') parser.add_argument('--skip_existing', help='If set, skip per-(K, sel_thresh, condition) computations whose output .txt files already exist. Useful for resuming preempted jobs.', action='store_true') diff --git a/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/CRT.sh b/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/CRT.sh index 9276b49..26f0f92 100644 --- a/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/CRT.sh +++ b/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/CRT.sh @@ -82,6 +82,11 @@ python3 /oak/stanford/groups/engreitz/Users/ymo/Tools/PerturbNMF/src/Stage2_Eval # Recommended: --resampling fixed_count --outcome usage # --matched_ntc_null off by default; append to also write a cell-count-matched # NTC null ({K}_CRT_matched_null_*.txt; ~2x runtime) + # --ntc_replace off by default; append to build NTC null groups with + # replacement across groups (as SCEPTRE does); use when + # --number_guide is large and the NTC pool is small + # --max_ntc_groups default None; cap on NTC null groups per ensemble + # (with --ntc_replace, None = one per real target) # --save_dir default //Evaluation/{K}_{thresh}/ # --covariates default None (e.g. --covariates pct_counts_mt doublet_scores) # --log_covariates default None (e.g. --log_covariates total_counts guides_per_cell) diff --git a/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/CRT_parallel.sh b/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/CRT_parallel.sh index 6a9b9b7..75afc04 100644 --- a/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/CRT_parallel.sh +++ b/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/CRT_parallel.sh @@ -113,6 +113,11 @@ python3 /oak/stanford/groups/engreitz/Users/ymo/Tools/PerturbNMF/src/Stage2_Eval # Recommended: --resampling fixed_count --outcome usage # --matched_ntc_null off by default; append to also write a cell-count-matched # NTC null ({K}_CRT_matched_null_*.txt; ~2x runtime) + # --ntc_replace off by default; append to build NTC null groups with + # replacement across groups (as SCEPTRE does); use when + # --number_guide is large and the NTC pool is small + # --max_ntc_groups default None; cap on NTC null groups per ensemble + # (with --ntc_replace, None = one per real target) # --save_dir default //Evaluation/{K}_{thresh}/ # --covariates default None (e.g. --covariates pct_counts_mt doublet_scores) # --log_covariates default None (e.g. --log_covariates total_counts guides_per_cell) diff --git a/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/README.md b/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/README.md index 1be0f07..88eb089 100644 --- a/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/README.md +++ b/src/Stage2_Evaluation/B_Calibration/Slurm_version/CRT/README.md @@ -133,7 +133,7 @@ The NTC null is meaningful only if each synthetic NTC "pseudo-gene" looks like a 3. **Real-gene "bin signature"** (`_real_gene_bin_signatures`). For each real gene, take its guides, sort them by prevalence, keep the first `--number_guide` of them, and record the multiset of their bin indices. This length-`--number_guide` vector of bins is the gene's *bin signature* — a compact fingerprint of "this gene is targeted by guides of these prevalence levels." -4. **Match NTC guides to a real gene's signature** (`make_ntc_groups_matched_by_freq`). Pool the NTC guides by bin. Repeatedly draw a real-gene bin signature at random and try to assemble a pseudo-gene by pulling NTC guides from exactly the bins that signature calls for (e.g. a signature `[2,2,5,7,7,9]` pulls two NTC guides from bin 2, one from bin 5, two from bin 7, one from bin 9). Guides are drawn without replacement within a replicate; if a bin is short of NTC guides the draw is skipped and another signature is tried. Each successful group is labelled `ntc_{idx}` and, by construction, carries the **same size and the same per-guide prevalence composition as some real gene** — so its union of tagged cells, propensity model, and resulting CRT p-value are directly comparable to a real gene's. +4. **Match NTC guides to a real gene's signature** (`make_ntc_groups_matched_by_freq`). Pool the NTC guides by bin. Repeatedly draw a real-gene bin signature at random and try to assemble a pseudo-gene by pulling NTC guides from exactly the bins that signature calls for (e.g. a signature `[2,2,5,7,7,9]` pulls two NTC guides from bin 2, one from bin 5, two from bin 7, one from bin 9). Guides are drawn without replacement within a replicate (default), or with `--ntc_replace` returned to the pool after each group so a guide can appear in several groups but never twice in one, with identical groups rejected. The latter follows SCEPTRE's negative-control group construction (`sample_combinations_v2`, [Katsevich-Lab/sceptre](https://github.com/Katsevich-Lab/sceptre)) and keeps the null from running out of groups when guides per target are many relative to NTC guides. If a bin is short of NTC guides the draw is skipped and another signature is tried. Each successful group is labelled `ntc_{idx}` and, by construction, carries the **same size and the same per-guide prevalence composition as some real gene** — so its union of tagged cells, propensity model, and resulting CRT p-value are directly comparable to a real gene's. 5. **Ensemble over seeds** (`make_ntc_groups_ensemble`). Repeat the matched partitioning `n_ensemble` times with different seeds so the null is averaged over many random NTC→real-gene assignments rather than one arbitrary partition. @@ -178,6 +178,8 @@ The p-values from step 5 are exactly the **NTC null p-values** written to `{K}_C | resampling | str | "bernoulli" | How null treated sets are drawn: "bernoulli" (each cell ~ Bernoulli(propensity); set size varies) or "fixed_count" (treated count fixed within each categorical-covariate stratum; see *Fixed-count resampling*). Output file names do not change with this flag — use `--save_dir` to keep runs apart | | outcome | str | "clr" | CRT outcome $Y$: "clr" (CLR of usage floored at 1e-8; effect column `log2FC`) or "usage" (per-cell usage share, row-normalized, no log/floor; effect column `usage_share_diff`, a difference in mean usage share). See *Outcome*. Output file names do not change — use `--save_dir` to keep runs apart | | matched_ntc_null | flag | off | Also build a cell-count-matched NTC null: one NTC pseudo-target per real target with exactly its cell count (whole NTC guides in random order, topped up from the next guide), scored with the same CRT. Writes `{K}_CRT_matched_null_{covar_tag}_{condition}.txt` and prints the share of null p < 0.05 by cells per target. Roughly doubles runtime. Not recomputed under `--skip_existing` | +| ntc_replace | flag | off | Build NTC null groups with replacement across groups: a guide can appear in several groups (never twice in one), and identical groups are rejected. This is how SCEPTRE builds its negative-control gRNA groups. Without it, guides are consumed, so each ensemble gets at most ~n_ntc / `--number_guide` groups (fewer after frequency-bin matching), which is too few when `--number_guide` is large. With no `--max_ntc_groups`, builds one group per real target per ensemble (~10x the real-target CRT runtime) | +| max_ntc_groups | int | None | Cap on NTC null groups per ensemble. None = as many as the pool allows (with `--ntc_replace`: one per real target) | | save_dir | str | None | Directory to save results and figures. If not provided, defaults to `//Evaluation/_/` | | skip_existing | flag | off | If set, skip the CRT recompute for any (K, sel_thresh, condition) whose **both** result files (real `.txt` and fake `.txt`) already exist, and instead **regenerate the QQ `.png` from the cached raw p-values**. Use to resume a preempted job or to re-plot without recomputing. | diff --git a/src/Stage2_Evaluation/B_Calibration/src/CRT/ntc_groups.py b/src/Stage2_Evaluation/B_Calibration/src/CRT/ntc_groups.py index 9b8239e..2793ce1 100644 --- a/src/Stage2_Evaluation/B_Calibration/src/CRT/ntc_groups.py +++ b/src/Stage2_Evaluation/B_Calibration/src/CRT/ntc_groups.py @@ -117,9 +117,25 @@ def make_ntc_groups_matched_by_freq( max_groups: Optional[int] = None, drop_remainder: bool = True, max_attempts: int = 100, + replace: bool = False, ) -> Dict[str, List[str]]: """ - Returns dict group_id -> list of NTC guide names (no overlap within replicate). + Returns dict group_id -> list of NTC guide names. + + replace=False (default): guides are consumed as groups are built, so no guide is + shared between groups within a replicate. + + replace=True: guides are returned to the pool after each group, so a guide can + appear in several groups (never twice in one group), and identical groups are + rejected. This mirrors how SCEPTRE builds its negative-control ("undercover") + gRNA groups: each group is an independent without-replacement draw of + calibration_group_size NTC gRNAs, rejection-sampled so that no two groups are + the same set (sample_combinations_v2 in src/negative_control_functions.cpp, + https://github.com/Katsevich-Lab/sceptre). Without replacement, the NTC pool is + used up after about n_ntc / group_size groups, and the per-bin matching shrinks + that further: with ~200 NTC guides and 15 guides per target, only 3-6 groups per + replicate could be formed, too few for a stable null. If max_groups is None, + replace=True builds one group per real-gene signature (i.e. per real target). """ rng = np.random.default_rng(seed) @@ -130,10 +146,16 @@ def make_ntc_groups_matched_by_freq( continue bin_to_guides.setdefault(bin_id, []).append(guide) - for guides in bin_to_guides.values(): - rng.shuffle(guides) + if not replace: + for guides in bin_to_guides.values(): + rng.shuffle(guides) + + # With replacement the pool never runs out, so the loop needs an explicit cap + if replace and max_groups is None: + max_groups = len(real_gene_bin_sigs) groups: Dict[str, List[str]] = {} + seen_groups: Set[Tuple[str, ...]] = set() attempts = 0 group_idx = 0 @@ -168,10 +190,24 @@ def make_ntc_groups_matched_by_freq( continue selected: List[str] = [] - for bin_id, need in counts.items(): - pool = bin_to_guides[bin_id] - selected.extend(pool[:need]) - del pool[:need] + if replace: + for bin_id, need in counts.items(): + pool = bin_to_guides[bin_id] + selected.extend(rng.choice(pool, size=need, replace=False).tolist()) + + # Reject a group identical to one already drawn (as SCEPTRE does) + key = tuple(sorted(selected)) + if key in seen_groups: + attempts += 1 + if attempts >= max_attempts and drop_remainder: + break + continue + seen_groups.add(key) + else: + for bin_id, need in counts.items(): + pool = bin_to_guides[bin_id] + selected.extend(pool[:need]) + del pool[:need] groups[f"ntc_{group_idx}"] = selected group_idx += 1 @@ -190,9 +226,11 @@ def make_ntc_groups_ensemble( group_size: int = 6, max_groups: Optional[int] = None, drop_remainder: bool = True, + replace: bool = False, ) -> List[Dict[str, List[str]]]: """ Returns list of group dicts, one per ensemble replicate. + See make_ntc_groups_matched_by_freq for replace. """ groups_ens: List[Dict[str, List[str]]] = [] for e in range(n_ensemble): @@ -205,6 +243,7 @@ def make_ntc_groups_ensemble( seed=seed0 + e, max_groups=max_groups, drop_remainder=drop_remainder, + replace=replace, ) groups_ens.append(groups) return groups_ens diff --git a/tests/Script/Stage2_Evaluation/test_crt_ntc_replace.py b/tests/Script/Stage2_Evaluation/test_crt_ntc_replace.py new file mode 100644 index 0000000..72aa18c --- /dev/null +++ b/tests/Script/Stage2_Evaluation/test_crt_ntc_replace.py @@ -0,0 +1,120 @@ +""" +Unit tests for NTC null-group sampling with replacement (replace=True in +make_ntc_groups_matched_by_freq / make_ntc_groups_ensemble). Synthetic inputs only. + +Usage: + eval "$(conda shell.bash hook)" && conda activate NMF_Benchmarking + cd + python -m pytest tests/Script/Stage2_Evaluation/test_crt_ntc_replace.py -v + +Behavior verified: + (a) replace=False (default) is unchanged: no guide shared between groups, and the + pool runs out after about n_ntc / group_size groups + (b) replace=True: no guide twice in one group, guides reused across groups, + no two identical groups, each group matches a real-gene bin signature + (c) replace=True caps at max_groups, or at one group per real-gene signature when + max_groups is None + (d) the ensemble wrapper passes replace through +""" + +import importlib.util +import os +import sys +from collections import Counter + +import pytest + +# Load the vendored CRT package under its own name (see test_crt_resampling.py) +CRT_PKG_DIR = os.path.abspath(os.path.join( + os.path.dirname(__file__), "..", "..", "..", "src", + "Stage2_Evaluation", "B_Calibration", "src", "CRT")) +if "crt_calibration" not in sys.modules: + spec = importlib.util.spec_from_file_location( + "crt_calibration", os.path.join(CRT_PKG_DIR, "__init__.py"), + submodule_search_locations=[CRT_PKG_DIR]) + crt_pkg = importlib.util.module_from_spec(spec) + sys.modules["crt_calibration"] = crt_pkg + spec.loader.exec_module(crt_pkg) + +from crt_calibration.ntc_groups import ( # noqa: E402 + make_ntc_groups_ensemble, + make_ntc_groups_matched_by_freq, +) + + +# 200 NTC guides in 4 frequency bins; 15 guides per real target (the regime where +# sampling without replacement runs out of groups) +N_NTC = 200 +N_BINS = 4 +GROUP_SIZE = 15 +N_REAL_TARGETS = 40 + + +@pytest.fixture +def ntc_inputs(): + ntc_guides = [f"ntc_g{i}" for i in range(N_NTC)] + guide_to_bin = {g: i % N_BINS for i, g in enumerate(ntc_guides)} + ntc_freq = {g: 0.01 for g in ntc_guides} + # Real-gene signatures: sorted bin ids, one per real target + real_sigs = [ + sorted((t + j) % N_BINS for j in range(GROUP_SIZE)) + for t in range(N_REAL_TARGETS) + ] + return ntc_guides, ntc_freq, real_sigs, guide_to_bin + + +def _make(ntc_inputs, **kwargs): + ntc_guides, ntc_freq, real_sigs, guide_to_bin = ntc_inputs + return make_ntc_groups_matched_by_freq( + ntc_guides=ntc_guides, ntc_freq=ntc_freq, real_gene_bin_sigs=real_sigs, + guide_to_bin=guide_to_bin, group_size=GROUP_SIZE, seed=7, **kwargs) + + +def test_default_consumes_guides(ntc_inputs): + groups = _make(ntc_inputs) + all_guides = [g for grp in groups.values() for g in grp] + assert len(all_guides) == len(set(all_guides)) + assert len(groups) <= N_NTC // GROUP_SIZE + + +def test_replace_no_repeat_within_group(ntc_inputs): + groups = _make(ntc_inputs, replace=True) + for grp in groups.values(): + assert len(grp) == GROUP_SIZE + assert len(set(grp)) == GROUP_SIZE + + +def test_replace_reuses_guides_across_groups(ntc_inputs): + groups = _make(ntc_inputs, replace=True) + counts = Counter(g for grp in groups.values() for g in grp) + assert max(counts.values()) > 1 + assert len(groups) > N_NTC // GROUP_SIZE + + +def test_replace_no_identical_groups(ntc_inputs): + groups = _make(ntc_inputs, replace=True) + keys = [tuple(sorted(grp)) for grp in groups.values()] + assert len(keys) == len(set(keys)) + + +def test_replace_groups_match_a_real_signature(ntc_inputs): + _, _, real_sigs, guide_to_bin = ntc_inputs + sigs = {tuple(s) for s in real_sigs} + groups = _make(ntc_inputs, replace=True) + for grp in groups.values(): + assert tuple(sorted(guide_to_bin[g] for g in grp)) in sigs + + +def test_replace_group_count_caps(ntc_inputs): + assert len(_make(ntc_inputs, replace=True)) == N_REAL_TARGETS + assert len(_make(ntc_inputs, replace=True, max_groups=9)) == 9 + + +def test_ensemble_passes_replace(ntc_inputs): + ntc_guides, ntc_freq, real_sigs, guide_to_bin = ntc_inputs + ens = make_ntc_groups_ensemble( + ntc_guides=ntc_guides, ntc_freq=ntc_freq, real_gene_bin_sigs=real_sigs, + guide_to_bin=guide_to_bin, n_ensemble=3, seed0=7, group_size=GROUP_SIZE, + max_groups=9, replace=True) + assert [len(groups) for groups in ens] == [9, 9, 9] + assert ens[0] != ens[1]