Metadata-Version: 2.4
Name: pysceptre
Version: 0.2.0
Summary: Standalone Python port of sceptre's CRT discovery-analysis statistical engine (complement control group)
Author: Siddarth Wekhande, Eugenio Mattei, Broad Institute of MIT and Harvard
License-Expression: GPL-3.0-only
Project-URL: Homepage, https://github.com/broadinstitute/pysceptre
Project-URL: Repository, https://github.com/broadinstitute/pysceptre
Project-URL: Issues, https://github.com/broadinstitute/pysceptre/issues
Keywords: crispr,single-cell,sceptre,conditional-randomization-test,differential-expression,bioinformatics
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
License-File: THIRD_PARTY_LICENSES
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Requires-Dist: pandas>=2.0
Requires-Dist: threadpoolctl>=3.0
Provides-Extra: dev
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: ruff>=0.6; extra == "dev"
Requires-Dist: pre-commit>=3.0; extra == "dev"
Provides-Extra: fast
Requires-Dist: numba>=0.58; extra == "fast"
Provides-Extra: io
Requires-Dist: mudata>=0.3; extra == "io"
Requires-Dist: anndata>=0.10; extra == "io"
Requires-Dist: pyarrow>=14; extra == "io"
Requires-Dist: matplotlib>=3.7; extra == "io"
Provides-Extra: examples
Requires-Dist: pysceptre[io]; extra == "examples"
Requires-Dist: scanpy>=1.10; extra == "examples"
Requires-Dist: igraph>=0.11; extra == "examples"
Requires-Dist: leidenalg>=0.10; extra == "examples"
Provides-Extra: docs
Requires-Dist: mkdocs-material>=9.5; extra == "docs"
Requires-Dist: mkdocstrings[python]>=0.26; extra == "docs"
Dynamic: license-file

# pysceptre

A standalone Python port of the statistical engine behind
[`sceptre`](https://timothy-barry.github.io/sceptre-book/)'s discovery
analysis for single-cell CRISPR screens -- specifically the **complement
control group + CRT (conditional randomization test) resampling** path used
for high-MOI data.

It covers three of sceptre's analysis steps on that path -- the **discovery
analysis**, the **calibration check** and the **power check** -- and is *not*
a general
reimplementation. Targeting one validated path is what lets it batch the
linear-algebra work sceptre does per-gene/per-target in R/C++ loops into
vectorized numpy calls, cutting real-dataset runtimes from hours to tens of
minutes. See [Scope and limitations](#scope-and-limitations) for exactly what
is and isn't covered.

## Why this exists

<!-- --8<-- [start:why] -->

`sceptre`'s own discovery-analysis implementation is correct but processes
one (gene, gRNA-target) pair, and one resampling draw, largely with R and
per-cell C++ loops. At real dataset scale (hundreds of thousands of cells,
thousands of gRNA targets, tens of thousands of pairs) that adds up. Every
GLM fit in this pipeline shares one design matrix across many response
columns (all genes share the covariate matrix; all gRNA targets share it
too), which is exactly the shape numpy's batched matrix operations are
built for -- so the rewrite fits one batched IRLS call per resampling stage
instead of one per gene or target.

<!-- --8<-- [end:why] -->

## Installation

<!-- --8<-- [start:installation] -->

```bash
pip install -e ".[dev,fast]"
```

- `dev` installs `pytest` (needed for the test suite) plus `ruff` and
  `pre-commit` for linting and formatting.
- `fast` installs [`numba`](https://numba.readthedocs.io/), which JIT-compiles
  the CRT sampler's cell-grouping step (a counting sort). Without it,
  `pysceptre` falls back to a slower pure-numpy `argsort`-based version
  automatically -- everything still works, just slower at real dataset scale.
- `io` installs `mudata`, `anndata`, `pyarrow` and `matplotlib`. These are
  needed only by the dataset-export and benchmarking scripts under
  `scripts/`, **not** by the package: `run_discovery_analysis` takes in-memory
  arrays, so calling it does not require them.

Requires Python >= 3.10. Core dependencies: `numpy`, `scipy`, `pandas`,
`threadpoolctl`.

### Datasets

The analysis functions take plain arrays, so any route that produces them
works. The `scripts/` directory also carries one out of R end to end:
`export_sceptre_dataset.R` extracts a post-QC `sceptre_object`, and
`make_h5mu.py` converts that to a **MuData `.h5mu`** with two assays --
`rna` (cells x genes counts) and `grna` (cells x assignment units, whose
`var` marks each unit as a per-target union or an individual non-targeting
gRNA). `scripts/sceptre_io.load_export` reads one back, and with
`backed=True` serves genes from disk on demand rather than holding the whole
matrix.

<!-- --8<-- [end:installation] -->

## Quick start

<!-- --8<-- [start:quickstart] -->

```python
from pysceptre.pipeline.api import run_discovery_analysis

result = run_discovery_analysis(
    response_matrix=response_matrix,      # (n_genes, n_cells) dense array or scipy.sparse
    gene_ids=gene_ids,                    # row labels for response_matrix, in order
    covariate_matrix=covariate_matrix,    # (n_cells, p) numeric design matrix
    grna_target_cells=grna_target_cells,  # dict[target_id -> 0-based treated-cell indices]
    pairs=pairs,                          # DataFrame['response_id', 'grna_target']
    side="both",                          # the default, and the side validated against R
    seed=0,
)
# result: DataFrame with one row per pair --
#   response_id, grna_target, p_value, fold_change, se_fold_change,
#   pct_change_es, pct_change_es_ci_low, pct_change_es_ci_high, z_orig, stage
```

[`examples/scanpy_interop.py`](https://github.com/broadinstitute/pysceptre/blob/main/examples/scanpy_interop.py) runs a whole screen
inside an ordinary `scanpy` workflow -- QC and covariates out of `obs`, gRNA
assignments out of the gRNA modality's `var`, discovery and calibration, then
results back onto the same object and on to PCA and clustering, with no R at
any point. Run it with `uv run --extra examples python
examples/scanpy_interop.py`.

See [TUTORIAL.md](https://github.com/broadinstitute/pysceptre/blob/main/TUTORIAL.md) for a complete, runnable walkthrough
(including how to build each input from scratch) and for guidance on
picking `target_chunk_size` for your dataset.

<!-- --8<-- [end:quickstart] -->

## API reference

### `pysceptre.pipeline.api.run_discovery_analysis`

The one function most users need.

```python
run_discovery_analysis(
    response_matrix,
    gene_ids: list[str],
    covariate_matrix: np.ndarray,
    grna_target_cells: dict[str, np.ndarray],
    pairs: pd.DataFrame,
    *,
    side: str = "both",
    resampling_approximation: str = "skew_normal",
    seed: int | None = None,
    target_chunk_size: int = 200,
    chunk_memory_gb: float = 1.0,
) -> pd.DataFrame
```

| Parameter | Type | Description |
|---|---|---|
| `response_matrix` | `(n_genes, n_cells)` dense `ndarray`, `scipy.sparse` matrix, or a backed reader | Gene expression counts. Rows must correspond 1:1 with `gene_ids`, in order. Passing a genome-wide matrix is fine: only genes appearing in `pairs` are fitted, so untested rows cost storage but not compute. (Earlier versions fitted every row and this table told you to pre-filter; that is no longer necessary.) |
| `gene_ids` | `list[str]` | Row labels for `response_matrix`, in the same order as its rows. Matched against `pairs['response_id']`. |
| `covariate_matrix` | `(n_cells, p)` `ndarray` | Already formula-expanded numeric design matrix (intercept column, `log(umis)`, batch dummies, etc. -- whatever R's `model.matrix()` would have produced). `pysceptre` does not parse an R-style formula DSL; build this matrix yourself, or extract it directly from an existing `sceptre_object`'s `@covariate_matrix` slot. |
| `grna_target_cells` | `dict[str, np.ndarray]` | Maps each gRNA target to the **0-based** indices (into `covariate_matrix`'s cell axis) of cells treated with that target. This is the "union" grna-integration-strategy convention: one entry per target, not per individual gRNA. |
| `pairs` | `pd.DataFrame` with columns `response_id`, `grna_target` | The QC-passed (gene, target) pairs to test. `pysceptre` does not run `assign_grnas()`/`run_qc()` itself -- feed it pairs that have already passed QC upstream. |
| `side` | `"left"` \| `"both"` \| `"right"` | Test sidedness, matching sceptre's own convention. Use `"left"` for expected-repression screens (e.g. CRISPRi enhancer knockdown), `"both"` for a two-sided test. |
| `resampling_mechanism` | `"crt"` \| `"permutations"` | Matches sceptre's own option for high-MOI data. The CRT (default) draws each target's synthetic treated set from that target's own fitted probabilities; permutations draw one set of random subsets, sized by the largest target, and reuse it for every target. **The choice is a real trade, and yours to make** -- see [Reproducibility](#reproducibility-and-incremental-analysis), because permutations cannot offer the invariance the CRT does. R pairs permutations with `B3 = 24999` against the CRT's `0`, so sampling is cheaper -- and the per-target logistic fit is skipped entirely, since only the CRT draws from it -- but the escalation batch is five times larger. |
| `resampling_approximation` | `"skew_normal"` \| `"no_approximation"` | `"skew_normal"` (default, matching sceptre): pairs whose initial empirical p-value (`B1=499` draws) is `<= 0.02` get a skew-normal tail fit from a further `B2=4999` draws, giving p-values far smaller than `1/(B1+1)` could resolve. `"no_approximation"` fits no curve and instead draws a third, larger empirical batch, sized by R's own rule: `B3 = ceil(mult * n_pairs / multiple_testing_alpha)`, `mult = 10` two-sided and `5` one-sided. That grows linearly in the number of pairs and is much slower -- see [Scope and limitations](#scope-and-limitations). Any other value raises `ValueError`. |
| `seed` | `int \| None` | Seeds the `numpy.random.Generator` used for all CRT draws in the run. Note this does **not** reproduce sceptre's own R/C++ RNG stream bit-for-bit (different algorithm and seeding scheme) -- see [Scope and limitations](#scope-and-limitations). |
| `target_chunk_size` | `int` | How many gRNA targets to fit and CRT-draw at once. An **upper bound, not a mandate** -- it is reduced automatically to respect `chunk_memory_gb`, so no value here can exhaust memory. Default `200`. |
| `n_jobs` | `int` | Workers for the per-pair tests, which are ~80% of the runtime. `1` (default) runs serially; a negative value uses every core. **Results do not depend on it** -- only the genes inside an already-drawn target chunk are distributed, so the resampling draws are made in the same order at any worker count, and output is bit-identical. Processes on Linux, threads elsewhere (`fork` after macOS's Accelerate BLAS can deadlock), so the ceiling is lower off Linux. Memory grows by about one gene's working arrays per worker, not by `chunk_memory_gb` per worker. |
| `chunk_memory_gb` | `float` | Budget for the arrays a *chunk* holds, which sizes how many genes or targets are processed together. **Not** a cap on the process's memory -- the input, retained state and allocator overhead sit outside it. **You should not normally need to change this.** The default is both the fastest and the leanest setting measured: a larger budget produces chunks past the point where batching still pays, costing memory for no throughput (4 GB gave 8.42 GB peak against 3.78 GB at 1 GB, for the same runtime). Default `1.0`. |

**Returns** a `pd.DataFrame`, one row per input pair, with columns:

| Column | Meaning |
|---|---|
| `response_id`, `grna_target` | Echoed from `pairs`. |
| `p_value` | The test p-value (see `stage`). |
| `fold_change` | Estimated fold change of the treated group vs. complement control (deterministic given the data -- no resampling randomness). |
| `pct_change_es` | Effect size as a percent change from baseline, `(fold_change - 1) * 100`. Named `_es` rather than `pct_change` because `DataFrame.pct_change` is a pandas *method*: `result.pct_change` would silently return the method instead of the column, and comparisons against it fail with a confusing `TypeError` rather than a `KeyError`. |
| `pct_change_es_ci_low`, `pct_change_es_ci_high` | Two-sided 95% Wald interval on `pct_change_es`, from `se_fold_change`. Deterministic, unlike the resampled p-value, so it is a useful independent read on pairs sitting near a significance threshold. Note it is a normal approximation, not sceptre's test. |
| `log2(fold_change)` | Not returned -- a pure transform of a column already present. Compute it if you need it. |
| `se_fold_change` | Standard error of `fold_change`, on the same ratio scale, so `fold_change ± 1.96 × se_fold_change` is a Wald interval around 1 (no effect). Deterministic, unlike the resampled p-value. |
| `z_orig` | The observed test statistic (before resampling). |
| `stage` | `1` = reported from the initial `B1=499`-draw empirical p-value (not significant enough to escalate). `2` = escalated to a skew-normal tail fit on `B2=4999` further draws. `3` = an empirical p-value from a further batch, reached either because `resampling_approximation="no_approximation"` (so no curve is fit, and the `B3` draws are used) or because a skew-normal fit was attempted and *rejected* (see `check_sn_tail`/`check_for_outliers` in `test_statistic/skew_normal.py`), in which case the already-drawn `B2=4999` statistics are used. |

### `pysceptre.pipeline.api.run_calibration_check`

Runs the discovery test over **synthetic negative-control targets**, built by
regrouping individual non-targeting (NTC) gRNAs. No target is real, so a
correctly calibrated method returns p-values uniform on (0, 1) -- and that
uniformity, not agreement with any other implementation, is what the check
measures.

```python
from pysceptre import run_calibration_check

calib = run_calibration_check(
    response_matrix=response_matrix,
    gene_ids=gene_ids,
    covariate_matrix=covariate_matrix,
    ntc_grna_cells=ntc_grna_cells,   # dict[NTC gRNA id -> 0-based cell indices]
    n_calibration_pairs=len(pairs),
    calibration_group_size=15,
    n_nonzero_trt_thresh=7,
    n_nonzero_cntrl_thresh=7,
    side="both",
    seed=0,
)
```

Returns the same columns as `run_discovery_analysis`, minus any `pass_qc`
column -- see below for why there isn't one.

| Argument | Notes |
|---|---|
| `ntc_grna_cells` | `dict[NTC gRNA id -> 0-based cell indices]`. Keyed by **individual gRNA**, not by target. A target-keyed mapping collapses every NTC into one entry, and in sceptre's own object omits them entirely (`"non-targeting"` is not a key in `grna_group_idxs`), leaving nothing to regroup. |
| `n_calibration_pairs` | How many pairs to test. R defaults this to the number of discovery pairs that passed QC. |
| `calibration_group_size` | NTC gRNAs per synthetic target. R's default is the median gRNAs per real target, capped at the NTC count; that median is not derivable from these arguments, so it is required here. |
| `n_nonzero_trt_thresh`, `n_nonzero_cntrl_thresh` | Pairwise QC thresholds, sceptre's defaults being `7`. Prefer your object's own values. |
| `pass_qc_rate` | R's `p_hat`, the fraction of discovery pairs clearing QC, which sizes how many synthetic groups get built. Only matters when the group count is above its floor of 100 -- but there it is decisive. |
| `negative_control_pairs` | Test exactly these pairs instead of constructing any, with `grna_target` entries being `&`-joined NTC gRNA ids. This is how you compare against an R result pair-by-pair. |

**QC works differently here, deliberately.** A discovery result reports QC
failures in-band (`pass_qc = False`, NaN p-value). A calibration check
*constructs* its pairs and only ever samples combinations that already clear
the thresholds, so every returned row passes and there is no `pass_qc`
column. Pairwise nonzero-count filtering is therefore inseparable from
building the pairs; cell-level and gRNA-level QC remain out of scope.

**R's own pair selection is not reproducible.** Nothing in sceptre's
calibration path calls `set.seed` and `sceptre_object` has no seed slot, so
re-running R gives a different pair set. Comparing pair-by-pair against an R
result means passing R's pairs back in via `negative_control_pairs`.

### `pysceptre.pipeline.api.run_power_check`

Runs the discovery test over **positive controls** -- pairs where an effect
is expected, usually a gRNA against the gene's own TSS. Where the calibration
check asks whether the pipeline invents effects it should not, this asks
whether it recovers effects it should. The two are read together.

```python
from pysceptre import run_power_check

power = run_power_check(
    response_matrix=response_matrix,
    gene_ids=gene_ids,
    covariate_matrix=covariate_matrix,
    grna_target_cells=grna_target_cells,
    positive_control_pairs=positive_control_pairs,  # response_id, grna_target
    side="both",
    seed=0,
)
```

| Argument | Notes |
|---|---|
| `positive_control_pairs` | The pairs to test. **Supply these**: which target perturbs which gene is a claim only the experiment can make. Omitted, sceptre's name-matching rule is used -- a target that is itself a gene id pairs with that gene -- which works when targets are named after genes and finds *nothing* when they are named after genomic intervals. On a real screen of the latter kind it matched 0 of 3,071 targets, so that case raises rather than quietly returning an empty result. |
| `n_nonzero_trt_thresh`, `n_nonzero_cntrl_thresh` | Pairwise QC thresholds. |

**QC is reported, not filtered** -- the opposite of the calibration check.
The result has one row per supplied pair, with `pass_qc`, `n_nonzero_trt` and
`n_nonzero_cntrl`, and NaN results where a pair did not meet the thresholds.
Dropping those would overstate power by hiding exactly the controls the
screen had too few cells to test.

**No multiple-testing correction is applied**, matching R, which returns no
`significant` column here. These are a diagnostic rather than discoveries.

### `pysceptre.compute_power`

A different question from the three above, and the only entry point here that
does not come from sceptre. `run_power_check` runs the real test on pairs
where an effect is expected; this estimates, in closed form and with no
resampling, what the screen *could* have detected: if this element really did
reduce this gene by X%, how likely was this screen to call it?

It is a port of
[PerturbPlan](https://github.com/Katsevich-Lab/perturbplan)'s
`compute_power_posthoc()` (MIT -- see
[THIRD_PARTY_LICENSES](https://github.com/broadinstitute/pysceptre/blob/main/THIRD_PARTY_LICENSES)),
validated against that package's own output to a relative 1e-9.

```python
from pysceptre import compute_power

power = compute_power(
    discovery_pairs,             # grna_target, response_id
    cells_per_grna,              # grna_id, grna_target, num_cells
    baseline_expression_stats,   # response_id, expression_mean, expression_size
    fold_change_mean=0.85,       # a MULTIPLIER: a 15% knockdown is 0.85
    fold_change_sd=0.13,
    cutoff=alpha / 2,            # this screen's own nominal threshold
    num_total_cells=n_cells,
    side="left",
)
```

| Argument | Notes |
|---|---|
| `cells_per_grna` | One row per **individual** gRNA. Per-gRNA granularity is required, not a convenience: the across-gRNA variance term needs the sum of *squared* per-gRNA counts, which no target-level total can supply. Note that the treated count is the **sum** over a target's gRNAs, not the size of the union of perturbed cells -- the two differ in high MOI, and the sum is what was validated. |
| `baseline_expression_stats` | `expression_size` is the NB size, i.e. theta = `1 / dispersion`, not the dispersion. The validated `expression_mean` is the size-factor-normalised mean. |
| `cutoff`, `fold_change_mean`, `fold_change_sd` | All three are **required, with no defaults**, and `cutoff` must be the analysed screen's own nominal threshold rather than a borrowed one. |
| `n_nonzero_trt_thresh`, `n_nonzero_cntrl_thresh` | Default to `0`, not to sceptre's `7`, which makes the QC factor 1 -- correct for pairs that already passed QC. Raise them to score pairs that were never tested. |

Because none of its inputs is a property of the *pair* -- expression and
dispersion belong to the gene, cell counts to the element -- it also answers
the question for pairs the screen never tested.

Read the accuracy limits in
[Design decisions](https://broadinstitute.github.io/pysceptre/design/#analytical-per-pair-power)
before using it on a single pair: the estimate is good enough to plan a screen
and to triage its negatives, not to close a question about one element-gene
pair.

<!-- --8<-- [start:tutorial-power] -->

#### Reporting power alongside a discovery result

`run_discovery_analysis` returns no power column, and `compute_power` is not
wired into it. That is deliberate: the estimator needs inputs the discovery
path has no business knowing about, and a wrapper is where a wrong cell set or
the per-target union instead of the per-gRNA sum would creep in unseen. The
join is five steps, and each one is a place to check you meant it.

```python
from pysceptre import compute_power, run_discovery_analysis
from pysceptre.analytical_power import (
    baseline_expression_stats_from_fits,
    bh_nominal_cutoff,
    cells_per_grna_from_assignments,
)
from pysceptre.pipeline.discovery import fit_all_genes

ALPHA = 0.1

# 1. the analysis
result = run_discovery_analysis(
    response_matrix, gene_ids, covariate_matrix, grna_target_cells, pairs,
    side="both", multiple_testing_alpha=ALPHA, seed=0,   # "both" is the default
)

# 2. the threshold this run actually applied, then halved. side="both" gives a
#    two-sided p-value, and compute_power's default side="left" wants the
#    one-sided threshold that corresponds to it. Do NOT halve if you ran the
#    discovery one-sided: those p-values are already one-sided.
cutoff = bh_nominal_cutoff(result["p_value"], alpha=ALPHA) / 2

# 3. the estimator's inputs, both on sceptre's own scale
fits = fit_all_genes(response_matrix, gene_ids, covariate_matrix)
baseline = baseline_expression_stats_from_fits(covariate_matrix, fits)
cells_per_grna = cells_per_grna_from_assignments(
    grna_target_data_frame, targeting_grna_cells,   # both from the export
)

# 4. power for the same pairs, at a 15% knockdown
power = compute_power(
    pairs, cells_per_grna, baseline,
    fold_change_mean=0.85, fold_change_sd=0.13,
    cutoff=cutoff, num_total_cells=covariate_matrix.shape[0],
)

# 5. one table
report = result.merge(
    power[["response_id", "grna_target", "power"]],
    on=["response_id", "grna_target"], how="left",
)
```

Four things worth knowing before you run it.

**Step 3 refits every gene, and the discovery run already did.** It discards
its fits rather than returning them, so this repeats the most expensive part
of the analysis. If that matters, call `fit_all_genes` once yourself and reuse
the result; the fit depends only on the counts and the covariate matrix, not
on the pairs.

**Step 2 raises if nothing is significant**, rather than returning a threshold
that would hand every pair zero power. On a screen with no discoveries, pass
`cutoff` yourself from a plain alpha and say so in whatever you report.

**Step 2 needs a real pair list to mean anything.** The threshold is the
largest p-value BH calls significant, so on a handful of pairs it is simply
whatever the smallest p-value happens to be: five pairs containing one strong
hit gave `1.2e-215`, and every power estimate came back 0. That is the rule R
uses, reproduced, not a defect -- but it only behaves like a threshold at
screen scale, where day0's 34,886 pairs give `7.26e-4`. On a toy example pass
`cutoff` explicitly.

**The halving in step 2 belongs to the two-sided run, not to
`compute_power`.** `side="both"` is sceptre's default and pysceptre's, and it
is the only side the discovery path is validated against R on. Its p-values
are two-sided, so the one-sided threshold `compute_power(side="left")` wants
is half of it. Run the discovery one-sided and its p-values are already
one-sided: use the threshold **unhalved**. Halving twice is the
easiest way to get a confidently wrong power estimate here, and it is
invisible in the output.

**The pairs that failed QC are the interesting ones.** They carry a NaN
p-value because they were never tested, and they are exactly where a power
estimate says something a p-value cannot. Pass them to `compute_power` too,
with `n_nonzero_trt_thresh` and `n_nonzero_cntrl_thresh` set to the
thresholds the QC used, so the estimate includes the probability the pair
would have failed QC at all.

<!-- --8<-- [end:tutorial-power] -->

### Lower-level building blocks

`run_discovery_analysis` is a thin wrapper around
`pysceptre.pipeline.discovery.run_discovery_ntcells_complement`, which
exposes a couple of additional knobs not surfaced at the top level (fixed
`B1=499, B2=4999, B3=0`, `fit_parametric_curve: bool`, `side_code: int`
instead of `side: str`). Most users won't need to go lower than
`run_discovery_analysis`, but each pipeline stage is also independently
importable and unit-tested, for anyone extending or debugging the pipeline:

| Module | What it does |
|---|---|
| `pysceptre.glm.irls` | Batched IRLS for Poisson (log link, gene fits) and binomial (logit link, gRNA-target fits) GLMs -- `fit_poisson_glm_batch`, `fit_binomial_glm_batch`. Ports R's `stats::glm.fit` algorithm and convergence criteria exactly. |
| `pysceptre.glm.nb_theta` | Negative-binomial dispersion (`theta`) estimation given a fitted mean -- `estimate_theta`. Exact port of sceptre's `estimate_theta`. |
| `pysceptre.precompute.pieces` | Per-gene precomputation pieces (the `D` matrix and friends) needed by the test statistic -- `compute_precomputation_pieces`. |
| `pysceptre.crt.sampler` | The CRT resampling draw itself -- `crt_index_sampler_fast` (sparse, real-scale-feasible; numba-accelerated if available) and `crt_index_sampler_naive` (dense reference implementation, used only for cross-checking in tests). |
| `pysceptre.test_statistic.score_stat` | The O(n_treated)-per-resample score-type test statistic -- `compute_observed_full_statistic`, `compute_null_full_statistics`. |
| `pysceptre.test_statistic.empirical_p` | Empirical p-value from a null distribution -- `compute_empirical_p_value`. |
| `pysceptre.test_statistic.skew_normal` | Skew-normal tail-fit escalation -- `fit_and_evaluate_skew_normal`. Exact port of `fit_skew_normal_funct`/`check_sn_tail`/`check_for_outliers`/`fit_and_evaluate_skew_normal`. |
| `pysceptre.test_statistic.fold_change` | Fold-change estimation -- `estimate_log_fold_change`. |
| `pysceptre.test_statistic.resampling` | Ties the above into the `B1 -> B2 -> B3` staged escalation for one pair -- `run_low_level_test_full`. |
| `pysceptre.pipeline.discovery` | Orchestration: `fit_all_genes`, `fit_all_targets`, `run_discovery_ntcells_complement`. |

## Reproducibility and incremental analysis

<!-- --8<-- [start:reproducibility] -->

A pair's result depends on `(seed, response_id, grna_target)` and the data,
and on nothing else about the run. Each target draws its CRT resamples from
its own stream keyed on the target's *name*, so results do not depend on the
order targets are processed, the chunk size, the worker count, or **which
other pairs are in the analysis**.

That makes an incremental workflow safe: run a subset, check it, then run the
full set and reuse what you already have. The pairs in common come back
identical rather than merely similar.

**This applies to the CRT, and cannot apply to permutations.** With
`resampling_mechanism="permutations"` the draws are made once and shared by
every target, sized by the largest target present, so adding a target bigger
than the current largest changes every result in the run. That is inherent to
sharing one draw set -- it is what makes permutations cheap -- not a defect
that could be fixed. Permutation runs remain fully deterministic for a fixed
pair list, and independent of chunk size and worker count; they are simply
not invariant to changing the analysis. Use the CRT if you intend to extend
an analysis and reuse earlier results.

Two limits on the CRT path, both worth knowing exactly:

| change between runs | shared pairs |
|---|---|
| more or fewer targets or pairs; different order; different `target_chunk_size` or `n_jobs` | **bitwise identical** |
| a different set of *genes* | identical to ~1e-16 |

The second is batched linear algebra, not randomness: gene GLMs are fitted in
batches, and changing the batch width changes the order of floating-point
additions. Measured at one unit in the last place (1.11e-16 on `fold_change`,
3 of 72 pairs) on OpenBLAS, and exactly zero on Apple Accelerate. It cannot
be removed without fitting genes one at a time, which is the batching this
package exists to do.

**Adjusted p-values are a different matter and must be recomputed.** A
Benjamini-Hochberg adjustment depends on every p-value in the set, so adding
pairs changes the adjusted value, and possibly the call, for pairs already
tested. That is multiple testing behaving correctly. `run_discovery_analysis`
returns raw p-values and applies no correction, so cache those and run the
adjustment over the union each time.

<!-- --8<-- [end:reproducibility] -->

## Scope and limitations

<!-- --8<-- [start:scope] -->

- **Complement control group only, high-MOI.** Non-complement control groups
  and low-MOI are out of scope.
- **The CRT is the validated resampling mechanism; permutations are
  implemented but not validated against R.**
  `resampling_mechanism="permutations"` works and is exercised end to end, but
  its tests check internal consistency, R's `B3` sizing rule and usability --
  not value-for-value agreement with R, which the CRT path does check. Use it
  knowing that, and don't report a permutation result as R-validated.
- **No `assign_grnas()` / `run_qc()`**, with one carve-out: the calibration
  check applies the *pairwise* nonzero-count thresholds, because it builds its
  own pairs and cannot select them otherwise. Cell-level and gRNA-level QC
  are still out of scope. Feed `run_discovery_analysis` pairs
  that have already passed QC (e.g. from a real `sceptre_object`'s
  `@discovery_pairs_with_info`, filtered to `pass_qc == TRUE`). This package
  is the statistical engine only.
- **No formula DSL.** `covariate_matrix` must already be a plain numeric
  design matrix; there's no `model.matrix()`-equivalent formula parser here.
- **RNG is not bit-for-bit reproducible against R.** sceptre seeds
  `boost::mt19937` with a fixed literal seed on every resampling call;
  `pysceptre` uses `numpy.random.Generator` with a different algorithm and
  seeding scheme. Validated against R by matching *distributions* (the
  resulting null-statistic and p-value distributions), not exact draws --
  see [Validation](https://github.com/broadinstitute/pysceptre/blob/main/README.md#validation).
- **`B1`/`B2`/`B3` are not exposed at the top-level API**, but they are no
  longer fixed: `run_discovery_analysis` derives them from
  `resampling_approximation` exactly as R does -- `B1=499` always, then
  `(B2, B3) = (4999, 0)` for `skew_normal` and `(0, ceil(mult * n_pairs /
  multiple_testing_alpha))` for `no_approximation`. `B3=0` on the
  `skew_normal` path is parity with R, which only uses `B3=24999` for the
  `permutations` mechanism -- which this package does implement, and sizes the
  same way.
  `run_discovery_ntcells_complement` and `run_low_level_test_full` accept all
  three directly if you need to override them.
- **`no_approximation` is expensive, and can be coarser at small scale.** Its
  `B3` grows linearly in pair count: a 33,066-pair one-sided run needs
  1,653,300 draws *per target* (5.2 GB), so the chunk collapses to a
  single target and the run warns that it may still exhaust memory.
  Conversely, at 4 pairs one-sided `B3 = 200`, *below* `B1 = 499`. Both are
  R's behavior, reproduced rather than corrected.
- **All three gRNA integration strategies.** `"union"` (the default),
  `"singleton"` and `"bonferroni"`, matching sceptre's own option. Under the
  latter two, `grna_target_cells` is keyed by guide rather than by target and
  `grna_target_data_frame` supplies the design; the pair expansion is
  validated against sceptre's own on a real screen, 36,450 pairs to 515,972.
  Nothing statistical differs between them -- only which cells count as
  treated, and what happens to the results afterwards.
- **`compute_power` is not a sceptre path, and its limits are its
  own.** It estimates in closed form what a screen *could* have detected,
  which is a different question from the three analyses above, and it is a
  port of [PerturbPlan](https://github.com/Katsevich-Lab/perturbplan) (MIT)
  rather than of sceptre. It is validated against that package's own output,
  not against sceptre's. Within it: complement control group only, an explicit
  `cutoff` only, per-gRNA cell counts required rather than the per-target
  union, and no minimum-detectable-effect-size path. It is accurate enough to
  plan a screen and triage its negatives, and **not** accurate enough to
  settle a question about one element-gene pair -- see
  [Design decisions](https://broadinstitute.github.io/pysceptre/design/#analytical-per-pair-power).
- **Two ways to build its baseline statistics, and they are not
  interchangeable.** `baseline_expression_stats_from_fits` takes the mean and
  theta from the same negative-binomial fit the discovery test uses, so the
  estimate and the test it predicts are on one expression scale. **Prefer it.**
  `baseline_expression_stats` instead computes a size-factor-normalised mean,
  which is what the published comparison used; measured on one screen that
  mean sits about **16 % below** the scale sceptre's own model works on, which
  makes the estimate conservative rather than wrong-shaped. Use it only to
  reproduce those numbers.

<!-- --8<-- [end:scope] -->

## Performance

<!-- --8<-- [start:performance] -->

Benchmarked on the `day0_grna20` single-cell CRISPR screen: 567,690 cells
after QC, 237 genes appearing in pairs, 3,071 gRNA targets, 34,886 QC-passing
pairs, two-sided, CRT resampling. The headline figures are below, each
reported with the hardware it was measured on; the full tables and the
comparison against R sceptre are in the manuscript repository,
`broadinstitute/pysceptre-paper`, which stays private until the preprint is
posted.

Reproduce with `scripts/benchmark_vs_r.R` (R side, which also exports the
exact inputs) and `scripts/benchmark_pysceptre.py` (pysceptre side). The
pinned R-plus-pysceptre image those numbers were measured in lives with the
comparison itself, in the private development archive, rather than here:
`docker/` in this repository holds only the pysceptre runtime image, which
carries no R.

### Sizing the machine

**The two mechanisms want different machines, so `n_jobs` is worth setting
deliberately rather than to the core count.** Measured on day0, one process
per configuration, Apple M4 Max; `cores` is `(user + sys) / wall`, the mean
number actually busy.

| `resampling_mechanism` | `n_jobs` | wall | peak RSS | cores |
|---|---|---|---|---|
| `"crt"` (default) | 1 | 555.6 s | 3.63 GB | 1.09 |
| | 4 | 165.3 s | 6.20 GB | 5.00 |
| | 8 | **138.2 s** | 6.81 GB | 6.74 |
| `"permutations"` | 1 | 125.3 s | 3.47 GB | 1.06 |
| | 8 | **30.8 s** | 5.85 GB | 6.53 |

**Both mechanisms now use a large machine, so give them one.** The CRT
reaches 6.74 of 8 cores and gains 1.20x from four workers to eight;
permutations reach 6.53 and 30.8 s. Neither is memory-constrained at these
sizes -- the largest peak here is 6.8 GB, inside what a 4-vCPU cloud
instance ships with, and cloud machine types bundle memory with cores
anyway, so asking for fewer cores buys less memory rather than a cheaper
machine at the same memory.

*This guidance changed.* The CRT used to cap near 3.90 cores with four
workers within 4.4% of eight, because its per-chunk logistic fit ran on one
thread and nothing else could proceed past it. Chunks are now prepared
several deep, so several fits run concurrently, and the ceiling moved. If
you tuned `n_jobs` down for the CRT on the old advice, undo it.

**Memory rises with workers and with pipeline depth, not with
`chunk_memory_gb`.** Each worker holds one gene's working arrays and the
pipeline holds `_PREFETCH_DEPTH` chunks of draws, so peak tracks those two;
the chunk budget is a weaker lever than it looks. Raising it from 1 GB to
8 GB buys 13% on the CRT for 2.3x the memory, and past that it gets
*slower* -- 16 GB ran 369 s against 8 GB's 353 s. Leave it alone unless you
have measured otherwise on your own data. Results are identical at every
setting either way.

<!-- --8<-- [end:performance] -->

## Validation

Two independent validation layers, both against a real, installed `sceptre`
R package (pinned upstream commit), not just internal self-consistency:

1. **Synthetic ground truth** (`tests/validation/`): small hand-built
   matrices run through both sceptre's internal (`:::`) R functions
   directly and the corresponding pysceptre function, compared value-for-
   value. Covers GLM fits, NB dispersion estimation, precomputation pieces,
   the test statistic, empirical p-values, skew-normal fitting, and the CRT
   sampler's *distribution* (not exact draws -- see
   [Scope and limitations](#scope-and-limitations)). Run via
   `scripts/dump_r_ground_truth.R` + `pytest tests/validation/`.

2. **Real dataset: discovery analysis** (day0_grna20: 567,690 cells /
   292 genes / 3,026 targets / 34,886 QC-passed pairs). The same screen was
   run through `pysceptre.run_discovery_analysis` and R's
   `sceptre::run_discovery_analysis()`, with R exporting the exact object it
   analysed so the two cannot disagree about their inputs:

   | Metric | Value |
   |---|---|
   | Fold-change agreement (Pearson r) | 1.000000 (max difference 2.6e-12) |
   | p-value agreement (Spearman rho) | 0.9865 |
   | p-value agreement (Pearson r on -log10 p) | 0.9940 |

   Significance calls at BH 0.1, with R as the reference:

   | | R: significant | R: not |
   |---|---|---|
   | **pysceptre: significant** | 251 | 10 |
   | **pysceptre: not significant** | 3 | 34,622 |

   Sensitivity 0.9882, specificity 0.9997, 13 discordant pairs of 34,886.
   Reported as a 2x2 table rather than a single index, because a Jaccard
   coefficient collapses it and cannot distinguish calling extra hits from
   missing them, which are different failures.

   Fold change carries the weight here: no resampling enters it, so agreement
   to 2.6e-12 establishes that both implementations fit the same models to
   the same cells. The p-value correlation is then bounded below by Monte
   Carlo noise -- pysceptre agrees with R about as closely as it agrees with
   itself across seeds. Every one of the 13 discordant pairs has an interval
   excluding no effect and a p-value within a factor of four of the
   threshold, so the disagreement is about which side of a cutoff a resampled
   p-value fell on, never about whether an effect was detected.

3. **Calibration check** (day0_grna20: 567,690 cells / 292 genes / 2,031
   non-targeting gRNAs / 34,886 negative-control pairs), with R's own pairs
   injected so both sides test identical ones:

   | Metric | Value |
   |---|---|
   | Pairs merged | 34,886 / 34,886 |
   | Fold-change agreement (Pearson r) | 1.0 (max abs difference 4.5e-11) |
   | p-value agreement (Spearman rho) | 0.9863 |
   | KS vs `U(0,1)` | 0.02592 (pysceptre) vs 0.02594 (R) |
   | False discoveries, BH at 0.1 | 7 (pysceptre) vs 6 (R) |

   Before any p-value, a deterministic checkpoint: rebuilding R's synthetic
   groups from their names and recomputing the pairwise counts reproduces R
   **exactly** -- `n_nonzero_trt` and `n_nonzero_cntrl` both 34,886/34,886,
   zero discrepancies over ~70,000 integer comparisons.

   The negative-control p-values are **not perfectly uniform** -- KS 0.0259
   against a 5% critical value of 0.0073, with 6.1% below 0.05 where 5% is
   expected. R deviates by the same amount to four decimal places, so this is
   sceptre's behaviour on this data and not an artefact of the port; pysceptre
   reproduces the reference's mild anti-conservatism rather than adding any.

4. **Power check** (day0_grna20, 266 positive-control pairs), against R's
   stored `power_result`:

   | Metric | Value |
   |---|---|
   | rows | 266 / 266 |
   | `pass_qc` agreement | 266 / 266 |
   | `n_nonzero_trt`, `n_nonzero_cntrl` | exact, 266 / 266 each |
   | Fold-change agreement (Pearson r) | 1.000000 (max difference 1.3e-12) |
   | p-value agreement (Spearman rho), 246 testable pairs | 0.9892 |
   | median effect | -35.2% |

   The 20 pairs that fail pairwise QC are reported with NaN results by both
   implementations, and both agree on exactly which 20. A median effect of
   -35.2% is the check working: positive controls target a gene's own TSS,
   so they should repress it.

## License and attribution

`pysceptre` is licensed under the GNU General Public License v3.0 only
(SPDX: `GPL-3.0-only`) -- full text in [LICENSE](LICENSE). That is inherited,
not chosen: this package ports algorithms and formulas from
[`sceptre`](https://github.com/Katsevich-Lab/sceptre), which is `GPL-3` --
version 3 exactly, not "or later". Each module's docstring names what it
ports and from which upstream file.

The original SCEPTRE is by Timothy Barry and Eugene Katsevich. The
statistical method is theirs, not ours -- if you use this in published work,
cite `sceptre`: Barry et al. (2021), *SCEPTRE improves calibration and
sensitivity in single-cell CRISPR screen analysis*, Genome Biology 22(1),
[doi:10.1186/s13059-021-02545-2](https://doi.org/10.1186/s13059-021-02545-2).
