CRT: option to sample NTC null groups with replacement (as SCEPTRE does)

.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_<covariates>_{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: `<out_dir>/<run_name>/Evaluation/<K>_<sel_thresh>/`) |
 | `--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. |
 

.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_<covariate_token>_{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: `<out_dir>/<run_name>/Evaluation/<K>_<sel_thresh>/`) |
 | `--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_<covariate_token>_{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`.

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.

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: <out_dir>/<run_name>/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')

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 <out_dir>/<run_name>/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)

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 <out_dir>/<run_name>/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)

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 `<out_dir>/<run_name>/Evaluation/<K>_<sel_thresh>/` |
 | 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. |
 

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