pysceptre

A standalone Python port of the statistical engine behind sceptre'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 for exactly what is and isn't covered.

Why this exists

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.

Installation

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, 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.

Quick start

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 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 for a complete, runnable walkthrough (including how to build each input from scratch) and for guidance on picking target_chunk_size for your dataset.

API reference

pysceptre.pipeline.api.run_discovery_analysis

The one function most users need.

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

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

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.

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.

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.

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,
)

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's compute_power_posthoc() (MIT -- see THIRD_PARTY_LICENSES), validated against that package's own output to a relative 1e-9.

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",
)

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 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.

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.

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.

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:

Reproducibility and incremental analysis

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:

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.

Scope and limitations

• 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.
• 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 (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.
• 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.

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.

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.

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). 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. That is inherited, not chosen: this package ports algorithms and formulas from 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.