API reference

Generated from the source. Twelve functions make up the public API, and all twelve are importable from the top level (from pysceptre import run_discovery_analysis) as well as from their defining module.

Reading these signatures

pysceptre.pipeline.api is the documented path and the top-level re-export is a convenience; both are supported and both will keep working.

Public API

pysceptre.pipeline.api.run_discovery_analysis

run_discovery_analysis(response_matrix, gene_ids: list[str], covariate_matrix: ndarray, grna_target_cells: dict[str, ndarray], pairs: DataFrame, *, side: str = 'both', moi: str = 'high', control_group: str | None = None, ntc_grna_cells: dict[str, ndarray] | None = None, grna_integration_strategy: str = 'union', grna_target_data_frame: DataFrame | None = None, drop_duplicate_design_rows: bool = False, resampling_approximation: str = 'skew_normal', multiple_testing_alpha: float = 0.1, seed: int | None = None, target_chunk_size: int = _DEFAULT_TARGET_CHUNK_SIZE, chunk_memory_gb: float = _DEFAULT_CHUNK_MEMORY_GB, n_jobs: int = 1, resampling_mechanism: str | None = None, grna_target_weights: dict[str, ndarray] | None = None, glm_solver: str = 'outer') -> DataFrame

response_matrix: (n_genes, n_cells) dense ndarray or scipy.sparse matrix. gene_ids: row labels for response_matrix, in order. covariate_matrix: (n_cells, p) already formula-expanded design matrix. grna_target_cells: dict[target -> 0-based treated-cell indices]. pairs: DataFrame['response_id', 'grna_target'] -- QC-passed pairs to test. moi: "high" (default) or "low", sceptre's import_data(moi=). It sets the defaults of control_group and resampling_mechanism, as R's set_analysis_parameters does, and is the only way to reach the NT cells. control_group: "complement" or "nt_cells"; None takes the MOI's default, "complement" in high MOI and "nt_cells" in low MOI. With "nt_cells" each pair is tested on its target's cells against the NT cells only, and the gene's GLM is refit for every pair, as in R. ntc_grna_cells: dict[NTC gRNA id -> 0-based cell indices], required with control_group="nt_cells" and refused otherwise. The NT cells are their union, which must share no cell with a tested target, and no cell may sit under two NT gRNAs: sceptre's low-MOI QC removes every cell with more than one gRNA, and this path relies on that. multiple_testing_alpha: only used to size the no_approximation resampling budget, exactly as R's run_qc does. pysceptre does not apply any multiple-testing correction to the returned p-values. target_chunk_size: how many gRNA targets' logistic fits + CRT draws to batch/hold in memory at once (see pipeline/discovery.py) -- lower this if you hit memory pressure, raise it for a modest speed gain if you have memory to spare. resampling_mechanism: "crt" or "permutations"; None takes the MOI's default, "crt" in high MOI and "permutations" in low MOI, as in R. The CRT draws each target's synthetic treated set from that target's own fitted probabilities; permutations draw one set of random subsets, sized by the targets present, and reuse it for every target.

The choice is a real trade and is left to the caller.
**Permutations cannot be reproducible across a change of pair list**:
the shared draws are sized by the largest target present, so adding a
target bigger than the current largest moves every result in the run.
The CRT path has no such dependence -- each target seeds its own
stream from its own name -- so a subset can be analysed, checked and
extended with the earlier pairs reused unchanged. Permutations also
pair with `B3 = 24999` in R against the CRT's `0`, so sampling is
cheaper but the escalation batch is five times larger.

n_jobs: worker processes (Linux) or threads (elsewhere) used for the per-pair tests, which are ~80% of the runtime. 1 disables parallelism, a negative value uses every core. Results do not depend on it: only the genes within an already-drawn target chunk are distributed, so the resampling draws are made in the same order whatever the worker count. Memory grows by roughly one gene's working arrays per worker, not by chunk_memory_gb per worker. chunk_memory_gb: budget for the gRNA-target chunk's arrays. It sizes target chunking only; gene fits do not use it. With glm_solver="structured" it cannot change a result. With the default "outer" a chunk's targets are fit as one batch, and on some BLAS kernels a target's fitted probabilities differ in the last bits with its chunk-mates, which very rarely moves a resample (docs/design.md, "Target fits no longer depend on their chunk-mates"). 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 the fastest and leanest setting measured. Reductions to target_chunk_size are warned about; like any chunk width, they can change a result only under "outer", as above. grna_target_weights: {target: weights}, aligned with grna_target_cells, runs the dose test instead of sceptre's: each treated cell enters the score statistic with its weight in place of 1 (dose_weights builds both dicts from gRNA counts). Weights must be finite and positive; all-ones weights give sceptre's result exactly. Complement control group, CRT and grna_integration_strategy="union" only. Not from sceptre: docs/design.md, "The dose test". glm_solver: how the gene and target GLMs are fit. "outer" (the default) is 0.5.0's IRLS over the per-row outer products of the covariate matrix; "structured" forms the same normal equations from the design's indicator columns by per-level sums, which on a design with many one-hot covariates is far cheaper and needs no (n_cells, p^2) array. The same GLM, not the same bits: a gene whose counts are close to Poisson can get a different theta, as it can on another CPU. Complement control group only. See docs/design.md, "Indicator columns in the IRLS".

pysceptre.pipeline.api.run_calibration_check

run_calibration_check(response_matrix, gene_ids: list[str], covariate_matrix: ndarray, ntc_grna_cells: dict[str, ndarray], *, n_calibration_pairs: int, calibration_group_size: int, n_nonzero_trt_thresh: int = 7, n_nonzero_cntrl_thresh: int = 7, pass_qc_rate: float = 1.0, negative_control_pairs: DataFrame | None = None, side: str = 'both', moi: str = 'high', control_group: str | None = None, resampling_approximation: str = 'skew_normal', multiple_testing_alpha: float = 0.1, seed: int | None = None, target_chunk_size: int = _DEFAULT_TARGET_CHUNK_SIZE, chunk_memory_gb: float = _DEFAULT_CHUNK_MEMORY_GB, n_jobs: int = 1, resampling_mechanism: str | None = None, ntc_grna_weights: dict[str, ndarray] | None = None, glm_solver: str = 'outer') -> DataFrame

Run sceptre's calibration check: the discovery test over negative controls.

Synthetic negative-control targets are built by regrouping individual non-targeting gRNAs, then tested with the same engine as run_discovery_analysis. Because no target is real, a correctly calibrated method returns p-values that are uniform on (0, 1); that uniformity, not agreement with any other implementation, is what the check measures.

dict[NTC gRNA id -> 0-based cell indices]. This must be

keyed by individual gRNA, not by target -- a target-keyed mapping collapses every non-targeting gRNA into one entry (and in sceptre's own object, omits them entirely), leaving nothing to regroup.

moi, control_group, resampling_mechanism: as in run_discovery_analysis, with the same MOI-dependent defaults. With control_group="nt_cells" the whole check runs on the NT cells alone, as in R: the synthetic targets are tested against the remaining NT cells, so it needs at least two NT gRNAs and no cell under two of them. Under permutations the shared draws are as wide as the calibration_group_size largest NT gRNAs together, R's rule. n_calibration_pairs: how many pairs to test. R defaults this to the number of discovery pairs that passed QC. calibration_group_size: how many NTC gRNAs per synthetic target. R's default is the median number of gRNAs per real target, capped at the number of NTC gRNAs; that median is not derivable from this function's arguments, so it is required here and the dataset export records it. n_nonzero_trt_thresh / n_nonzero_cntrl_thresh: pairwise QC thresholds. Note these filter construction rather than being reported: every returned pair passes, so there is no pass_qc column, unlike a discovery result. negative_control_pairs: skip construction and test exactly these pairs, whose grna_target entries must be "&"-joined NTC gRNA ids. This is the validation path -- R's own pair selection is unseeded and varies run to run, so a pair-by-pair comparison is only meaningful when both sides are given the same pairs. ntc_grna_weights: {NTC gRNA id: dose weights}, aligned with ntc_grna_cells, to check the dose test instead of sceptre's: a synthetic target's cell takes its largest weight over the group's gRNAs (DoseWeights.ntc_grna_weights). Complement control group and CRT only, as in run_discovery_analysis. glm_solver: as in run_discovery_analysis; with control_group="nt_cells" anything but "outer" raises ValueError, as it does there.

Returns the same columns as run_discovery_analysis.

pysceptre.pipeline.api.run_power_check

run_power_check(response_matrix, gene_ids: list[str], covariate_matrix: ndarray, grna_target_cells: dict[str, ndarray], *, positive_control_pairs: DataFrame | None = None, n_nonzero_trt_thresh: int = 7, n_nonzero_cntrl_thresh: int = 7, side: str = 'both', moi: str = 'high', control_group: str | None = None, ntc_grna_cells: dict[str, ndarray] | None = None, resampling_approximation: str = 'skew_normal', multiple_testing_alpha: float = 0.1, seed: int | None = None, target_chunk_size: int = _DEFAULT_TARGET_CHUNK_SIZE, chunk_memory_gb: float = _DEFAULT_CHUNK_MEMORY_GB, n_jobs: int = 1, resampling_mechanism: str | None = None, grna_target_weights: dict[str, ndarray] | None = None, glm_solver: str = 'outer') -> DataFrame

Run sceptre's power check: the discovery test over positive controls.

A positive control is a (gene, target) pair where an effect is expected -- usually a gRNA against the gene's own TSS. The check asks whether the pipeline recovers effects it should, and is read next to run_calibration_check, which asks whether it invents effects it should not.

the pairs to test, as

DataFrame['response_id', 'grna_target']. Supply these. They are a biological claim about which target perturbs which gene, and only the experiment knows it. When omitted, R's name-matching rule is used as a fallback -- 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. A screen of the latter kind raises rather than silently returning an empty result.

n_nonzero_trt_thresh / n_nonzero_cntrl_thresh: pairwise QC thresholds. Unlike the calibration check, failures are reported, not filtered: the returned frame has a pass_qc column and NaN results for pairs that did not meet them. Dropping them would overstate power by hiding the controls the screen had too few cells to test. With control_group="nt_cells" the control count is the gene's nonzero NT cells, as in R. moi, control_group, ntc_grna_cells, resampling_mechanism, grna_target_weights, glm_solver: as in run_discovery_analysis, with the same MOI-dependent defaults.

Returns one row per supplied pair, with pass_qc, n_nonzero_trt and n_nonzero_cntrl alongside the usual columns. No multiple-testing correction is applied, matching R: these are a diagnostic, not discoveries, and adjusting them against each other answers no question.

Analytical power

A different question from the three above, and the one function here does not run a test at all. run_power_check is sceptre's positive-control diagnostic -- it runs the real test on pairs where an effect is expected. compute_power estimates in closed form what the screen could have detected, per pair, including for pairs the screen never tested.

Three arguments have no defaults on purpose

cutoff, fold_change_mean and fold_change_sd are all required, and cutoff must be the analysed screen's own nominal threshold. Borrowing another design's was measured at more than double the error with a systematically optimistic bias. See Design decisions.

pysceptre.analytical_power.estimate.compute_power

compute_power(discovery_pairs: DataFrame, cells_per_grna: DataFrame, baseline_expression_stats: DataFrame, *, fold_change_mean: float, fold_change_sd: float, cutoff: float, num_total_cells: int, side: str = 'left', n_nonzero_trt_thresh: int = 0, n_nonzero_cntrl_thresh: int = 0) -> DataFrame

Analytical power for each pair in discovery_pairs.

Parameters:

Name	Type	Description	Default
discovery_pairs	DataFrame	columns grna_target and response_id, the pairs to score. They need not be pairs the screen tested.	required
cells_per_grna	DataFrame	one row per individual gRNA, columns grna_target and num_cells (and conventionally grna_id). Per-gRNA granularity is required, not optional: the across-gRNA variance term needs the sum of squared counts, which no target-level total can supply.	required
baseline_expression_stats	DataFrame	one row per gene, columns response_id, expression_mean and expression_size. expression_size is the NB size -- theta, i.e. 1 / dispersion. With a grna_target column it is instead one row per pair, each matched to its pair by (grna_target, response_id). That is the shape inputs.matched_expression_stats returns. expression_mean is E[observed count in a cell], on the raw count scale, and the formula says so twice: var_nb is the variance of the negative binomial the counts themselves follow, and qc_failure_prob asks how many cells have a nonzero count. Only a distribution over actual counts has a zero probability. Build it with inputs.baseline_expression_stats_from_fits, which is on that scale by construction. A size-factor-normalised mean is not on this scale and is the mistake to avoid: inputs.baseline_expression_stats computes one, and on day0 it sits 16% low. Passing it understates power, and when either nonzero threshold is above 0 it does so twice -- through the statistic and again through the QC factor.	required
fold_change_mean	float	a multiplier, not an effect size: a 15% knockdown is 0.85.	required
fold_change_sd	float	across-gRNA sd of the fold change. Required, with no default; see docs/design.md, "Analytical per-pair power".	required
cutoff	float	the nominal p-value threshold this screen's test is read at. Required, and it must be this screen's threshold rather than a borrowed one, because it carries the multiple-testing correction; see docs/design.md, "Analytical per-pair power".	required
num_total_cells	int	cells in the analysis. The complement control group is num_total_cells - num_trt_cells.	required
side	str	"left" (the default, paired with cutoff = alpha / 2), "right", or "both" (paired with cutoff = alpha). The default is the configuration the estimator was validated in.	'left'
n_nonzero_trt_thresh	int	sceptre's pairwise nonzero-count threshold for the treated cells. Defaults to 0, which makes the QC factor 1 -- correct for pairs that already passed QC. Raise it to score pairs that were never tested.	0
n_nonzero_cntrl_thresh	int	the same, for the complement group.	0

Returns:

Type	Description
DataFrame	discovery_pairs with power appended, in the order supplied, plus
DataFrame	the inputs each row's estimate was built from (num_trt_cells,
DataFrame	num_trt_cells_sq, num_cntrl_cells, expression_mean,
DataFrame	expression_size, mean_test_stat, sd_test_stat, qc_failure_prob)
DataFrame	so a surprising number can be traced without recomputing it.

Raises:

Type	Description
KeyError	a pair names a target with no gRNA counts, or a gene (or, with per-pair statistics, a pair) with no baseline statistics. Never filled in silently: a missing input would otherwise become a NaN power indistinguishable from a real one.
ValueError	an argument is out of range, a key is duplicated, or the summed treated count leaves the complement group empty.

pysceptre.analytical_power.estimate.target_cell_counts

target_cell_counts(cells_per_grna: DataFrame) -> DataFrame

Aggregate per-gRNA cell counts to the per-target terms the formula needs.

cells_per_grna needs columns grna_target and num_cells, one row per individual gRNA. Rows whose target is "non-targeting" are dropped, as R does -- NTCs are not a target and the complement control group does not use them.

Returns grna_target, num_trt_cells (the sum over the target's gRNAs) and num_trt_cells_sq (the sum of their squares). The sum is not the size of the union of perturbed cells, and in high MOI the two differ. See docs/design.md, "Analytical per-pair power".

Building its inputs

Three things compute_power needs that no other part of pysceptre produces. They take arrays and frames rather than a loaded export, because nothing under src/ imports the io extra; the file-reading glue lives in scripts/.

pysceptre.analytical_power.inputs.cells_per_grna_from_assignments

cells_per_grna_from_assignments(grna_target_data_frame: DataFrame, targeting_grna_cells: dict[str, ndarray], *, drop_duplicate_design_rows: bool = False) -> DataFrame

Per-gRNA cell counts in the shape compute_power wants.

Parameters:

Name	Type	Description	Default
grna_target_data_frame	DataFrame	the screen's design, one row per (grna_id, grna_target). A guide inside several overlapping candidate elements appears once per target, and every one of those rows is kept.	required
targeting_grna_cells	dict[str, ndarray]	grna_id -> that guide's cell indices. Guides absent from it are guides with no cells.	required

Returns:

Type	Description
DataFrame	grna_id, grna_target, num_cells, one row per design row.

Two things this does that a shorter version would not, both of which change the answer rather than the style:

It iterates the design, not the guides that have cells. A designed guide that ended up with none is a num_cells = 0 row, not an absent one. Dropped instead, num_trt_cells_sq would be computed over a smaller guide set than the screen used, which is the term no target-level total can supply in the first place.

It keeps a guide under every target it belongs to. The map is many-to-many, so a guide in two overlapping elements contributes its cells to both sums; that is what R does. Deduplicating by grna_id looks like hygiene and silently shrinks exactly those targets.

A guide listed twice against the same target is a different matter and is refused, because it would be counted twice into num_trt_cells and its squared sum. Unlike the discovery path there is no parity argument for keeping it: PerturbPlan never reads a design table. drop_duplicate_design_rows=True collapses them instead of raising.

pysceptre.analytical_power.inputs.baseline_expression_stats_from_fits

baseline_expression_stats_from_fits(covariate_matrix: ndarray, gene_fits: Mapping[str, GeneFit], *, gene_subset: list[str] | None = None) -> DataFrame

Baseline statistics on the scale sceptre's own model works on.

This is the recommended way to build them. Both numbers come out of the same negative-binomial fit the discovery test itself uses, so the power estimate and the test it predicts are on one expression scale by construction rather than by coincidence.

Parameters:

Name	Type	Description	Default
covariate_matrix	ndarray	(n_cells, p), the same design matrix the fits were produced with. Using a different one silently answers a different question.	required
gene_fits	Mapping[str, GeneFit]	gene_id -> an object with .fitted_coefs and .theta, which is exactly what discovery.py::fit_all_genes returns. A completed analysis already has them.	required
gene_subset	list[str] | None	return only these genes, in this order.	None

Returns:

Type	Description
DataFrame	response_id, expression_mean, expression_size.

expression_mean is mean(exp(Z @ coefs)), the average expected count the fitted model gives the gene. That is the quantity the score statistic's variance is built from, so it is what the closed form wants. expression_size is the same fit's theta.

Why not the size-factor-normalised mean. baseline_expression_stats computes that instead, because it is what the published comparison used. Measured on day0 it sits about 16% below this one, near enough a constant across genes (correlation of logs 0.9999).

That near-constancy makes it tempting to call the difference a conservative rescaling, and through the test statistic it is. Through the QC factor it is not. compute_power multiplies power by 1 - qc_failure_prob, which is built from P(count == 0), and that is nonlinear in the mean: a 16% dim input overstates the zero probability by 0.032 at the median on day0 and 0.073 at worst, inflating the QC discount for 28 of 237 genes at a typical target, the worst by 0.20.

The QC channel is inert at the default thresholds of 0, which is the setting for pairs that already passed QC -- so the "conservative rescaling" reading holds there and only there. It stops holding in exactly the case the thresholds exist for: scoring pairs the screen never tested. See docs/design.md, "Analytical per-pair power".

pysceptre.analytical_power.inputs.matched_expression_stats

matched_expression_stats(discovery_pairs: DataFrame, covariate_matrix: ndarray, gene_fits: Mapping[str, GeneFit], grna_target_cells: dict[str, ndarray], cells_per_grna: DataFrame, num_total_cells: int) -> DataFrame

Per-pair baseline statistics that carry the covariates into compute_power.

An alternative to baseline_expression_stats_from_fits. Instead of one mean per gene, each (target, gene) pair gets the single mean that gives the closed form's statistic the same variance as sceptre's per-cell fit does for that target's cells. A cell's information is w_j = mu_j / (1 + mu_j / theta) with mu_j = exp(Z_j @ coefs), and with f_p and f_c its average over the target's perturbed cells and over the complement,

f = (1/n_p + 1/n_c) / (1/(n_p f_p) + 1/(n_c f_c))
expression_mean = f / (1 - f / theta)

The last line inverts w, and it cannot fail: every w_j < theta, so f does too.

Parameters:

Name	Type	Description	Default
discovery_pairs	DataFrame	columns grna_target and response_id, the pairs to score.	required
covariate_matrix	ndarray	(n_cells, p), the same design matrix the fits were produced with.	required
gene_fits	Mapping[str, GeneFit]	gene_id -> an object with .fitted_coefs and .theta, as discovery.py::fit_all_genes returns.	required
grna_target_cells	dict[str, ndarray]	target -> the indices, into the rows of covariate_matrix, of the cells carrying any of its gRNAs (the union), as run_discovery_analysis takes it. f_p and f_c are averaged over these cells and the rest.	required
cells_per_grna	DataFrame	one row per gRNA, as compute_power takes it. n_p is target_cell_counts(cells_per_grna)["num_trt_cells"], the per-gRNA sum, and n_c = num_total_cells - n_p.	required
num_total_cells	int	the value that will be passed to compute_power.	required

Returns:

Type	Description
DataFrame	grna_target, response_id, expression_mean, expression_size
DataFrame	and info_ratio_pert, one row per pair in the order supplied.
DataFrame	expression_size is the fit's theta. info_ratio_pert is f_p over
DataFrame	the gene's all-cell mean of w, kept so a surprising pair can be
DataFrame	traced. The frame can be passed to compute_power as its
DataFrame	baseline_expression_stats.

The averages use the union and the counts use the sum on purpose; do not make them agree. See docs/design.md, "Covariates through the information-matched mean".

Raises:

Type	Description
KeyError	a pair names a gene with no fit, a target with no cells in grna_target_cells, or a target with no rows in cells_per_grna.
ValueError	the design and a fit disagree, a target has no cells or all of them, or a result is not finite and positive.

pysceptre.analytical_power.inputs.baseline_expression_stats

baseline_expression_stats(counts: spmatrix, gene_ids: list[str], thetas: dict[str, float] | ndarray, *, size_factors: ndarray | None = None, gene_subset: list[str] | None = None) -> DataFrame

Per-gene expression_mean and expression_size for compute_power.

Parameters:

Name	Type	Description	Default
counts	spmatrix	(n_genes, n_cells) sparse counts for all genes and cells.	required
gene_ids	list[str]	row labels of counts, in order.	required
thetas	dict[str, float] | ndarray	per-gene NB size. A mapping keyed by gene id, or an array aligned to gene_ids. discovery.py::fit_all_genes returns objects carrying exactly this as .theta, so a completed analysis already has them and nothing needs refitting.	required
size_factors	ndarray | None	precomputed per-cell factors, to skip a pass over the matrix. Computed by poscounts_size_factors when omitted.	None
gene_subset	list[str] | None	return only these genes. The statistics are still computed over the full matrix first, which is the point: the normalisation depends on the whole gene and cell set, so subsetting the input and subsetting the output are different operations with different answers.	None

Returns:

Type	Description
DataFrame	response_id, expression_mean, expression_size.

expression_mean is the size-factor-normalised mean: each cell's counts divided by its size factor, then averaged over cells. expression_size is the NB size, theta, not the dispersion.

This is the wrong scale for compute_power, and is kept only to reproduce the published comparison. Use baseline_expression_stats_from_fits. The mean here comes from a normalisation scheme sceptre does not use, and on day0 it sits about 16% below the mean sceptre's own model implies -- which is what compute_power needs, since its var_nb is the variance of the counts themselves and its QC factor asks how many cells hold a nonzero count. Calling this function emits a warning for that reason.

Two caveats on theta, because a number that transfers badly here is invisible in the output. fit_all_genes clamps theta to (0.01, 1000.0), so a gene at a bound carries the bound rather than its estimate; and theta is fitted under whatever covariate_matrix was passed, so it matches another implementation's only if the design matrices match.

pysceptre.analytical_power.inputs.poscounts_size_factors

poscounts_size_factors(counts: spmatrix) -> ndarray

DESeq2 "poscounts" size factors, one per cell, over the nonzeros only.

Parameters:

Name	Type	Description	Default
counts	spmatrix	(n_genes, n_cells) sparse counts. All genes and all cells, not a subset -- see the note below.	required

Returns:

Type	Description
ndarray	(n_cells,) positive size factors.

Reproduces the R this estimator's expression_mean was validated against (WattEG prepare_sim_input.R:157-223): a per-gene geometric mean over the nonzero counts, then each cell's factor is the median log-ratio of its own nonzeros against those means.

The gene and cell sets are part of the definition. The geometric mean runs over every gene and the median over every cell, so computing this on a subset gives different factors for the cells that remain, and therefore different normalised means for every gene. Both results look equally complete. Pass the whole matrix.

DESeq2 itself is not used because DESeqDataSetFromMatrix densifies; this walks the nonzeros.

Raises:

Type	Description
ValueError	a negative count, or a cell whose factor is undefined because it has no nonzero in any gene with a finite geometric mean.

pysceptre.analytical_power.inputs.bh_nominal_cutoff

bh_nominal_cutoff(p_values: ndarray, alpha: float) -> float

The nominal p-value threshold a BH-corrected analysis actually applied.

Parameters:

Name	Type	Description	Default
p_values	ndarray	the discovery p-values. NaN entries -- the pairs that failed pairwise QC and were never tested -- are dropped before the correction, as R does, rather than inflating its denominator.	required
alpha	float	the analysis's multiple_testing_alpha.	required

Returns:

Type	Description
float	The largest p-value BH calls significant, which is what
float	compute_power's cutoff wants (halved, for the default
float	side="left").

Ports WattEG's discovery_threshold(), including its refusal to return a number when nothing is significant. The code that preceded it returned -Inf there and every pair silently got zero power, so an empty result is raised rather than propagated.

pysceptre applies no multiple-testing correction itself, so this does both halves: the BH step and then reading the threshold off it.

Raises:

Type	Description
ValueError	no usable p-value, or none significant at alpha.

Specificity check

How many discovered links are more than background, measured on tests between elements and genes on other chromosomes. Not from sceptre; see Design decisions.

pysceptre.specificity.check.run_specificity_check

run_specificity_check(cis_result: DataFrame, trans_result: DataFrame, element_positions: DataFrame, gene_positions: DataFrame, cutoff: float, *, control_targets: Collection[str] = (), measured_genes: Collection[str] | None = None, distance_bins: Sequence[float] = DEFAULT_DISTANCE_BINS, far_distance: float = 100000.0, tss_window: float = 1000.0, broad_effect_p: float = 0.001, nominal: float = 0.05, exclude_gene_chroms: Collection[str] = ('chrY',), n_bootstrap: int = 2000, seed: int | Generator | None = None) -> SpecificityResult

Which discovered links are more than background.

The calibration check asks whether the pipeline invents effects and the power check whether it recovers effects it should. This asks how many of the links it did discover are more than the rate at which a test is called without regulation, measured on tests between elements and genes on other chromosomes. It is not sceptre's and has no R counterpart; see docs/design.md, "Specificity check", for its scope and evidence.

Parameters:

Name	Type	Description	Default
cis_result	DataFrame	the screen's discovery result on element-gene pairs on the same chromosome: grna_target, response_id, p_value, log_2_fold_change (or fold_change, from which it is derived), and pass_qc if it has one. Ids may be categorical. A screened result is accepted as for trans_result.	required
trans_result	DataFrame	a discovery result, same columns and the same test, for the screen's elements against its cis genes on other chromosomes. Pairs outside those sets are ignored. A screened result from run_trans_analysis is accepted: a row with screen "fold_change" or "early_stop" was tested and not resampled, so it counts as a test that is not a call.	required
element_positions	DataFrame	grna_target, chrom, centre, for every target in cis_result. Take these from the screen's own coordinate table, never from parsing target names.	required
gene_positions	DataFrame	response_id, chrom, tss. May be a whole annotation; the genes the check uses must each appear once.	required
cutoff	float	the nominal p-value threshold the cis analysis was called at, e.g. bh_nominal_cutoff(p_values, alpha).	required
control_targets	Collection[str]	targets that are not elements (TSS positive controls). Left out of the link and background counts, and still used as TSS knockdowns.	()
measured_genes	Collection[str] | None	the genes the screen measured, for finding TSS targets. Defaults to every gene in either result.	None
distance_bins	Sequence[float]	bin edges in bp for the link counts.	DEFAULT_DISTANCE_BINS
far_distance	float	links beyond this are the ones the detour check reads.	100000.0
tss_window	float	a target within this many bp of a measured TSS knocks that gene down.	1000.0
broad_effect_p	float	binomial threshold for the broad-effect flag.	0.001
nominal	float	the p-value for the nominal detour level.	0.05
exclude_gene_chroms	Collection[str]	gene chromosomes left out of the background.	('chrY',)
n_bootstrap	int	bootstrap draws over elements.	2000
seed	int | Generator | None	an int, or a Generator to continue drawing from. The draws with all elements come first, then the draws without the broad-effect elements, whether or not any element is flagged.	None

Returns:

Type	Description
SpecificityResult	A SpecificityResult.

Raises:

Type	Description
KeyError	a column, a target position or a gene position is missing.
ValueError	cutoff is out of range, a cis pair spans two chromosomes, a gene has two positions, or a screened row could be a call at cutoff (or, for a TSS target, at nominal).

Warns:

Type	Description
UserWarning	some cis genes have no background test. Their background rate counts as 0 in the matched background.

pysceptre.specificity.check.SpecificityResult dataclass

Everything run_specificity_check computes, as tables.

Attributes:

Name	Type	Description
by_distance	DataFrame	per distance bin, the cis rate, the matched background and links above background, with bootstrap intervals.
by_distance_without_broad	DataFrame	the same with the broad-effect elements removed from both the cis and the background tests.
broad_effect	DataFrame	per element, its background calls and the binomial test against the screen's background rate.
per_gene	DataFrame	per gene, its background tests, calls and rate.
links	DataFrame	the tested element-gene pairs, with distance and link call.
background	DataFrame	the tests the background is measured on.
tss_targets	DataFrame	targets within the TSS window of a measured gene.
gene_lookup	DataFrame	what each TSS knockdown lowers.
far_links	DataFrame	the detour check of every far link.
background_rate	float	the pooled background rate.

Its background

The trans pairs the check reads, tested with the cis run's own test. Not from sceptre; see Design decisions.

pysceptre.pipeline.api.run_trans_analysis

run_trans_analysis(response_matrix, gene_ids: list[str], covariate_matrix: ndarray, grna_target_cells: dict[str, ndarray], pairs: DataFrame, *, reference_pairs: DataFrame, seed: int, side: str = 'both', resampling_approximation: str = 'skew_normal', target_chunk_size: int = _DEFAULT_TARGET_CHUNK_SIZE, chunk_memory_gb: float = _DEFAULT_CHUNK_MEMORY_GB, expected_fit_chunks: tuple[int, int] | None = None, fit_chunk_ids=None, gene_fits=None, p_exact_below: float = 0.05, block_rows: int = 32, n_jobs: int = 1, backend: str | None = None, fit_concurrency: int = 3, glm_solver: str = 'outer') -> DataFrame

Test background pairs with a reference discovery run's test, for run_specificity_check.

Every pair in pairs is tested as the reference run, which was given reference_pairs and the same data, seed and settings, would test it: targets are fit in the reference run's own chunks, so a pair tested in full gets the result run_discovery_analysis gives it, bit for bit on the same CPU class and BLAS, with numba, and with the covariate matrix in the same memory layout. Not from sceptre; see docs/design.md, "The trans background".

Two screens skip work the specificity check cannot see, and are exact for every table it returns. A pair with fold_change >= 1 is not resampled (screen "fold_change"), and the first resampling stage stops once the p-value is known to be at least p_exact_below and the pair cannot escalate (screen "early_stop", with that bound in p_value_floor). A screened pair has p_value NaN, which here means "tested, not resampled". A mode column in pairs selects per pair "screened" (the default), "full" (no screen) or "audit" (screened, with the full result beside it in the audit_* columns).

Complement control group, CRT, skew_normal, union strategy only; seed is required, and so is numba. expected_fit_chunks is the reference layout's (width, count), refused when different; fit_chunk_ids selects reference chunks for a run split into shards. gene_fits reuses width-1 gene fits (see pysceptre.pipeline.gene_cache): a GeneFitCache is checked against the data as load_gene_cache checks it (the covariate matrix's values and memory layout, the cell count and the counts of every gene it supplies) and must have been made with this glm_solver; a plain dict of fits is not checked. glm_solver must be the reference run's (see run_discovery_analysis).

block_rows is the number of first-stage resamples evaluated between early-stop checks. It decides where a pair can stop early, so whether it is screened "early_stop" and its p_value_floor and n_null_evaluated; a pair that is not screened gets the full test's result whatever its value. n_jobs sets the workers for the per-gene tests and backend their kind, "fork" (processes) or "thread", the platform default when None; fit_concurrency is the number of target chunks fit at once. None of the three can change a result. iter_trans_analysis yields the same results one gene at a time, for runs too large to hold as one table.

Returns:

Type	Description
DataFrame	response_id, grna_target, p_value, fold_change,
DataFrame	se_fold_change, z_orig, stage, screen, p_value_floor,
DataFrame	n_null_evaluated, and the audit_* columns when a pair was audited;
DataFrame	one row per pair, in pairs order.

Raises:

Type	Description
ValueError	a setting that cannot reproduce the reference run, a pair outside the reference chunks, an unknown gene or mode, or a GeneFitCache that does not match the data or glm_solver.
RuntimeError	numba is not available; the trans background needs it.

Gene fits made once and reused by every shard of a trans run, passed as gene_fits. A cache records what its fits and the pieces built from them depend on, the covariate matrix's memory layout included, and is refused where that differs.

pysceptre.pipeline.gene_cache.fit_gene_cache

fit_gene_cache(response_matrix, gene_ids: list[str], covariate_matrix: ndarray, genes=None, *, n_jobs: int = 1, glm_solver: str = 'outer') -> GeneFitCache

Fit genes exactly as the complement engine fits them, with a fingerprint to check later.

The covariate matrix is validated and converted as run_discovery_analysis does, and the fits are the engine's own fit_all_genes call: one gene per fit and the engine's gene memory budget, with the shared x_outer_flat under the default glm_solver="outer" and the structured design under "structured". A gene's fit depends only on its counts and the covariates, so fitting any subset gives each gene the fit a run over all of them gives it. The fits also depend on the covariate matrix's memory layout, which the fingerprint records; the matrix is used as given.

Parameters:

Name	Type	Description	Default
response_matrix		(n_genes, n_cells) counts, dense or sparse, as for run_discovery_analysis.	required
gene_ids	list[str]	the row labels of response_matrix, in order.	required
covariate_matrix	ndarray	(n_cells, p) design matrix, in the memory layout the fits will be used with.	required
genes		the genes to fit; None fits every row. Fitted in gene_ids order.	None
n_jobs	int	workers for the fits; cannot change a fit.	1
glm_solver	str	the engine's glm_solver, recorded in the fingerprint; fits of one solver are never loaded for the other.	'outer'

Returns:

Type	Description
GeneFitCache	A GeneFitCache with one fit per gene, keyed by its id.

Raises:

Type	Description
ValueError	a covariate matrix run_discovery_analysis refuses, gene_ids not labelling every row, or genes naming an id that gene_ids lacks or holds twice.

pysceptre.pipeline.gene_cache.save_gene_cache

save_gene_cache(cache: GeneFitCache, path) -> None

Write cache to path as an .npz that np.load(..., allow_pickle=False) reads.

Arrays: gene_ids (unicode), coefs (G, p) float64, theta, min_eigenvalue and max_eigenvalue float64, theta_method and n_covariates int64, theta_clamped and glm_converged bool, and fingerprint, the fingerprint as a JSON string. Floats are stored as their float64 bits. The file is written to path exactly, without a suffix added.

Raises:

Type	Description
ValueError	a gene id that is not a string, or coefficients whose length is not the fingerprint's n_covariates.

pysceptre.pipeline.gene_cache.load_gene_cache

load_gene_cache(paths, *, response_matrix=None, gene_ids=None, covariate_matrix=None, check_environment: bool = True, glm_solver: str | None = None) -> GeneFitCache

Read gene fits written by save_gene_cache, refusing any that cannot be reused.

Several files are merged into one cache; they must hold disjoint genes fit on the same data with the same settings. Every file must carry the engine's current gene batch width (1) and theta bounds. Given the data, the covariate matrix must hash to the fingerprint's and have its memory layout, and every cached gene's counts must match its digest. With check_environment, every file's numeric_environment_key must equal this process's. A file whose fingerprint does not record the covariate layout is refused, to be refit. See docs/design.md, "The trans background".

Parameters:

Name	Type	Description	Default
paths		one path, or a sequence of paths.	required
response_matrix, gene_ids		the data the fits will be used with; checked together, so give both or neither.	required
covariate_matrix		the covariate matrix the fits will be used with, in the memory layout they will be used with.	None
check_environment	bool	refuse fits made in another numeric environment.	True
glm_solver	str | None	refuse fits made with another glm_solver; None accepts either, provided all files agree.	None

Returns:

Type	Description
GeneFitCache	A GeneFitCache of every file's fits, in file order. Its fingerprint
GeneFitCache	is the files' shared one, with the count digests merged and the
GeneFitCache	environment of the first file. Floats are bitwise those saved.

Raises:

Type	Description
ValueError	no path; a file save_gene_cache did not write, or one that does not record the covariate layout; another batch width, theta bounds or glm_solver; files fit on different data or settings, or sharing a gene; a covariate matrix (values or memory layout), cell count, gene or counts other than the fits were made with; or another numeric environment. The message names what differs.

pysceptre.pipeline.gene_cache.GeneFitCache dataclass

Gene fits and the fingerprint they were made under.

fits maps a gene id to its GenePrecomputation, as fit_all_genes returns it. fingerprint holds format; n_cells, n_covariates and dfr of the fit; covariate_sha256, of the covariate matrix as float64 C-order bytes; covariate_layout, that matrix's memory layout ("C" or "F"), which the outer route's fits and every route's per-pair pieces depend on in their last bits; theta_bounds, batch_width and glm_solver, the engine settings the fits were made with; count_digests, per gene the nnz, sum and sha1 of its nonzero counts; and environment, an environment_fingerprint().

Its pieces

pysceptre.specificity.background.cis_links

cis_links(cis_result: DataFrame, element_positions: DataFrame, gene_positions: DataFrame, cutoff: float, *, control_targets: Collection[str] = ()) -> DataFrame

The tested element-gene pairs, with their distance and whether each is a link.

Parameters:

Name	Type	Description	Default
cis_result	DataFrame	a discovery result (grna_target, response_id, p_value, log_2_fold_change or fold_change, optionally pass_qc) on pairs with the element and the gene on the same chromosome. May be a screened result (screen, p_value_floor): its screened rows are tests that are not links.	required
element_positions	DataFrame	grna_target, chrom, centre.	required
gene_positions	DataFrame	response_id, chrom, tss.	required
cutoff	float	the nominal p-value threshold the screen was called at.	required
control_targets	Collection[str]	targets that are not elements, such as TSS positive controls. Their pairs are left out.	()

Returns:

Type	Description
DataFrame	grna_target, response_id, p_value, log_2_fold_change, then
DataFrame	screen and p_value_floor when cis_result has them, distance
DataFrame	(element centre to TSS) and link (called with a negative effect),
DataFrame	one row per tested pair.

Raises:

Type	Description
KeyError	a target or gene has no position.
ValueError	a pair has its element and gene on different chromosomes, or a screened row of an element with a negative effect could be a link at cutoff: screened on its fold change, or stopped early with a p_value_floor that is missing or not above cutoff.

pysceptre.specificity.background.background_pairs

background_pairs(trans_result: DataFrame, links: DataFrame, element_positions: DataFrame, gene_positions: DataFrame, cutoff: float, *, exclude_gene_chroms: Collection[str] = ('chrY',)) -> DataFrame

The tests the background is measured on, and which of them are called.

The cis elements against the cis genes, gene on another chromosome, gene chromosome not in exclude_gene_chroms, tested. Restricting to the cis sets is what lets the two sides be matched gene for gene. trans_result may be a screened trans result (run_trans_analysis): its screened rows are tests that are not calls.

Returns:

Type	Description
DataFrame	grna_target, response_id, p_value, log_2_fold_change, then
DataFrame	screen and p_value_floor when trans_result has them, and
DataFrame	background (called at the cis cutoff with a negative effect).

Raises:

Type	Description
ValueError	a screened row with a negative effect could be a call at cutoff: screened on its fold change, or stopped early with a p_value_floor that is missing or not above cutoff.

pysceptre.specificity.background.above_background_by_distance

above_background_by_distance(links: DataFrame, background: DataFrame, *, distance_bins: Sequence[float] = DEFAULT_DISTANCE_BINS, n_bootstrap: int = 2000, seed: int | Generator | None = None) -> DataFrame

Per distance bin: cis rate, matched background and links above background.

The background of a bin is the mean, over its tested pairs, of each pair's gene's background rate, so a gene counts as often as it is tested there. Links above background is (cis rate - background) / cis rate. Intervals are a bootstrap over elements: an element's cis and background tests are drawn together, and every gene's rate is recomputed in each draw.

Parameters:

Name	Type	Description	Default
links	DataFrame	cis_links output.	required
background	DataFrame	background_pairs output for the same elements.	required
distance_bins	Sequence[float]	bin edges in bp, increasing. Bins are closed on the right, and the first also includes its left edge.	DEFAULT_DISTANCE_BINS
n_bootstrap	int	bootstrap draws.	2000
seed	int | Generator | None	an int, or a Generator to continue drawing from.	None

Returns:

Type	Description
DataFrame	One row per bin: distance_bin, pairs, links, then cis_rate,
DataFrame	background and above, each with _lo and _hi (the 2.5 and 97.5
DataFrame	percentiles of the draws).

Raises:

Type	Description
ValueError	a distance falls outside the bins, or the bins do not increase.

pysceptre.specificity.background.broad_effect_elements

broad_effect_elements(background: DataFrame, element_positions: DataFrame, gene_positions: DataFrame, *, p_threshold: float = 0.001) -> DataFrame

Elements with more background calls than the screen's background rate allows.

Each element's background calls are tested against the pooled rate over all background tests (one-sided binomial). Such an element makes its own background, which a per-gene rate cannot see.

Returns:

Type	Description
DataFrame	grna_target, trans_pairs, trans_calls, p_excess,
DataFrame	broad_effect (p_excess < p_threshold), nearest_gene and
DataFrame	tss_distance (the nearest TSS in gene_positions), one row per
DataFrame	element with background tests.

pysceptre.specificity.background.nearest_tss

nearest_tss(chrom: Series, position: Series, gene_positions: DataFrame) -> DataFrame

The nearest TSS in gene_positions on the same chromosome, for each position.

Returns nearest_gene and tss_distance, indexed like chrom; NaN on a chromosome with no gene. A tie goes to the upstream TSS.

pysceptre.specificity.lookup.tss_targets

tss_targets(targets: Collection[str], element_positions: DataFrame, gene_positions: DataFrame, measured_genes: Collection[str], *, window: float = 1000.0) -> DataFrame

Targets within window bp of a measured gene's TSS, and that gene.

A target near several measured TSSs is assigned the nearest.

Returns:

Type	Description
DataFrame	grna_target, tss_gene, tss_distance.

Raises:

Type	Description
KeyError	a target has no position.

pysceptre.specificity.lookup.gene_lookup

gene_lookup(tss: DataFrame, cis_result: DataFrame, trans_result: DataFrame, cutoff: float, *, nominal: float = 0.05) -> DataFrame

What each TSS knockdown lowers: gene G knocked down -> gene H.

Every tested pair of a TSS target, from either result, against any gene other than its own. Several TSS targets of one gene are pooled. A screened row of either result counts as tested, not lowering.

Returns:

Type	Description
DataFrame	knocked_down, gene, tested, lowers (called at cutoff with a
DataFrame	negative effect by any of G's TSS targets), lowers_nominal (the same
DataFrame	at p < nominal) and best_lfc (the most negative log2 fold change).

Raises:

Type	Description
ValueError	a screened row of a TSS target with a negative effect could lower at cutoff or at nominal: screened on its fold change, or stopped early with a p_value_floor that is missing, not above cutoff or below nominal.

pysceptre.specificity.lookup.detour_check

detour_check(links: DataFrame, tss: DataFrame, lookup: DataFrame, *, far_distance: float = 100000.0) -> DataFrame

For each far link E -> H, whether it can run through another gene E lowers.

The genes E lowers are its links at any distance, plus the gene whose TSS E sits on (knocked down even when that pair was not tested).

Returns:

Type	Description
DataFrame	The far links (distance > far_distance) with element_lowers and a
DataFrame	status: "detour" when some G that E lowers has a TSS knockdown
DataFrame	that lowers H at the cutoff, "detour (nominal)" when only at the
DataFrame	nominal level, "no detour" when such G were tested against H and
DataFrame	none lowers it, and "can't check" when no G was tested against H.
DataFrame	through names the G that explain a detour.

gRNA assignment

Which cells carry which gRNA, from raw gRNA UMI counts. assign_grnas is the one entry point, like sceptre's, and runs any of the four methods in pysceptre.assignment.ASSIGNMENT_METHODS. assign_grnas_mixture, assign_grnas_thresholding and assign_grnas_maximum are sceptre's three methods; assign_grnas_fishash is not sceptre's (a port of fishash). See Design decisions.

pysceptre.assignment.api.assign_grnas

assign_grnas(grna_matrix: spmatrix | ndarray, grna_ids: Sequence[str], *, method: str = 'default', moi: str | None = None, covariate_matrix: ndarray | None = None, **hyperparameters: Any) -> MixtureResult | ThresholdingResult | MaximumResult | FishashResult

Assign gRNAs to cells with one of the package's methods.

method="default" is sceptre's default: "maximum" for a low-MOI screen and "mixture" for a high-MOI one, so it needs moi. The methods are "mixture" (assign_grnas_mixture, which needs covariate_matrix), "thresholding" (assign_grnas_thresholding), "maximum" (assign_grnas_maximum, refused for a high-MOI screen as in sceptre) and "fishash" (assign_grnas_fishash, not from sceptre). Each method's options are passed through hyperparameters under that function's keyword names.

Parameters:

Name	Type	Description	Default
grna_matrix	spmatrix | ndarray	(n_grnas, n_cells) raw integer gRNA UMI counts, scipy.sparse or dense.	required
grna_ids	Sequence[str]	The n_grnas row names, distinct.	required
method	str	"default" or one of ASSIGNMENT_METHODS.	'default'
moi	str | None	"low" or "high"; required by "default".	None
covariate_matrix	ndarray | None	The mixture's (n_cells, p) design; mixture_design_matrix builds sceptre's default.	None
**hyperparameters	Any	Keyword options of the chosen method's function.	{}

Returns:

Type	Description
MixtureResult | ThresholdingResult | MaximumResult | FishashResult	The chosen method's result: MixtureResult, ThresholdingResult, MaximumResult or
MixtureResult | ThresholdingResult | MaximumResult | FishashResult	FishashResult.

Raises:

Type	Description
ValueError	an unknown method or moi, "default" without moi, "maximum" for a high-MOI screen, "mixture" without covariate_matrix, or covariate_matrix given to another method.
TypeError	a hyperparameter the chosen method does not take.

pysceptre.assignment.mixture.assign_grnas_mixture

assign_grnas_mixture(grna_matrix: spmatrix | ndarray, grna_ids: Sequence[str], covariate_matrix: ndarray, *, probability_threshold: float = 0.8, n_nonzero_cells_cutoff: int = 10, backup_threshold: float = 5, pi_guesses: Sequence[float] = PI_GUESSES, g_pert_guesses: Sequence[float] = G_PERT_GUESSES) -> MixtureResult

Assign gRNAs to cells with sceptre's mixture method.

Parameters:

Name	Type	Description	Default
grna_matrix	spmatrix | ndarray	(n_grnas, n_cells) raw integer gRNA UMI counts, scipy.sparse or dense.	required
grna_ids	Sequence[str]	The n_grnas row names, distinct.	required
covariate_matrix	ndarray	(n_cells, p) design, intercept included -- sceptre's default is built by mixture_design_matrix.	required
probability_threshold	float	Posterior a cell needs to be assigned.	0.8
n_nonzero_cells_cutoff	int	A gRNA with fewer cells of count >= 1 uses the backup rule.	10
backup_threshold	float	The backup rule assigns cells whose count reaches this.	5
pi_guesses	Sequence[float]	Starting mixing proportions, one per EM start.	PI_GUESSES
g_pert_guesses	Sequence[float]	Starting perturbation effects, one per start.	G_PERT_GUESSES

Returns:

Type	Description
MixtureResult	A MixtureResult.

Raises:

Type	Description
ValueError	invalid inputs, or a GLM fit with a non-finite fitted mean (the gRNA is named; sceptre's runners retry with a smaller design then).

pysceptre.assignment.mixture.MixtureResult dataclass

What assign_grnas_mixture returns.

Attributes:

Name	Type	Description
assigned	csr_matrix	(n_grnas, n_cells) boolean CSR matrix of assignments.
posterior	csr_matrix	(n_grnas, n_cells) CSR matrix of the chosen start's posterior for the second component, at each mixture-path gRNA's cells with a nonzero count (and any zero-count cell with a fitted mean above 200).
posterior_zero	ndarray	(n_grnas,) the posterior shared by a gRNA's other zero-count cells; NaN for gRNAs on the backup rule.
fits	DataFrame	One row per gRNA: grna_id, n_nonzero, method ("mixture" or "backup"), backup_reason ("", "n_nonzero" or "em"), em_converged, log_lik, best_start (0-based), pi, g_pert (the chosen start's last estimates), glm_converged, glm_n_iter, n_assigned.
grna_ids	tuple[str, ...]	Row names.

pysceptre.assignment.thresholding.assign_grnas_thresholding

assign_grnas_thresholding(grna_matrix: spmatrix | ndarray, grna_ids: Sequence[str], *, threshold: float = 5) -> ThresholdingResult

Assign gRNAs to cells with sceptre's thresholding method.

A gRNA is assigned to every cell in which its UMI count is at least threshold.

Parameters:

Name	Type	Description	Default
grna_matrix	spmatrix | ndarray	(n_grnas, n_cells) raw integer gRNA UMI counts, scipy.sparse or dense.	required
grna_ids	Sequence[str]	The n_grnas row names, distinct.	required
threshold	float	The count a gRNA needs in a cell; at least 1, as sceptre requires.	5

Returns:

Type	Description
ThresholdingResult	A ThresholdingResult.

Raises:

Type	Description
ValueError	invalid counts or ids, or a threshold that is not a finite number >= 1.

pysceptre.assignment.thresholding.ThresholdingResult dataclass

What assign_grnas_thresholding returns.

Attributes:

Name	Type	Description
assigned	csr_matrix	(n_grnas, n_cells) boolean CSR matrix of assignments.
threshold	float	The count a gRNA needed in a cell.
grna_ids	tuple[str, ...]	Row names.

pysceptre.assignment.maximum.assign_grnas_maximum

assign_grnas_maximum(grna_matrix: spmatrix | ndarray, grna_ids: Sequence[str], *, umi_fraction_threshold: float = 0.8, min_grna_n_umis_threshold: float = 5) -> MaximumResult

Assign gRNAs to cells with sceptre's maximum method, which sceptre allows in low MOI only.

Every cell is assigned the gRNA with the most UMIs in it, the first in row order on a tie. A cell is flagged in cells_w_zero_or_twoplus_grnas when that gRNA holds at most umi_fraction_threshold of the cell's gRNA UMIs, or when the cell has fewer than min_grna_n_umis_threshold gRNA UMIs. Drop the flagged cells before a low-MOI analysis, as sceptre's QC does.

Parameters:

Name	Type	Description	Default
grna_matrix	spmatrix | ndarray	(n_grnas, n_cells) raw integer gRNA UMI counts, scipy.sparse or dense.	required
grna_ids	Sequence[str]	The n_grnas row names, distinct.	required
umi_fraction_threshold	float	In (0, 1), as sceptre requires.	0.8
min_grna_n_umis_threshold	float	At least 0, as sceptre requires.	5

Returns:

Type	Description
MaximumResult	A MaximumResult.

Raises:

Type	Description
ValueError	invalid counts or ids, a matrix with no gRNA, umi_fraction_threshold outside (0, 1), or a negative or non-finite min_grna_n_umis_threshold.

pysceptre.assignment.maximum.MaximumResult dataclass

What assign_grnas_maximum returns.

Attributes:

Name	Type	Description
assigned	csr_matrix	(n_grnas, n_cells) boolean CSR matrix with exactly one assignment per cell, its top gRNA. A cell with no gRNA UMIs is assigned the first gRNA, as in sceptre; the UMI rule flags it.
cells_w_zero_or_twoplus_grnas	ndarray	Ascending 0-based int64 indices of the cells sceptre flags, and its low-MOI QC removes: the top gRNA holds at most umi_fraction_threshold of the cell's gRNA UMIs, or the cell has fewer than min_grna_n_umis_threshold gRNA UMIs.
max_grna	ndarray	(n_cells,) int64 row of each cell's top gRNA, the first on a tie.
max_grna_frac_umis	ndarray	(n_cells,) the top gRNA's share of the cell's gRNA UMIs; NaN for a cell with none.
grna_n_umis	ndarray	(n_cells,) the cell's gRNA UMIs.
grna_ids	tuple[str, ...]	Row names.

pysceptre.assignment.fishash.assign_grnas_fishash

assign_grnas_fishash(grna_matrix: spmatrix | ndarray, grna_ids: Sequence[str], *, padj_cutoff: float = 0.05, padj_method: str = 'GS', min_count: float = 2, min_frac: float = 0.0, refit: int = 10, exclude_empty: bool = True, background: spmatrix | ndarray | None = None) -> FishashResult

Assign guides to cells with fishash.

Port of fishash's fishash() with the same defaults. For each nonzero count it runs a one-sided Fisher exact test of whether guide and cell co-occur more than their totals predict, and assigns the entries that pass an FDR cut and the count and fraction floors. With refit > 0, each later pass masks the previous pass's assignments, refits the noise as a rank-one Poisson matrix (impute_masked_counts), and takes the off-cell margins from that fit; the loop stops early once the assignments stop changing.

Parameters:

Name	Type	Description	Default
grna_matrix	spmatrix | ndarray	(n_grnas, n_cells) raw integer UMI counts, scipy.sparse or dense.	required
grna_ids	Sequence[str]	The n_grnas row names, distinct.	required
padj_cutoff	float	FDR level, in (0, 1).	0.05
padj_method	str	"GS" (Guo and Sarkar's procedure, which treats tests within a cell as dependent and cells as independent), "BH" or "BY".	'GS'
min_count	float	An entry needs at least this many UMIs to be assigned (when positive).	2
min_frac	float	An entry needs at least this fraction of its cell's UMIs (when positive).	0.0
refit	int	Maximum number of refit passes after the first.	10
exclude_empty	bool	Count only guides and cells with a nonzero total in the number of tests.	True
background	spmatrix | ndarray | None	Optional (n_grnas, n_cells) matrix of noise counts used for the off-cell margins of the first pass; only allowed with refit = 0.	None

Returns:

Type	Description
FishashResult	A FishashResult.

Raises:

Type	Description
ValueError	invalid arguments, a matrix with no nonzero count, a background given with refit > 0 (R refuses it too), or a NaN p-value (where R would return NA).
RuntimeError	fishash's own BH/BY consistency check fails.

pysceptre.assignment.fishash.FishashResult dataclass

What assign_grnas_fishash returns.

Attributes:

Name	Type	Description
assigned	csr_matrix	(n_grnas, n_cells) boolean CSR matrix, True where a guide is assigned.
log_pval	csr_matrix	(n_grnas, n_cells) CSR matrix of log p-values at every nonzero count.
odds_ratio	csr_matrix	Empirical odds ratio of each tested 2x2 table, same pattern.
odds_ratio_regularized	csr_matrix	The odds ratio with 1 added to both off-diagonal counts.
demux_type	ndarray	(n_cells,) "singlet", "doublet" or "unknown" (no guide assigned).
assignment	ndarray	(n_cells,) comma-joined ids of the assigned guides, "" if none, in fishash's order (gRNA order under GS; least to most significant under BH and BY).
log_pval_cutoff	float	The final pass's cut on the log p-value.
num_iter	int	Number of passes run.
grna_ids	tuple[str, ...]	Row names.
iterations	tuple[FishashIteration, ...]	One FishashIteration per pass.
background_guide_freqs	ndarray | None	Guide frequencies of the noise fit the final pass used, or None when no refit pass ran.
background_cell_sizes	ndarray | None	Cell sizes of that fit, or None.

Its pieces

pysceptre.assignment.design.mixture_design_matrix

mixture_design_matrix(grna_matrix: spmatrix | ndarray, *, response_n_nonzero: ndarray | None = None, response_n_umis: ndarray | None = None, extra_covariates: DataFrame | None = None) -> tuple[ndarray, list[str]]

The design sceptre's mixture assignment uses by default, and its column names.

The covariates, in sceptre's order, are response_n_nonzero and response_n_umis (when given), grna_n_nonzero and grna_n_umis (computed from grna_matrix), then the columns of extra_covariates. A covariate whose name contains n_umis or n_nonzero enters as log(x), or log(x + 1) when any cell has zero; any other covariate with 15 or more distinct values is left out; the rest enter as they are, categorical ones as treatment dummies with the first level as reference. An intercept comes first.

Parameters:

Name	Type	Description	Default
grna_matrix	spmatrix | ndarray	(n_grnas, n_cells) gRNA UMI counts.	required
response_n_nonzero	ndarray | None	(n_cells,) number of genes detected per cell, or None.	None
response_n_umis	ndarray | None	(n_cells,) gene expression UMIs per cell, or None. Give both response covariates or neither; with neither, the design uses the gRNA covariates only.	None
extra_covariates	DataFrame | None	Optional per-cell covariates, one row per cell. Categorical levels order as pandas categories, or sorted as strings.	None

Returns:

Type	Description
tuple[ndarray, list[str]]	(X, column_names), X of shape (n_cells, p).

Raises:

Type	Description
ValueError	only one response covariate given, lengths disagree, an extra covariate uses a name sceptre reserves, or the design is not finite and full rank.

pysceptre.assignment.design.design_from_covariates

design_from_covariates(covariates: DataFrame) -> tuple[ndarray, list[str]]

sceptre's default formula applied to a covariate frame, as model.matrix would.

covariates is sceptre's covariate_data_frame, columns in its order. response_p_mito is left out; a column whose name contains n_umis or n_nonzero enters as log(x), or log(x + 1) when any value is 0; any other column with 15 or more distinct values is left out; the rest enter as they are, categorical ones as treatment dummies (first level as reference). An intercept comes first.

Returns:

Type	Description
tuple[ndarray, list[str]]	(X, column_names).

Raises:

Type	Description
ValueError	the design is not finite and full rank.

pysceptre.assignment.design.cell_count_covariates

cell_count_covariates(matrix: spmatrix | ndarray) -> tuple[ndarray, ndarray]

Per-cell (n_nonzero, n_umis) of a (n_features, n_cells) count matrix.

As sceptre's compute_cell_covariates: n_nonzero counts the entries above 0.5 and n_umis is the column sum.

pysceptre.assignment.fishash.impute_masked_counts

impute_masked_counts(counts: spmatrix | ndarray, mask: spmatrix | ndarray, *, eps: float = 0.0001, max_iter: int = 10) -> ImputedCounts

Replace masked entries of a count matrix by a rank-one Poisson fit to the rest.

Port of fishash's impute_masked_counts. The unmasked counts are fitted as guide_freqs[g] * cell_sizes[c] by alternating closed-form updates, until both factors move less than eps on the log scale or max_iter iterations have run.

Parameters:

Name	Type	Description	Default
counts	spmatrix | ndarray	(n_grnas, n_cells) integer counts, scipy.sparse or dense.	required
mask	spmatrix | ndarray	Same shape; nonzero (or True) marks an entry to replace. May mark zero counts.	required
eps	float	Convergence tolerance on the largest change in log factor.	0.0001
max_iter	int	Maximum alternating iterations.	10

Returns:

Type	Description
ImputedCounts	An ImputedCounts.

Raises:

Type	Description
ValueError	shapes differ, or counts is not a count matrix.

pysceptre.assignment.hypergeom.log_phyper

log_phyper(q: ArrayLike, m: ArrayLike, n: ArrayLike, k: ArrayLike, *, lower_tail: bool = True) -> ndarray

R's phyper(q, m, n, k, lower.tail, log.p = TRUE), elementwise.

X is the number of white balls in k draws without replacement from an urn with m white and n black balls (R's argument names).

Parameters:

Name	Type	Description	Default
q	ArrayLike	Quantiles; floored after adding 1e-7, as R does.	required
m	ArrayLike	Number of white balls. Rounded to the nearest integer, half to even.	required
n	ArrayLike	Number of black balls. Rounded separately from m.	required
k	ArrayLike	Number of draws. Rounded.	required
lower_tail	bool	If True, log P(X <= q); otherwise log P(X > q).	True

Returns:

Type	Description
ndarray	Log probabilities, broadcast to the arguments' common shape. NaN where R returns NaN
ndarray	(a NaN argument, or rounded arguments outside the domain).

pysceptre.assignment.cells.cells_by_grna

cells_by_grna(assigned: spmatrix | ndarray, grna_ids: Sequence[str]) -> dict[str, ndarray]

Each gRNA's assigned cells.

Parameters:

Name	Type	Description	Default
assigned	spmatrix | ndarray	(n_grnas, n_cells) boolean matrix, scipy.sparse or dense, as returned in FishashResult.assigned or MixtureResult.assigned.	required
grna_ids	Sequence[str]	The row names.	required

Returns:

Type	Description
dict[str, ndarray]	{grna_id: ascending 0-based int64 cell indices}, in grna_ids order, an empty array
dict[str, ndarray]	for a gRNA assigned to no cell.

pysceptre.assignment.cells.cells_by_target

cells_by_target(grna_cells: dict[str, ndarray], grna_target_data_frame: DataFrame) -> tuple[dict[str, ndarray], dict[str, ndarray]]

Group per-gRNA cell sets into the two dicts run_discovery_analysis takes.

Parameters:

Name	Type	Description	Default
grna_cells	dict[str, ndarray]	{grna_id: cell indices}, as from cells_by_grna.	required
grna_target_data_frame	DataFrame	Columns grna_id and grna_target; a target of "non-targeting" marks a non-targeting gRNA.	required

Returns:

Type	Description
dict[str, ndarray]	(grna_target_cells, ntc_grna_cells). The first maps each targeting target to the union
dict[str, ndarray]	of its gRNAs' cells (ascending), keys in order of first appearance in the data frame; the
tuple[dict[str, ndarray], dict[str, ndarray]]	second maps each non-targeting gRNA to its own cells. A gRNA in the data frame but not in
tuple[dict[str, ndarray], dict[str, ndarray]]	grna_cells contributes no cells.

Raises:

Type	Description
KeyError	a required column is missing.

pysceptre.assignment.cells.cells_w_zero_or_twoplus_grnas

cells_w_zero_or_twoplus_grnas(assigned: spmatrix | ndarray) -> ndarray

The cells sceptre's low-MOI QC removes after a thresholding, mixture or fishash assignment.

assign_grnas_maximum returns its own set, which follows other rules.

Parameters:

Name	Type	Description	Default
assigned	spmatrix | ndarray	(n_grnas, n_cells) boolean matrix, scipy.sparse or dense.	required

Returns:

Type	Description
ndarray	Ascending 0-based int64 indices of the cells assigned no gRNA, or two or more.

The dose test

Not from sceptre. Each cell carries a target with a weight set by its gRNA UMI count instead of being called; the three analyses take the weights through grna_target_weights (ntc_grna_weights for the calibration check). See Design decisions.

pysceptre.assignment.dose.dose_weights

dose_weights(grna_matrix: spmatrix | ndarray, grna_ids: Sequence[str], grna_target_data_frame: DataFrame, *, floor: float | str = 'auto', ceiling: float = 500.0) -> DoseWeights

Each target's cells and dose weights, from raw gRNA UMI counts.

A cell's count for a target is its largest count over the target's gRNAs; the cell is kept when that count is above floor and weighted by dose_ramp. A non-targeting gRNA is kept on its own, as run_calibration_check regroups them.

Parameters:

Name	Type	Description	Default
grna_matrix	spmatrix | ndarray	(n_grnas, n_cells) raw integer UMI counts, scipy.sparse or dense.	required
grna_ids	Sequence[str]	The row names.	required
grna_target_data_frame	DataFrame	Columns grna_id and grna_target; a target of "non-targeting" marks a non-targeting gRNA. A gRNA listed here but absent from grna_ids contributes no cells.	required
floor	float | str	Counts at or below this get no weight. "auto" (default) estimates it from grna_matrix with estimate_dose_floor's defaults, and at least 1.	'auto'
ceiling	float	Counts at or above this get weight 1.	500.0

Returns:

Type	Description
DoseWeights	A DoseWeights, whose floor is the one used. Targets keep their order of first
DoseWeights	appearance in the data frame.

Raises:

Type	Description
KeyError	a required column is missing.
ValueError	non-integer or negative counts, or not 0 < floor < ceiling.

pysceptre.assignment.dose.DoseWeights dataclass

Cells and weights for run_discovery_analysis, run_power_check and run_calibration_check.

Attributes:

Name	Type	Description
grna_target_cells	dict[str, ndarray]	{target: ascending cell indices}, the cells whose largest count over the target's gRNAs is above the floor.
grna_target_weights	dict[str, ndarray]	{target: weights}, aligned with grna_target_cells.
ntc_grna_cells	dict[str, ndarray]	{non-targeting gRNA: ascending cell indices}, each gRNA on its own.
ntc_grna_weights	dict[str, ndarray]	{non-targeting gRNA: weights}, aligned with ntc_grna_cells.
floor	float	Counts at or below this get no weight.
ceiling	float	Counts at or above this get weight 1.

pysceptre.assignment.dose.dose_ramp

dose_ramp(counts: ndarray, *, floor: float, ceiling: float) -> ndarray

The dose weight of each count: clip(log(c / floor) / log(ceiling / floor), 0, 1).

Zero at or below floor, one at or above ceiling, linear in log(c) between.

Raises:

Type	Description
ValueError	unless 0 < floor < ceiling.

pysceptre.assignment.dose.estimate_dose_floor

estimate_dose_floor(grna_matrix: spmatrix | ndarray, *, real_min: int = 50, max_count: int = 30, threshold: float = 0.5) -> DoseFloor

The dose test's floor from the gRNA counts alone: where single-UMI-like noise ends.

For each count k up to max_count, the entries with exactly k UMIs are spread over the gRNAs as a mixture of the single-UMI entries' spread and the spread of entries with at least real_min UMIs; the mixing share is fitted by maximum likelihood. The floor is the largest k whose share is at least threshold.

Parameters:

Name	Type	Description	Default
grna_matrix	spmatrix | ndarray	(n_grnas, n_cells) raw integer UMI counts, scipy.sparse or dense.	required
real_min	int	Entries with at least this many UMIs stand for real guides.	50
max_count	int	The largest count examined.	30
threshold	float	The noise share a count's entries must reach to fall at or below the floor.	0.5

Raises:

Type	Description
ValueError	non-integer or negative counts, no single-UMI entries, or none with at least real_min UMIs.

pysceptre.assignment.dose.DoseFloor dataclass

The noise share of gRNA entries at each count, and the floor it gives.

Attributes:

Name	Type	Description
floor	int	The largest count whose entries are at least threshold noise; 0 if none is.
noise_share	Series	Indexed by count k, the share of the entries with exactly k UMIs whose spread over gRNAs matches the single-UMI entries rather than the real ones.
n_entries	Series	Indexed by count k, how many (gRNA, cell) entries have exactly k UMIs.

Lower-level building blocks

Most users will not need to go below run_discovery_analysis. Each pipeline stage is independently importable and unit-tested, for anyone extending or debugging the pipeline.

The GLM layer

pysceptre.glm.irls.fit_poisson_glm_batch

fit_poisson_glm_batch(X: ndarray, Y: ndarray, *, eps: float = _EPS, maxit: int = _MAXIT, X_outer_flat: ndarray | None = None, mu_floor: float = _MU_FLOOR) -> GlmFitBatchResult

X: (n, p) shared design matrix. Y: (k, n) or (n,) response(s), one per row.

X_outer_flat: optional, from x_outer_flat(X). Pass it when fitting repeatedly against the same design matrix; it is rebuilt per call otherwise.

mu_floor: lower bound on a fitted mean. R's glm.fit bounds it at the machine epsilon; the default here is 1e-10. See docs/design.md, "One Poisson GLM per gRNA".

pysceptre.glm.irls.fit_binomial_glm_batch

fit_binomial_glm_batch(X: ndarray, Y: ndarray, *, eps: float = _EPS, maxit: int = _MAXIT, X_outer_flat: ndarray | None = None) -> GlmFitBatchResult

X: (n, p) shared design matrix. Y: (k, n) or (n,) 0/1 indicator(s), one per row.

X_outer_flat: optional, from x_outer_flat(X). Pass it when fitting repeatedly against the same design matrix; it is rebuilt per call otherwise.

pysceptre.glm.nb_theta.estimate_theta

estimate_theta(y: ndarray, mu: ndarray, dfr: float, limit: int = 50, eps: float = eps ** 0.25) -> tuple[float, int]

Returns (theta_estimate, method) with method in {1: MLE, 2: MM, 3: pilot}.

Precomputation and resampling

pysceptre.precompute.pieces.compute_precomputation_pieces

compute_precomputation_pieces(expression_vector: ndarray, covariate_matrix: ndarray, fitted_coefs: ndarray, theta: float, full_test_stat: bool = True) -> PrecomputationPieces

pysceptre.crt.sampler.crt_index_sampler

crt_index_sampler(fitted_probabilities: ndarray, B: int, rng: Generator, n_trt: int) -> list[ndarray]

The CRT draw for a target with n_trt treated cells among fitted_probabilities.

crt_index_sampler_exact when the target is more than EXACT_SAMPLER_ABOVE_SHARE of the cells, crt_index_sampler_fast otherwise. The choice depends only on the two counts, never on fitted values, so it cannot change with how a run is chunked.

pysceptre.crt.sampler.crt_index_sampler_fast

crt_index_sampler_fast(fitted_probabilities: ndarray, B: int, rng: Generator) -> list[ndarray]

Returns a length-B list of 0-based integer arrays: synthetic_idxs[b] holds the cells treated in synthetic draw b, and may list a cell more than once. See the module docstring for the sparse-generation strategy (O(B * n_trt) rather than O(B * n_cells)) and for when that repeat matters.

pysceptre.crt.sampler.crt_index_sampler_exact

crt_index_sampler_exact(fitted_probabilities: ndarray, B: int, rng: Generator) -> list[ndarray]

crt_index_sampler_fast without its approximation, for large probabilities.

Same return shape: a length-B list of 0-based cell-index arrays, each ascending. Cell j is in each draw independently with probability fitted_probabilities[j], which is the law of R's binomial count plus without-replacement placement, and no cell is ever listed twice in a draw. Work stays proportional to the number of inclusions: each cell's draw positions are the partial sums of geometric gaps.

For the NT-cells control group, where a target is a large share of its combined cells. See docs/design.md, "The CRT against the NT cells".

pysceptre.crt.permutations.permutation_draws

permutation_draws(n_cells: int, m: int, b: int, rng: Generator) -> ndarray

b random subsets of m distinct cells, as a (b, m) index array.

One Generator.choice(..., replace=False) per draw. That looks like the naive option and is in fact the fast one: numpy samples without replacement in C, so the Python loop costs one call per draw rather than one per cell. Measured 0.28 s for 30,497 draws of 692 cells from 567,690.

Two cleverer approaches were tried and rejected. Rejection sampling -- draw with replacement, redraw rows containing a repeat -- is elegant when m is tiny relative to n_cells, and collapses when it is not: at m = 120 from 1,500 cells the expected collision count per row is 4.8 and fewer than 1% of rows survive, so it cannot serve as a general sampler. Ranking rng.random((b, n_cells)) is O(b * n_cells), which is 17 billion values at real scale for 21 million wanted.

Each row is a random ordering, deliberately not sorted. That is what lets every target share these draws: a prefix of a uniformly random ordering of distinct cells is a uniformly random subset of that length, so a target with n_trt cells takes row[:n_trt]. Sorting here would break that silently -- prefixes would be the lowest-numbered cells rather than random ones -- so a test asserts the rows are not sorted.

Reimplemented rather than ported: sceptre's fisher_yates_samlper is C++ reached through .Call, and the draws could not match bit-for-bit anyway, since the two languages seed different generators.

pysceptre.crt.permutations.nested_permutation_draws

nested_permutation_draws(n_control: int, m: int, M: int, b: int, rng: Generator) -> ndarray

b draws for the NT-cells control group, as a (b, M) index array.

For every k in [m, M], the first k entries of each row are a uniformly random k-subset of {0, ..., n_control + k - 1}. A target with n_trt treated cells is tested against n_control NT cells, so its combined vector has n_control + n_trt positions and its draws are the first n_trt entries of each row: one shared array serves every target whose size lies in [m, M].

Same law as sceptre's hybrid_fisher_iwor_sampler(N, m, M, B), built backwards from a uniform ordering instead of forwards from a shuffle. See docs/design.md, "Permutations against the NT cells".

The test statistic

pysceptre.test_statistic.score_stat.compute_observed_full_statistic

compute_observed_full_statistic(a: ndarray, w: ndarray, D: ndarray, trt_idxs: ndarray, trt_weights: ndarray | None = None) -> float

The score statistic of the treated cells, or of the dose test's weighted cells.

With trt_weights (one per treated cell; the dose test), the sums become sum t a, sum t^2 w and sum_k (sum t D[k])^2; weights all equal to 1 give the unweighted statistic exactly. docs/design.md, "The dose test".

pysceptre.test_statistic.score_stat.compute_null_full_statistics

compute_null_full_statistics(a: ndarray, w: ndarray, D: ndarray, synthetic_idxs: list[ndarray]) -> ndarray

synthetic_idxs: a list of B (possibly ragged, i.e. differently-sized) 0-based integer index arrays -- one synthetic treated-cell set per resample. This targets sceptre's use_all_cells=True case (the only one the complement-control-group CRT path uses): each resample's own set size is used directly, rather than a fixed external n_trt.

Convenience wrapper around compute_null_full_statistics_flat for callers with a single synthetic_idxs list to consume once; prefer flattening once and calling compute_null_full_statistics_flat directly when the same synthetic_idxs will be reused across many (a, w, D) triples.

pysceptre.test_statistic.empirical_p.compute_empirical_p_value

compute_empirical_p_value(null_statistics: ndarray, z_orig: float, side: int) -> float

pysceptre.test_statistic.skew_normal.fit_and_evaluate_skew_normal

fit_and_evaluate_skew_normal(z_orig: float, null_statistics: ndarray, side_code: int) -> SkewNormalEvalResult

pysceptre.test_statistic.fold_change.estimate_log_fold_change

estimate_log_fold_change(y: ndarray, mu: ndarray, trt_idxs: ndarray, trt_weights: ndarray | None = None) -> tuple[float, float]

trt_idxs: 0-based indices of treated cells. Returns (fold_change, se_fold_change).

With trt_weights (the dose test), the fold change at full weight, with the effect taken as linear in the weight: 1 + sum t (y - mu) / sum t^2 mu. Weights all equal to 1 give the unweighted estimate exactly. docs/design.md, "The dose test".

pysceptre.test_statistic.resampling.run_low_level_test_full

run_low_level_test_full(y: ndarray, mu: ndarray, a: ndarray, w: ndarray, D: ndarray, trt_idxs: ndarray, synthetic_idxs, *, stacked: ndarray | None = None, B1: int = 499, B2: int = 4999, B3: int = 0, fit_parametric_curve: bool = True, side_code: int = 0, p_thresh: float = P_THRESH, return_resampling_dist: bool = False, null_statistics_fn: Callable[[int, int], ndarray | None] | None = None, trt_weights: ndarray | None = None) -> PairResult

0-based observed treated-cell indices. With none, or with a

statistic that is not finite, the p-value is NaN rather than a test.

the B1+B2+B3 draws, consumed in three consecutive slices

[0:B1], [B1:B1+B2], [B1+B2:B1+B2+B3]. Normally a StagedDraws, which materializes each slice only when that stage is reached and memoizes it, so a target's draws are built once per stage rather than once per pair -- and the stages a run never reaches are never built at all. A CSR matrix or a list of ragged index arrays is also accepted and wrapped, since those are the natural things to write by hand in a test.

stacked: stack_pieces(a, w, D), built once per gene and reused across that gene's targets. Computed here when omitted. null_statistics_fn: an optional faster route to a stage's null statistics, taking (lo, hi). Permutations pass PermutationPrefixSums.statistics, which reads a column out of one per-gene prefix-sum array instead of running a per-target matmul -- possible only because every target shares the permutation rows. Returning None for a stage means "no cheaper route here", and the draw matrix is used for it, so a provider can decline the stages whose arrays would be too large.

Orchestration

pysceptre.pipeline.discovery.run_discovery_ntcells_complement

run_discovery_ntcells_complement(response_matrix, gene_ids: list[str], covariate_matrix: ndarray, grna_target_cells: dict[str, ndarray], pairs: DataFrame, *, B1: int = 499, B2: int = 4999, B3: int = 0, fit_parametric_curve: bool = True, side_code: int = 0, seed: int | None = None, target_chunk_size: int = _DEFAULT_TARGET_CHUNK_SIZE, chunk_memory_gb: float = _DEFAULT_CHUNK_MEMORY_GB, n_jobs: int = 1, resampling_mechanism: str = 'crt', permutations: ndarray | None = None, permutation_width: int | None = None, grna_target_weights: dict[str, ndarray] | None = None, glm_solver: str = 'outer') -> DataFrame

pairs: DataFrame with columns 'response_id', 'grna_target' -- the QC-passed pairs to test. Returns a 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.

Under resampling_mechanism="permutations" the shared draws are a (B1 + B2 + B3, width) array of cell indices whose rows are random orderings, and each target reads the first n_trt entries of each row. permutation_width sets the width, which defaults to the largest target in grna_target_cells; permutations supplies the draws themselves instead of drawing them, for example R's own.

Targets are fit and CRT-drawn in chunks of target_chunk_size rather than all at once: each target's B1+B2+B3 synthetic index sets are individually small, but holding every target's draws in memory simultaneously does not scale -- at real dataset sizes (~3,000 targets x ~5,500 resamples) this was measured to OOM-kill the process. Chunking keeps peak memory bounded to O(chunk_size) while still batching the (bulk of the) per-target logistic fit and CRT draw across many targets at once for speed, not falling back to a slow one-target-at-a-time loop.

target_chunk_size is an upper bound, not a mandate: it is reduced automatically so the working arrays stay under chunk_memory_gb, so no chunk size a caller passes can blow up memory. Both per-chunk costs are linear in the chunk size -- the dense binomial-fit arrays and the CRT draws -- which is why one budget covers both. That matters most for resampling_approximation="no_approximation", whose B3 grows as n_pairs / multiple_testing_alpha and reaches 1.65M draws (5.2 GB) per target at real dataset scale.

pysceptre.pipeline.discovery.run_discovery_nt_cells

run_discovery_nt_cells(response_matrix, gene_ids: list[str], covariate_matrix: ndarray, grna_target_cells: dict[str, ndarray], nt_cells: ndarray, pairs: DataFrame, *, B1: int = 499, B2: int = 4999, B3: int = 0, fit_parametric_curve: bool = True, side_code: int = 0, seed: int | None = None, target_chunk_size: int = _DEFAULT_TARGET_CHUNK_SIZE, chunk_memory_gb: float = _DEFAULT_CHUNK_MEMORY_GB, n_jobs: int = 1, resampling_mechanism: str = 'crt', permutations: ndarray | None = None) -> DataFrame

The discovery test with the NT cells as the control group, sceptre's low-MOI default.

Each pair is tested on its target's cells together with the NT cells, so the gene's response GLM is fit once per pair rather than once per gene. nt_cells is the NT pool as 0-based cell indices, in the order R keeps all_nt_idxs; it must not share a cell with any tested target.

Under the CRT each target draws its resamples from a logistic fit over its own combined cells. Under permutations one shared (B1 + B2 + B3, M) array serves every target, from nested_permutation_draws sized by the smallest and largest nonempty target in grna_target_cells; permutations supplies that array instead (R's own hybrid_fisher_iwor_sampler output, say).

Takes and returns the same things as run_discovery_ntcells_complement otherwise. See docs/design.md, "Low MOI and the NT-cells control group".

pysceptre.pipeline.api.resolve_analysis_settings

resolve_analysis_settings(moi: str, control_group: str | None, resampling_mechanism: str | None) -> tuple[str, str]

(control_group, resampling_mechanism) with R's defaults filled in.

Ports the defaults of sceptre's set_analysis_parameters: high MOI uses the complement control group and the CRT; low MOI uses the NT cells and permutations. Either can be overridden, except that the NT cells are refused in high MOI. R silently replaces any high-MOI control group with the complement; refusing is stricter, so an explicit setting is never ignored.

Raises:

Type	Description
ValueError	an unknown value, or control_group="nt_cells" with moi="high".

pysceptre.pipeline.api.nt_cell_pool

nt_cell_pool(ntc_grna_cells: dict[str, ndarray], n_cells: int) -> ndarray

The NT cells, as R's all_nt_idxs: each NT gRNA's cells in iteration order.

Each gRNA's cells are sorted and concatenated in the dict's order, which is R's order when the dict is built in R's gRNA order. The order decides which cell each resample position refers to, so a different order is a different Monte Carlo draw: every p-value moves at the resampling resolution, and fitted values in their last bits. Build the dict in a fixed order.

Raises:

Type	Description
ValueError	no NT cells, an index out of range, or a cell under two NT gRNAs -- sceptre's low-MOI QC removes every cell with more than one gRNA, and the NT-cells analysis assumes it.

pysceptre.pipeline.discovery.parallel_backend

parallel_backend() -> str

The platform default backend for a worker pool: "fork" or "thread".

fork on Linux, thread everywhere else. PYSCEPTRE_BACKEND overrides it to either value; a request for fork off Linux is refused with a warning, because forking after Apple's Accelerate BLAS can deadlock.

Callers that build one pool per chunk should ask gene_job_backend instead, which may override this. See docs/design.md, "Choosing a parallel backend".

pysceptre.pipeline.discovery.gene_job_backend

gene_job_backend(n_chunks: int) -> str | None

Backend for the per-chunk gene pool, or None for the platform default.

Returns "thread" when the platform default is fork and the run has more than _THREAD_ABOVE_N_CHUNKS chunks, since one fork per chunk then costs more than the GIL does. None everywhere else, including on platforms whose default is already threads. See docs/design.md, "Choosing a parallel backend".

Analytical power primitives

The pieces compute_power is assembled from, each a direct port of the correspondingly named function in PerturbPlan (MIT -- see THIRD_PARTY_LICENSES).

pysceptre.analytical_power.closed_form.test_stat_distribution

test_stat_distribution(*, num_trt_cells: ndarray, num_cntrl_cells: ndarray, num_trt_cells_sq: ndarray, expression_mean: ndarray, expression_size: ndarray, fold_change_mean: ndarray, fold_change_sd: ndarray) -> tuple[ndarray, ndarray]

Mean and sd of the score statistic under the alternative.

fold_change_mean is a multiplier: a 15% knockdown is 0.85, not 0.15.

num_trt_cells is the sum of per-gRNA cell counts for the target, and num_trt_cells_sq the sum of their squares -- not the size of the union of perturbed cells. See docs/design.md, "Analytical per-pair power".

Returns (mean, sd), both broadcast to the shape of the inputs.

pysceptre.analytical_power.closed_form.qc_failure_prob

qc_failure_prob(*, fold_change_mean: ndarray, expression_mean: ndarray, expression_size: ndarray, num_trt_cells: ndarray, num_cntrl_cells: ndarray, n_nonzero_trt_thresh: int, n_nonzero_cntrl_thresh: int) -> ndarray

P(the pair fails sceptre's pairwise nonzero-count QC).

Both thresholds at 0 make this 0 everywhere, which is the setting for pairs that already passed QC in a real analysis.

pysceptre.analytical_power.closed_form.rejection_prob

rejection_prob(mean: ndarray, sd: ndarray, *, side: str, cutoff: float) -> ndarray

P(the test rejects), for a normal statistic with the given mean and sd.

cutoff is the nominal p-value threshold the test is read at, and it is not interchangeable with a screen's alpha: the validated configuration is side="left" at cutoff = alpha / 2, sceptre's p-value being two-sided while a knockdown is a one-sided claim.

pysceptre.analytical_power.closed_form.zero_prob

zero_prob(fold_change_mean: ndarray, expression_mean: ndarray, expression_size: ndarray) -> ndarray

P(count == 0) for an NB gene under a multiplicative fold change.

pysceptre.analytical_power.closed_form.var_nb

var_nb(mean: ndarray, size: ndarray) -> ndarray

Negative-binomial variance for a given mean and NB size.

size is theta, not the dispersion: pass 1 / dispersion if that is what you hold.