User guide¶

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

A low-MOI screen takes one more input, each non-targeting gRNA's cells, and moi="low" selects sceptre's low-MOI defaults: the cells carrying a non-targeting gRNA as the control group, and permutation resampling.

result = run_discovery_analysis(
    response_matrix, gene_ids, covariate_matrix, grna_target_cells, pairs,
    moi="low",
    ntc_grna_cells=ntc_grna_cells,        # dict[NT gRNA id -> 0-based cell indices]
    seed=0,
)

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.

The other two sceptre analyses¶

The calibration check and the power check take the same inputs as the discovery analysis and return the same shape of frame. They differ in where their pairs come from, and that difference shows up in the result:

	pairs	QC failures
`run_discovery_analysis`	you supply them	reported in-band: `pass_qc = False`, `p_value = NaN`
`run_calibration_check`	constructed from non-targeting gRNAs	cannot occur -- pairs are chosen to pass, so there is no `pass_qc` column
`run_power_check`	you supply positive controls	reported in-band, like discovery

The asymmetry is deliberate. A positive control is a claim about one specific pair, so silently dropping the ones with too few cells would overstate power; a calibration pair is interchangeable, so one that would fail is simply never drawn. See Design decisions.

Signatures and parameters for all three are in the API reference.

Low MOI¶

In a low-MOI screen each cell carries at most one gRNA, so the cells carrying a non-targeting (NT) gRNA are untreated and can be the control group on their own. sceptre calls that the NT-cells control group and makes it the low-MOI default, together with permutation resampling. moi="low" selects both, as R's set_analysis_parameters does, and either can be overridden:

`moi`	default `control_group`	default `resampling_mechanism`
`"high"` (default)	`"complement"`	`"crt"`
`"low"`	`"nt_cells"`	`"permutations"`

control_group="nt_cells" is refused in high MOI, where R would quietly use the complement instead. Low MOI with control_group="complement" runs exactly the high-MOI engine.

What the NT cells need. Pass ntc_grna_cells, each individual NT gRNA's cells. The analysis assumes what sceptre's low-MOI QC guarantees: every cell carries at most one gRNA, because run_qc removes cells with zero or two or more. pysceptre does not run that QC, so it refuses what it can see breaking the assumption -- a cell under two NT gRNAs, or a tested target sharing a cell with the NT cells -- and those cells have to be removed upstream. If pysceptre made the assignment, pysceptre.assignment.cells_w_zero_or_twoplus_grnas, or a maximum assignment's field of that name, lists the cells sceptre's QC would remove. A cell in two targets is accepted: one guide can belong to two overlapping elements.

Each pair gets its own fit. Against the NT cells a pair is tested on its target's cells together with the NT cells, so the gene's GLM is refit for every pair, as in R, rather than once per gene. The fits are independent and spread across workers with n_jobs. A covariate constant on one target's cells and on every NT cell makes that pair's GLM unfittable even when the full design is full rank; R stops there, and pysceptre refuses it before fitting anything, naming the target.

The calibration check runs on the NT cells alone, as in R: a synthetic target is tested against the rest of the NT cells, so at least two NT gRNAs are needed. The power check counts a pair's control cells as the gene's nonzero NT cells.

The covariates are yours to build. sceptre's default formula leaves out the gRNA-count covariates (grna_n_nonzero, grna_n_umis) in low MOI and keeps them in high MOI. There is no formula DSL here, so build covariate_matrix the same way to match an R run.

The reasoning, and what it was checked against, is in Design decisions.

A worked example is in Tutorials.

Choosing what counts as one treated unit¶

run_discovery_analysis takes grna_integration_strategy, matching sceptre's own option:

	one test per	result identifies
`"union"` (default)	gRNA target, on the union of its guides' cells	`grna_target`
`"singleton"`	individual guide	`grna_id` and `grna_target`
`"bonferroni"`	guide, then aggregated back to the target	`grna_target`

Nothing statistical differs between them. The same CRT, the same score statistic, the same escalation; only which cells count as treated, and in bonferroni's case what happens to the results afterwards. The two per-guide strategies want grna_target_cells keyed by guide and the design table in grna_target_data_frame, and they cost roughly a guides-per-target multiple of the work. Tutorials has a worked example.

A fourth entry point, which is not sceptre's¶

compute_power estimates in closed form what a screen could have detected, per pair and without simulation, so it speaks to the pairs a p-value cannot: the ones that failed QC and were never tested. It is a port of PerturbPlan rather than of sceptre and is validated against that package instead, so read its limits in Design decisions before using it on a single pair. run_discovery_analysis does not call it, and Tutorials shows the join. matched_expression_stats is an optional per-pair input that carries the covariates into the estimate; see Design decisions.

A fifth, also not sceptre's¶

run_specificity_check asks how many of the links a screen discovered are more than background: how often a test is called without regulation, measured on tests between elements and genes on other chromosomes, matched gene for gene and bootstrapped over elements. It needs a discovery result on the cis pairs, one from the same test on those elements against genes elsewhere, and positions for both. It has no R counterpart; read Design decisions for what it has been checked against.

Assigning gRNAs to cells¶

The analyses take each target's cells as given. assign_grnas computes them from raw gRNA UMI counts, like sceptre's function of the same name, with one of four methods: sceptre's mixture, thresholding and maximum methods, each validated against sceptre's own, or fishash's one-sided Fisher test, which is not sceptre's and is validated against fishash's R. Left at its default it takes sceptre's choice, the maximum method in low MOI and the mixture in high MOI. Each method is also a function of its own: assign_grnas_mixture, assign_grnas_thresholding, assign_grnas_maximum and assign_grnas_fishash. pysceptre.assignment.cells_by_target turns any result into the grna_target_cells and ntc_grna_cells the analyses take.

In low MOI, sceptre's QC then removes cells, and that removal is left to you: a maximum assignment returns the cells it flags, and after any other method pysceptre.assignment.cells_w_zero_or_twoplus_grnas lists the cells assigned no gRNA or two or more. Read Design decisions for what each method normalizes for, and Thresholding and maximum for those two methods' edge cases.

If a design matrix is refused¶

covariate_matrix is checked before any fitting and a rank-deficient design is refused by name, because there is no formula DSL here to catch an aliased contrast for you. The usual causes are a factor's full dummy set alongside an intercept, a covariate included twice on different scales, or a column that is constant within another's levels. Hold out a reference level, or drop the column the error names.

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.

Against the NT cells the same split holds. Each pair is fit on its own target's cells and the NT cells, so under the CRT a pair's result depends on its gene, its target, the NT cells and the seed, and not on the rest of the analysis. Permutation draws there are sized by the number of NT cells and by the smallest and largest target, so adding a target outside that range changes every result.

One exception, by construction: under resampling_approximation="no_approximation" the third batch is sized by the number of pairs, so changing the pair list re-draws every resample, on both mechanisms.

The limits on the CRT path, worth knowing exactly:

change between runs	shared pairs
more or fewer pairs of the same targets, or genes; a different pair order; different `n_jobs`	bitwise identical
different target chunks: other targets, another order of `grna_target_cells`, or another `target_chunk_size`	bitwise identical with `glm_solver="structured"`; with the default `"outer"`, on some BLAS kernels a target's fitted probabilities differ in the last bits, which very rarely moves a resample
the covariate matrix row-major in one run and column-major in the other	the same to rounding, not bitwise, except as below
numba in one run and not in the other	the same to rounding, not always bitwise: the same resamples, some with their cells listed in another order
another CPU type	the same to rounding, not bitwise, except as below

None of these is randomness. Gene GLMs are fitted one gene at a time, so the set of genes never matters. Target GLMs are fitted a chunk at a time on the default route, and a batched BLAS product need not give one column the same last bits in every batch; glm_solver="structured" fits each target alone (Design decisions). BLAS rounds a product differently for the two memory layouts, and pysceptre uses the covariate matrix in the layout it is given, so give every run that must agree bit for bit the same layout, for example np.ascontiguousarray(covariate_matrix). Without numba the CRT sampler lists each resample's cells in another order, and a null statistic adds them up in that order. numpy and OpenBLAS choose their kernels by CPU model, so the same code gives different last bits on different machines (Design decisions).

The exception is a gene close to Poisson. In the two rows marked "except as below" the gene fits can themselves differ in their last bits (always on another CPU type; for another covariate layout, under the default glm_solver="outer", since "structured" fits from a row-major copy), and for a gene with little or no overdispersion which branch the dispersion estimate ends on can turn on those bits. Such a gene's theta, and with it its p-values, can then differ by much more than rounding (Design decisions).

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.

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.