spacr.power_simulate

Simulator half of the spaCR power analysis — a Python port of spaCRPower.

Ported from spaCRPower (R/simulate_screen.R), Copyright (c) 2025 Matthew O’Meara (maom@umich.edu, ORCID 0000-0002-3128-5331), released under the MIT licence:

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

What this is for

A pooled Toxoplasma CRISPR screen has three noisy layers stacked on top of each other — which genotypes landed in which well, how many of a well’s cells the microscope saw and how well the classifier called them, and how many reads each genotype got out of the amplicon pool. The question a screener has to answer before spending the microscope time is “with this library, this many wells, this many cells per well and a classifier at 0.80/0.12, would I actually find a hit?”. That is a power analysis, and it is only answerable by simulation because no layer is analytically tractable once the other two are stacked on it.

This module is the simulator half: it produces the tidy (well, gene) table that the model half (a Poisson regression of per-well positive counts on per-gene log10 read fraction) consumes by column name. It deliberately does not fit anything.

Four stages plus an orchestrator, mirroring R/simulate_screen.R:

• simulate_library() — the genes, which of them are hits, and how unevenly the library is represented.

• simulate_spot_plate() — which genes landed in which wells.

• simulate_imaging_plate() — how many cells of each genotype were imaged in each well, and how many of them the classifier called positive.

• simulate_sequencing_plate() — the read count for each genotype in each well, after PCR amplification and a finite sequencing depth.

• simulate_screen() — all four, joined on (well, gene).

Reproducibility is not optional

Every public sampler takes exactly one of seed= or rng= and refuses to run with neither. The global numpy random state is never touched. A power analysis is a number somebody defends in a grant review or a methods section; one that cannot be re-run bit-for-bit is not evidence, so drawing silently from OS entropy is treated as a configuration error rather than a convenience.

Deviations from the R package, and why

The R sources are the specification, but four of their behaviours are defects that a line-by-line translation would faithfully reproduce. They are listed here, and each is repeated at the function that departs from upstream. See also proposals/SIM_PORT_PLAN.md §3 in this repository, which catalogues them.

1. imaging_n_cells_per_well did not exist. R/fit_model.R reads well_data$imaging_n_cells_per_well[1] as the Poisson offset, but simulate_imaging_plate emits only ..._mu, ..._var and ..._gene_per_well; R’s $ partial matching is ambiguous across those three and yields NULL, so the offset column vanishes and the simulate → fit path does not run as committed. simulate_imaging_plate() emits the column for real, as the realised per-well cell total.

2. Reads per well were divided by the gene count, not the well count. round(n_reads_total / nrow(well_data)) groups by well, so nrow is the library size; with 452 genes and n_reads_total = 128318 that is 284 reads per well against a real screen near 3e4. The parameter was also documented as “total reads” in one vignette and “geometric mean of well reads” in another. simulate_sequencing_plate() takes an unambiguous n_reads_per_well and derives n_reads_total.

3. The hypergeometric draw size could exceed the urn. Upstream guards rmvhyper with k = min(n_cells_in_well * pcr_factor, n_reads_per_well), but the urn is round(cells * pcr) computed element-wise, whose sum is not round(sum(cells) * pcr). The guard here is computed from the actual colour vector.

4. The COM-Poisson branch is unusable. COMPoissonReg::rcmp is called but not declared in the R package’s DESCRIPTION, and has no maintained Python equivalent inside spaCR’s dependency set. sample_count_mean_variance() replaces the nu dispersion knob with a mean/variance pair dispatched to negative-binomial / Poisson / binomial, which spans the same over-, equi- and under-dispersed range in closed form with no new dependency. Third and higher moments differ from COM-Poisson; nothing upstream used them.

One further departure is a choice, not a bug fix, and is exposed as a parameter: upstream splits a well’s imaged cells uniformly across the genes present (prob = gene_in_well is a 0/1 vector), ignoring how abundant each gene is. simulate_imaging_plate() defaults to imaging_split='abundance' and keeps 'uniform' for exact R parity.

Two things neither package modelled

The list above is about fidelity to the R. These two are about fidelity to a real screen, and both make the answer worse — which is why they are worth having and why they are off by default. A simulator whose baseline shifted under a version bump would make every power figure already quoted from it wrong, so turning either on is an explicit act.

• Sequencing error (misassign_reads(), sequencing_error_rate=). Both packages treat the read fraction as an exact record of which genotypes were in a well. Substitution errors, index hopping, PCR chimeras and mismatch-tolerant demultiplexing all credit reads to the wrong gene. Simulating it says the direct dilution is small — but that it silently disables the unidentified-gene check, turning genes that were correctly reported as untested into confident non-hits. On the reference design that costs ten times more than the dilution does. See misassign_reads().

• Well dropout from too few imaged cells (drop_low_cell_wells(), min_cells_per_well=). A well where the microscope found three cells enters the fit as one observation next to a well with four hundred. Its positive fraction can only be 0, 1/3, 2/3 or 1, and its standard error is several times the classifier’s whole signal gap. The Poisson offset stops it dominating the scale of the fit; it does not stop its read-fraction covariate being paired with noise. Dropping such wells is what an analyst does by hand, and it costs wells — so which way the trade comes out is a thing to simulate rather than assert.

Column names follow the R package exactly so the model half can consume them by name — including the upstream misspelling n_barcodes_per_genes_per_well, which is kept rather than quietly corrected, because renaming it would break the join for anyone reading against the R vignettes.

seealso:

spacr.sim for the older, unrelated in-house screen simulator.

Exceptions

AbundanceClippedWarning	A gene_abundance x well_abundance product exceeded 1 and was clipped.
ImpossibleMomentsError	The requested (mean, variance) pair is not attainable by the distribution.
MalformedPlateError	A plate frame handed to a downstream stage is not a valid plate.
PowerSimulationError	Base class for every error this simulator raises deliberately.
ScreenDesignError	The screen design is degenerate or out of range.
SequencingScaleError	The amplified barcode pool is too large to sample exactly.

Functions

drop_low_cell_wells(→ pandas.DataFrame)	Remove wells whose imaged cell count is too low to be informative.
misassign_reads(→ numpy.ndarray)	Credit a fraction of each well's reads to the wrong gene.
rbeta_mean_variance(→ numpy.ndarray)	Draw n beta variates parameterised by mean and variance.
rdirichlet_stable(→ numpy.ndarray)	Draw one Dirichlet vector, in log space so small alpha does not underflow.
resolve_rng(→ numpy.random.Generator)	Return the generator to draw from, insisting the caller chose one.
rgamma_mean_variance(→ numpy.ndarray)	Draw n gamma variates parameterised by mean and variance.
rnbinom_mean_variance(→ numpy.ndarray)	Draw n negative-binomial variates parameterised by mean and variance.
sample_count_mean_variance(→ numpy.ndarray)	Draw n non-negative integer counts with the requested mean and variance.
simulate_imaging_plate(→ pandas.DataFrame)	Simulate the imaged cells per genotype per well and the classifier's calls.
simulate_library(→ pandas.DataFrame)	Simulate the perturbation library: which genes exist, how abundant, which are hits.
simulate_screen(→ pandas.DataFrame)	Run all four stages and return one tidy (well, gene) table.
simulate_sequencing_plate(→ pandas.DataFrame)	Simulate the barcode read counts for each genotype in each well.
simulate_spot_plate(→ pandas.DataFrame)	Simulate which genes landed in which wells.

Module Contents

exception spacr.power_simulate.AbundanceClippedWarning

Bases: UserWarning

A gene_abundance x well_abundance product exceeded 1 and was clipped.

That product is used as a Bernoulli probability. R’s rbinom returns NA for prob > 1 with a warning nobody reads, which poisons every downstream count. Clipping keeps the run alive, but it means the realised genes-per-well is below the requested well_abundance_factor_mu for the affected genes, so the warning carries the clip count and the worst offender.

Initialize self. See help(type(self)) for accurate signature.

exception spacr.power_simulate.ImpossibleMomentsError

Bases: PowerSimulationError

The requested (mean, variance) pair is not attainable by the distribution.

Raised instead of clamping. A beta distribution cannot have variance above mean * (1 - mean), and a negative binomial cannot have variance below its mean; quietly moving the request onto the nearest feasible point would hand back samples from a distribution the caller did not ask for and never find out about — which is precisely how a power analysis ends up defending a number it did not compute.

Initialize self. See help(type(self)) for accurate signature.

exception spacr.power_simulate.MalformedPlateError

Bases: PowerSimulationError

A plate frame handed to a downstream stage is not a valid plate.

Missing a column the stage needs, or not exactly one row per (well, gene) pair. The downstream stages pivot the frame to a gene-by-well matrix, and a duplicated or missing pair would either silently drop observations or fill them with NaN that propagates into counts.

Initialize self. See help(type(self)) for accurate signature.

exception spacr.power_simulate.PowerSimulationError

Bases: spacr.errors.SpacrError

Base class for every error this simulator raises deliberately.

Subclasses spacr.errors.SpacrError so callers can catch spaCR’s own failures separately from an incidental ValueError out of numpy.

Initialize self. See help(type(self)) for accurate signature.

exception spacr.power_simulate.ScreenDesignError

Bases: PowerSimulationError

The screen design is degenerate or out of range.

Zero genes, zero wells, a hit rate outside [0, 1], a non-positive abundance concentration. R’s 1:n_wells idiom turns n_wells = 0 into the two-element vector c(1, 0) and simulates two wells; there is no silently-empty frame to return here, so we refuse the design instead.

Initialize self. See help(type(self)) for accurate signature.

exception spacr.power_simulate.SequencingScaleError

Bases: PowerSimulationError

The amplified barcode pool is too large to sample exactly.

See MAX_HYPERGEOMETRIC_URN. Almost always means pcr_factor_mu was given on the linear scale when it is a log-scale parameter.

Initialize self. See help(type(self)) for accurate signature.

spacr.power_simulate.drop_low_cell_wells(screen: pandas.DataFrame, min_cells_per_well: int, *, drop: bool = True) → pandas.DataFrame

Remove wells whose imaged cell count is too low to be informative.

The second thing neither spaCRPower nor the port did. A well where the microscope found three cells produces a positive fraction that can only take the values 0, 1/3, 2/3 and 1; its binomial standard error is around 0.27, which is three times the entire gap between the classifier’s hit-cell and background rates. It carries essentially no information about which genotypes were in it — and it enters the fit as one more observation alongside a well with four hundred cells.

The Poisson offset does part of the job: a well with three cells has a small expected count, so it does not dominate the scale of the fit. What the offset does not do is stop the well’s read-fraction covariate from being paired with a response that is almost pure noise, and a screen with a long tail of thin wells is a screen whose covariate-response relationship is being averaged against nothing.

Dropping them is what an analyst does by hand, and it costs something: fewer wells is less power. Simulating both sides of that trade is the reason this is a parameter rather than a fixed rule.

Whole wells go, never single (gene, well) rows. A partially dropped well would leave the well’s positive total and its cell total describing different sets of genes, which is a table the model half is entitled to assume cannot exist.

Parameters:

• screen – joined screen table; needs well and one of imaging_n_cells_per_well / imaging_n_cells_per_gene_per_well.

• min_cells_per_well – wells with fewer imaged cells than this are removed. 0 disables the filter and returns the input unchanged.

• drop – False annotates with well_kept and removes nothing, for a caller that wants to see what would go.

Returns:

a new frame carrying a boolean well_kept column, filtered when drop. frame.attrs['n_wells_dropped'] and ['n_wells_before'] record the cost.

Raises:

• MalformedPlateError – If neither cell-count column is present, or a well’s imaging_n_cells_per_well is not constant within the well.

• ScreenDesignError – If min_cells_per_well is negative.

spacr.power_simulate.misassign_reads(reads: numpy.ndarray, sequencing_error_rate: float, *, rng: numpy.random.Generator | None = None, seed: int | None = None) → numpy.ndarray

Credit a fraction of each well’s reads to the wrong gene.

Neither spaCRPower nor the port had this stage, and its absence is optimistic in a specific way: it makes the read fraction an exact measure of which genotypes were in a well. In a real amplicon screen it is not. A barcode arrives at the wrong gene through base-call substitutions, index hopping on a patterned flow cell, PCR chimeras across the pooled reaction, and mismatch-tolerant demultiplexing — and every one of those mechanisms moves reads toward the middle, because the destination is drawn from the whole library rather than from the wells the source genotype was actually in.

Mis-assignment gives a gene phantom reads in wells it never entered, so the covariate the model regresses positive counts on is a shrunk-toward-uniform version of the truth. The obvious consequence is regression dilution — an attenuated coefficient and a smaller apparent effect — and measuring it in this simulator says that at realistic rates it is small: on the 452-gene reference design, 0.5 % mis-assignment moves the hit/non-hit separation from 0.799 to 0.794 among the genes that were identifiable to begin with. Deep sequencing averages most of it away.

The consequence that is not obvious, and is much larger, is what it does to the genes that were never testable. A gene that landed in every well, or in none, has a constant read-fraction column; its coefficient is confounded with the intercept and power_model.prepare_model_data reports it as unidentified rather than as a non-hit. That check is the one honest thing standing between a thin design and a page of confident negative results — and mis-assignment defeats it, because phantom reads give every gene a covariate that varies. On the same reference design, 0.5 % error takes the scored library from 317 genes to all 452: the 135 genes that were correctly flagged “untested” become scored non-hits with a covariate made entirely of noise, and the screen-wide separation falls from 0.799 to 0.723 — fourteen times the direct dilution, and all of it from a safeguard being switched off rather than from any loss of signal.

So the reason to simulate this is not to price the noise. It is that a screen with sequencing error and a screen without it disagree about how many genes were tested, and only one of them is telling the truth.

The model, per well:

1. Each gene’s reads survive independently with probability 1 - sequencing_error_rate.

2. Everything lost is pooled and redistributed uniformly across the whole library, including genes absent from the well and including the gene it came from.

Uniform, not abundance-weighted, because the destination of a misread barcode is set by which barcodes are one edit away from it, not by how much of the source was in the tube. And because self-assignment is allowed, the effective mis-assignment rate is sequencing_error_rate * (1 - 1/n_genes); at a 452-gene library that is a 0.2 % correction, and keeping it makes the read total exactly conserved, which is a much more useful invariant to be able to assert.

Cross-well hopping — a read landing in the wrong sample entirely — is not modelled. It needs a plate-level index layout the simulator does not carry, and within-well gene confusion is the mechanism that dilutes the effect being measured.

Parameters:

• reads – (n_genes, n_wells) integer read counts.

• sequencing_error_rate – probability a read is credited elsewhere, in [0, 1]. 0.0 returns the input unchanged.

• rng – Generator to draw from; mutually exclusive with seed.

• seed – Seed for a fresh generator; mutually exclusive with rng.

Returns:

a new (n_genes, n_wells) array. Every column sums to exactly what it summed to before.

Raises:

ScreenDesignError – If the rate is outside [0, 1], the array is not 2-D or holds a negative count, or neither/both of rng and seed are given.

Example:

Total reads are conserved, and at a full error rate nothing of the original assignment survives except by chance:

>>> counts = np.array([[100, 0], [0, 100]])
>>> out = misassign_reads(counts, 1.0, seed=3)
>>> [int(column.sum()) for column in out.T]
[100, 100]

spacr.power_simulate.rbeta_mean_variance(n: int, mean: float, var: float, *, rng: numpy.random.Generator | None = None, seed: int | None = None) → numpy.ndarray

Draw n beta variates parameterised by mean and variance.

Inverting mean = a / (a + b) and var = ab / ((a+b)**2 (a+b+1)) gives a = mean * (mean(1-mean)/var - 1) and b = (1-mean) * (same).

var == 0 is accepted and returns a constant array at mean. That is the exact limit of the family as a, b -> inf with the mean held, and it is the only way to express a perfect classifier — class_pos_mu=1.0, class_pos_var=0.0 — which is the single most useful configuration for testing the imaging stage, since every hit cell must then be called positive.

Parameters:

• n – Number of variates.

• mean – Target mean, must lie in [0, 1].

• var – Target variance, must lie in [0, mean * (1 - mean)).

• rng – Generator to draw from; mutually exclusive with seed.

• seed – Seed for a fresh generator; mutually exclusive with rng.

Returns:

Float array of shape (n,) with values in [0, 1].

Raises:

• ImpossibleMomentsError – If var is negative, or if var >= mean * (1 - mean). That bound is the variance of the Bernoulli with the same mean and is the supremum over all distributions on [0, 1]; no beta reaches it. Note this makes mean of exactly 0 or 1 admissible only with var == 0. Clamping instead of raising would return a distribution with the wrong spread and, for a classifier accuracy, silently change the effect size the whole power analysis is measuring.

• ScreenDesignError – If mean is outside [0, 1] or not finite, n is negative, or neither/both of rng and seed given.

Example:

>>> p = rbeta_mean_variance(200000, mean=0.8, var=0.01, seed=2)
>>> bool(abs(p.mean() - 0.8) < 0.005 and abs(p.var() - 0.01) < 0.001)
True
>>> rbeta_mean_variance(3, mean=1.0, var=0.0, seed=2).tolist()
[1.0, 1.0, 1.0]

spacr.power_simulate.rdirichlet_stable(alpha: float | Sequence[float] | numpy.ndarray, n_categories: int | None = None, *, rng: numpy.random.Generator | None = None, seed: int | None = None) → numpy.ndarray

Draw one Dirichlet vector, in log space so small alpha does not underflow.

numpy.random.Generator.dirichlet draws independent gammas and normalises them. At small concentration the individual gamma draws underflow to exactly 0.0, so a gene gets an abundance of exactly zero — measured here, at alpha = 0.05 with 452 categories, numpy produces exact zeros in roughly one draw in ten, and this routine produces none. Zero abundance is not merely unlucky: it removes the gene from every well, and the downstream read fraction for that gene becomes 0 / 0 in any well where nothing else amplified.

The fix is Marsaglia and Tsang’s boosting identity — for a < 1, G(a) =d= G(a + 1) * U**(1/a) — evaluated in logs, followed by a log-sum-exp normalisation. Underflow then requires the log weight to overflow, which is ~700 times further away.

Parameters:

• alpha – Either a scalar concentration (then n_categories is required) or a 1-D array of per-category concentrations.

• n_categories – Number of categories when alpha is a scalar.

• rng – Generator to draw from; mutually exclusive with seed.

• seed – Seed for a fresh generator; mutually exclusive with rng.

Returns:

Float array of shape (n_categories,) summing to 1.

Raises:

ScreenDesignError – If any concentration is not finite and positive, if n_categories is missing for a scalar alpha or is not positive, or if neither/both of rng and seed given.

Example:

>>> w = rdirichlet_stable(0.6, 452, seed=3)
>>> bool(abs(w.sum() - 1.0) < 1e-12 and (w > 0).all())
True

spacr.power_simulate.resolve_rng(rng: numpy.random.Generator | None = None, seed: int | None = None) → numpy.random.Generator

Return the generator to draw from, insisting the caller chose one.

Exactly one of rng or seed must be supplied.

Parameters:

• rng – An existing numpy.random.Generator to advance. Use this to thread one stream through several stages so the whole screen is one reproducible draw.

• seed – Integer seed for a fresh numpy.random.default_rng().

Returns:

The generator to sample from.

Raises:

ScreenDesignError – If neither or both were supplied. Neither would mean seeding from OS entropy, which makes the result impossible to reproduce; both would make it ambiguous which one actually applied, and a caller who passed a seed would reasonably but wrongly believe the run was pinned by it.

Example:

>>> a = resolve_rng(seed=0).normal(size=3)
>>> b = resolve_rng(seed=0).normal(size=3)
>>> bool((a == b).all())
True

spacr.power_simulate.rgamma_mean_variance(n: int, mean: float, var: float, *, rng: numpy.random.Generator | None = None, seed: int | None = None) → numpy.ndarray

Draw n gamma variates parameterised by mean and variance.

A gamma with shape k and rate r has mean = k / r and var = k / r**2; inverting gives r = mean / var and k = mean**2 / var.

The R-to-numpy trap this function exists to close: R’s rgamma takes a rate, numpy’s takes a scale, and scale = 1 / rate = var / mean. A port that passes the rate to numpy gets a distribution wrong by a factor of var**2 / mean**2 in the variance while still looking entirely plausible.

Parameters:

• n – Number of variates.

• mean – Target mean, must be finite and > 0.

• var – Target variance, must be finite and > 0.

• rng – Generator to draw from; mutually exclusive with seed.

• seed – Seed for a fresh generator; mutually exclusive with rng.

Returns:

Float array of shape (n,).

Raises:

ScreenDesignError – If mean or var is not finite and positive, or if n is negative, or if neither/both of rng and seed given.

Example:

>>> x = rgamma_mean_variance(200000, mean=4.0, var=2.0, seed=0)
>>> bool(abs(x.mean() - 4.0) < 0.05 and abs(x.var() - 2.0) < 0.05)
True

spacr.power_simulate.rnbinom_mean_variance(n: int, mean: float, var: float, *, rng: numpy.random.Generator | None = None, seed: int | None = None) → numpy.ndarray

Draw n negative-binomial variates parameterised by mean and variance.

With size = mean**2 / (var - mean) and prob = mean / var the distribution has exactly the requested first two moments. numpy’s negative_binomial(n, p) and R’s rnbinom(size, prob) agree on the convention (both have mean n(1-p)/p), so only the moment inversion needed porting.

Parameters:

• n – Number of variates.

• mean – Target mean, must be finite and > 0.

• var – Target variance. Must be strictly greater than mean.

• rng – Generator to draw from; mutually exclusive with seed.

• seed – Seed for a fresh generator; mutually exclusive with rng.

Returns:

Integer array of shape (n,).

Raises:

• ImpossibleMomentsError – If var <= mean. The negative binomial is over-dispersed by construction: at var == mean the size parameter is +inf (the Poisson limit, which numpy cannot sample as an NB), and below it the size is negative. R’s assertthat allows var == mean and then hands Inf to rnbinom, which returns NA — a whole column of silently missing counts. Use sample_count_mean_variance() if you want the equi- and under-dispersed cases handled for you.

• ScreenDesignError – If mean is not finite and positive, n is negative, or neither/both of rng and seed given.

Example:

>>> x = rnbinom_mean_variance(200000, mean=10.0, var=40.0, seed=1)
>>> bool(abs(x.mean() - 10.0) < 0.2 and abs(x.var() - 40.0) < 2.0)
True

spacr.power_simulate.sample_count_mean_variance(n: int, mean: float, var: float | None = None, *, rng: numpy.random.Generator | None = None, seed: int | None = None) → numpy.ndarray

Draw n non-negative integer counts with the requested mean and variance.

Dispatches on the dispersion, because no single two-parameter count family covers the whole range:

dispersion	distribution
var > mean	NegativeBinomial(mean**2/(var-mean), mean/var)
var == mean	Poisson(mean)
0 <= var < mean	Binomial(n, mean/n) with n the nearest integer to mean**2 / (mean - var)

This replaces upstream’s COMPoissonReg::rcmp dispersion parameter nu, which is called but not declared in the R package’s DESCRIPTION and has no maintained Python equivalent inside spaCR’s dependency set. The third and higher moments differ from COM-Poisson; nothing upstream used them.

Where the under-dispersed case is approximate, and which way. A binomial has an integer number of trials, so it cannot hit an arbitrary (mean, var) pair. Given the rounded n, the success probability is set to mean / n so the mean is exact and the realised variance is mean * (1 - mean/n), the nearest value the family attains — within about 1/n relative of the request. The mean is the one preserved because it is the biologically meaningful quantity (cells per well) and it sets the scale of everything downstream, whereas a variance off by a percent changes only how wide the well-to-well spread is. Rounding n while keeping p = (mean - var) / mean — the naive inversion — biases the mean instead, by up to half a count, which is a systematic error in every well of a sweep.

Parameters:

• n – Number of variates.

• mean – Target mean, must be finite and > 0.

• var – Target variance. None means “Poisson”, i.e. var = mean.

• rng – Generator to draw from; mutually exclusive with seed.

• seed – Seed for a fresh generator; mutually exclusive with rng.

Returns:

Integer array of shape (n,).

Raises:

• ImpossibleMomentsError – If var is negative. A count distribution cannot have negative variance, and the caller has almost certainly passed a standard deviation or a coefficient of variation by mistake.

• ScreenDesignError – If mean is not finite and positive, n is negative, or neither/both of rng and seed given.

Example:

>>> under = sample_count_mean_variance(100000, mean=100.0, var=25.0, seed=4)
>>> bool(abs(under.mean() - 100.0) < 0.5 and abs(under.var() - 25.0) < 1.0)
True

spacr.power_simulate.simulate_imaging_plate(spot_plate: pandas.DataFrame, imaging_n_cells_per_well_mu: float, imaging_n_cells_per_well_var: float | None, class_pos_mu: float, class_pos_var: float, class_neg_mu: float, class_neg_var: float, *, imaging_split: str = 'abundance', rng: numpy.random.Generator | None = None, seed: int | None = None) → pandas.DataFrame

Simulate the imaged cells per genotype per well and the classifier’s calls.

Per well: draw the well’s total imaged cell count from sample_count_mean_variance(), split it multinomially across the genes present, then call each cell positive with a probability drawn per (gene, well) from a beta — class_pos_* for cells of a hit genotype, class_neg_* otherwise.

The gap between class_pos_mu and class_neg_mu is the signal the whole power analysis is trying to detect. The real MaxViT classifier this was fitted to sat at 0.80 / 0.12, which is a modest gap, and the point of the exercise is that a modest gap is survivable given enough wells.

imaging_split chooses how a well’s cells are divided between the genes in it:

• 'abundance' (default) weights by each gene’s library abundance, so a gene that is 20% of the well gets 20% of its cells.

• 'uniform' splits evenly, which is what the R package does — its prob = gene_in_well is a 0/1 vector, so abundance is ignored at this step. Kept for parity; it understates the imbalance between genotypes and therefore overstates power.

Parameters:

• spot_plate – Result of simulate_spot_plate(). Needs gene, well, gene_in_well and hit, plus gene_abundance when imaging_split='abundance'.

• imaging_n_cells_per_well_mu – Mean cells imaged per well, must be > 0.

• imaging_n_cells_per_well_var – Variance of that count. None means Poisson, and the emitted imaging_n_cells_per_well_var column then echoes the mean, because that is the variance the model used. Unlike the R original, variance below the mean is allowed and draws from a binomial.

• class_pos_mu – Mean probability a hit-genotype cell is called positive.

• class_pos_var – Variance of it across (gene, well); 0 is allowed and gives a deterministic classifier.

• class_neg_mu – Mean probability a non-hit cell is called positive.

• class_neg_var – Variance of it across (gene, well).

• imaging_split – 'abundance' or 'uniform', see above.

• rng – Generator to draw from; mutually exclusive with seed.

• seed – Seed for a fresh generator; mutually exclusive with rng.

Returns:

DataFrame with one row per (gene, well) and columns [gene, well, imaging_n_cells_per_well_mu, imaging_n_cells_per_well_var, imaging_n_cells_per_gene_per_well, imaging_n_cells_per_well, class_pos_mu, class_pos_var, class_neg_mu, class_neg_var, positive].

imaging_n_cells_per_well is the well total, repeated on every row of the well. It is not in the R output — R/fit_model.R reads it as the Poisson offset, R’s $ partial matching finds three columns with that prefix and returns NULL, and the fit silently loses its offset. It is emitted here because the model half needs it.

Raises:

• ScreenDesignError – If imaging_split is not one of the two supported values, the cell-count mean is not positive, or neither/both of rng and seed given.

• ImpossibleMomentsError – If a classifier (mu, var) pair is not a realisable beta — see rbeta_mean_variance().

• MalformedPlateError – If spot_plate is not one row per (gene, well) or is missing a required column.

Example:

A perfect classifier must call every hit cell positive and no other:

>>> lib = simulate_library(20, 50.0, 0.5, seed=8)
>>> spot = simulate_spot_plate(lib, 40, 2.0, 0.05, seed=9)
>>> img = simulate_imaging_plate(spot, 100.0, None,
...                              class_pos_mu=1.0, class_pos_var=0.0,
...                              class_neg_mu=0.0, class_neg_var=0.0, seed=10)
>>> merged = img.merge(spot[['gene', 'well', 'hit']], on=['gene', 'well'])
>>> bool((merged.loc[merged['hit'] == 0, 'positive'] == 0).all())
True

spacr.power_simulate.simulate_library(n_genes_in_library: int, gene_abundance_alpha: float, gene_hit_rate: float, *, rng: numpy.random.Generator | None = None, seed: int | None = None) → pandas.DataFrame

Simulate the perturbation library: which genes exist, how abundant, which are hits.

gene_abundance ~ Dirichlet(gene_abundance_alpha * 1_n) and hit_i ~ Bernoulli(gene_hit_rate), independently.

gene_abundance_alpha controls how even the library is, and it runs the opposite way from intuition: large alpha gives every gene about 1/n of the pool, alpha near zero concentrates the whole pool on one gene, and alpha = 1 corresponds to a Gini index of 0.5. The real T. gondii screen this simulator was fitted to came out at alpha = 0.6, i.e. quite skewed — which is why some genes in it never reached enough wells to be callable at all.

Parameters:

• n_genes_in_library – Number of genes, must be > 0.

• gene_abundance_alpha – Dirichlet concentration, must be > 0.

• gene_hit_rate – Probability each gene is a true hit, in [0, 1].

• rng – Generator to draw from; mutually exclusive with seed.

• seed – Seed for a fresh generator; mutually exclusive with rng.

Returns:

DataFrame with n_genes_in_library rows and columns [gene, gene_abundance_alpha, gene_hit_rate, gene_abundance, hit]. gene is a 1-based integer index (matching the R package), and the gene_abundance column sums to exactly 1.

Raises:

ScreenDesignError – If the library is empty, the concentration is not positive, the hit rate is outside [0, 1], or neither/both of rng and seed given. An empty library is refused rather than returned as an empty frame, because every downstream stage would then produce an empty frame too and the run would report “no hits found” instead of “you asked for no genes”.

Example:

>>> lib = simulate_library(500, gene_abundance_alpha=10.0,
...                        gene_hit_rate=0.1, seed=5)
>>> bool(abs(lib['gene_abundance'].sum() - 1.0) < 1e-12)
True
>>> int(lib['gene'].iloc[0]), int(lib['gene'].iloc[-1])
(1, 500)

spacr.power_simulate.simulate_screen(n_genes_in_library: int, gene_abundance_alpha: float, gene_hit_rate: float, n_wells_per_screen: int, well_abundance_factor_mu: float, well_abundance_factor_var: float, imaging_n_cells_per_well_mu: float, imaging_n_cells_per_well_var: float | None, class_pos_mu: float, class_pos_var: float, class_neg_mu: float, class_neg_var: float, sequencing_n_cells_per_well_lambda: float, pcr_factor_mu: float, pcr_factor_var: float, n_reads_per_well: float, *, sequencing_n_cells_per_well_var: float | None = None, read_depth_cv: float = 0.0, sequencing_error_rate: float = 0.0, min_cells_per_well: int = 0, imaging_split: str = 'abundance', rng: numpy.random.Generator | None = None, seed: int | None = None) → pandas.DataFrame

Run all four stages and return one tidy (well, gene) table.

The imaging and sequencing plates are both simulated from the same spot plate — which genotypes are in a well is one physical fact, observed twice. Each stage draws from its own spawned child stream, so changing the number of imaging cells does not shift the sequencing draws; that independence is what makes a parameter sweep interpretable, since otherwise every point on the sweep would differ by an unrelated re-randomisation as well as by the parameter.

Parameters:

• n_genes_in_library – See simulate_library().

• gene_abundance_alpha – See simulate_library().

• gene_hit_rate – See simulate_library().

• n_wells_per_screen – See simulate_spot_plate().

• well_abundance_factor_mu – See simulate_spot_plate().

• well_abundance_factor_var – See simulate_spot_plate().

• imaging_n_cells_per_well_mu – See simulate_imaging_plate().

• imaging_n_cells_per_well_var – See simulate_imaging_plate().

• class_pos_mu – See simulate_imaging_plate().

• class_pos_var – See simulate_imaging_plate().

• class_neg_mu – See simulate_imaging_plate().

• class_neg_var – See simulate_imaging_plate().

• sequencing_n_cells_per_well_lambda – See simulate_sequencing_plate().

• pcr_factor_mu – See simulate_sequencing_plate().

• pcr_factor_var – See simulate_sequencing_plate().

• n_reads_per_well – See simulate_sequencing_plate().

• sequencing_n_cells_per_well_var – See simulate_sequencing_plate().

• read_depth_cv – See simulate_sequencing_plate().

• sequencing_error_rate – See misassign_reads(). Off by default; turning it on always lowers power, because a mis-assigned read moves the covariate toward uniform and attenuates the coefficient.

• min_cells_per_well – See drop_low_cell_wells(). Off by default. Turning it on trades wells for well quality, and which way that comes out depends on how long the thin tail is – which is the thing worth simulating rather than arguing about.

• imaging_split – See simulate_imaging_plate().

• rng – Generator to draw from; mutually exclusive with seed.

• seed – Seed for a fresh generator; mutually exclusive with rng.

Returns:

DataFrame with one row per (gene, well), carrying every column of all four stages: the ground truth hit and gene_abundance, the observed positive and imaging_n_cells_per_well, and the observed n_reads_per_gene_per_well. This is the frame the model half consumes.

Raises:

• ScreenDesignError – Propagated from any stage; see each.

• ImpossibleMomentsError – Propagated from any stage; see each.

Example:

>>> screen = simulate_screen(
...     n_genes_in_library=40, gene_abundance_alpha=20.0, gene_hit_rate=0.1,
...     n_wells_per_screen=24, well_abundance_factor_mu=4.0,
...     well_abundance_factor_var=0.5,
...     imaging_n_cells_per_well_mu=120.0, imaging_n_cells_per_well_var=8000.0,
...     class_pos_mu=0.8, class_pos_var=0.01,
...     class_neg_mu=0.12, class_neg_var=0.005,
...     sequencing_n_cells_per_well_lambda=1000.0,
...     pcr_factor_mu=2.0, pcr_factor_var=1.0, n_reads_per_well=30000,
...     seed=14)
>>> len(screen)
960
>>> bool(screen['positive'].sum() > 0)
True

spacr.power_simulate.simulate_sequencing_plate(spot_plate: pandas.DataFrame, sequencing_n_cells_per_well_lambda: float, pcr_factor_mu: float, pcr_factor_var: float, n_reads_per_well: float, *, sequencing_n_cells_per_well_var: float | None = None, read_depth_cv: float = 0.0, sequencing_error_rate: float = 0.0, rng: numpy.random.Generator | None = None, seed: int | None = None) → pandas.DataFrame

Simulate the barcode read counts for each genotype in each well.

Four steps per well, the last of them optional:

1. Cells contributing DNA: gene_in_well * Count(lambda, var). This is normally far larger than the imaged count, because sequencing sees the whole well and the microscope sees a few fields.

2. Amplification: one lognormal PCR factor per well, pcr_factor ~ LogNormal(meanlog=pcr_factor_mu, sdlog=sqrt(pcr_factor_var)). Note both parameters are on the log scale despite their names; a pcr_factor_mu of 2.0 is a median amplification of exp(2) ~ 7.4. The factor is per well, not per gene, because a well is one PCR reaction — which is exactly why read counts within a well are not independent.

3. Sequencing: reads are drawn from the amplified barcode pool without replacement (multivariate hypergeometric), because a finite flow cell reading a finite library is sampling without replacement, and modelling it as multinomial overstates how much independent information deep wells carry.

4. Mis-assignment, when sequencing_error_rate > 0: a fraction of the reads is credited to the wrong gene. Off by default, because it is not in the R and a silently different baseline would make every number already quoted from this module wrong. See misassign_reads() for what it models and why it always costs power.

Two departures from the R original, both deliberate. Upstream computes n_reads_per_well = round(n_reads_total / nrow(well_data)) inside a group_by(well), so nrow is the library size: with 452 genes and n_reads_total = 128318 that is 284 reads per well, against a real screen near 3e4. The parameter is documented as “total reads in the screen” in one vignette and “geometric mean of well reads” in another. This function takes n_reads_per_well unambiguously and reports n_reads_total as the derived sum. Second, upstream’s guard k = min(n_cells_in_well * pcr_factor, n_reads_per_well) is computed from round(sum(cells) * pcr) while the urn is round(cells * pcr) element-wise; the two differ by rounding and the draw can ask for more balls than the urn holds. The guard here is computed from the urn itself.

Parameters:

• spot_plate – Result of simulate_spot_plate(); needs gene, well and gene_in_well.

• sequencing_n_cells_per_well_lambda – Mean cells per gene per well contributing DNA, must be > 0.

• pcr_factor_mu – Log-scale mean of the per-well amplification factor.

• pcr_factor_var – Log-scale variance of it, must be >= 0.

• n_reads_per_well – Target reads per well, must be >= 0.

• sequencing_n_cells_per_well_var – Variance of the per-gene cell count; None means Poisson, and the emitted column then echoes the lambda, because that is the variance the model used.

• read_depth_cv – Coefficient of variation of read depth between wells. 0.0 gives every well exactly n_reads_per_well; real screens are far from uniform, and shallow wells are where hits go to die.

• sequencing_error_rate – probability a read is credited to the wrong gene, in [0, 1]. 0.0 (the default) is the R behaviour; DEFAULT_SEQUENCING_ERROR_RATE is a realistic figure.

• rng – Generator to draw from; mutually exclusive with seed.

• seed – Seed for a fresh generator; mutually exclusive with rng.

Returns:

DataFrame with one row per (gene, well) and columns [gene, well, sequencing_n_cells_per_well_lambda, sequencing_n_cells_per_well_var, n_reads_per_well, n_reads_total, pcr_factor, sequencing_n_cells_per_gene_per_well, n_barcodes_per_genes_per_well, sequencing_error_rate, n_reads_true_per_gene_per_well, n_reads_per_gene_per_well]. n_barcodes_per_genes_per_well keeps the R package’s misspelling so the two are joinable by name.

n_reads_per_gene_per_well is what the analyst sees, i.e. after mis-assignment, because that is what the model half must consume for the dilution to be in the answer. n_reads_true_per_gene_per_well is the same quantity before it, so a test can plant the truth and measure how far the observation moved. With the error rate at zero the two are identical.

Raises:

• ScreenDesignError – If a parameter is out of range, or neither/both of rng and seed given.

• SequencingScaleError – If a well’s amplified barcode pool exceeds MAX_HYPERGEOMETRIC_URN, above which numpy’s exact sampler loses precision. Almost always means pcr_factor_mu was supplied on the linear scale.

• MalformedPlateError – If spot_plate is not one row per (gene, well) or is missing a required column.

Example:

Reads never exceed the requested depth, and a gene absent from a well never gets a read:

>>> lib = simulate_library(30, 20.0, 0.1, seed=11)
>>> spot = simulate_spot_plate(lib, 12, 3.0, 0.1, seed=12)
>>> seq = simulate_sequencing_plate(spot, 500.0, 1.0, 0.2, 5000, seed=13)
>>> per_well = seq.groupby('well')['n_reads_per_gene_per_well'].sum()
>>> bool((per_well <= 5000).all())
True

spacr.power_simulate.simulate_spot_plate(gene_library: pandas.DataFrame, n_wells_per_screen: int, well_abundance_factor_mu: float, well_abundance_factor_var: float, *, rng: numpy.random.Generator | None = None, seed: int | None = None) → pandas.DataFrame

Simulate which genes landed in which wells.

well_abundance_j ~ Gamma(mean=well_abundance_factor_mu, var=well_abundance_factor_var) and gene_in_well_ij ~ Bernoulli(gene_abundance_i * well_abundance_j).

well_abundance_factor_mu is the knob that trades genes-per-well against wells-per-gene, and it is the sweep the R package cared most about: with 452 genes, mu = 4.6 gives roughly 4.6 genes per well.

Probabilities above 1 are clipped, loudly. The Bernoulli probability is a product of two independently drawn quantities and nothing constrains it to [0, 1]; at alpha = 0.6 the most abundant gene’s share times a mu = 4.6 well factor exceeds 1 routinely. R’s rbinom returns NA there and carries on. This clips to 1 and emits one AbundanceClippedWarning carrying the number of clipped cells and the largest offending probability, because a clipped run has a realised genes-per-well below the one that was requested — the answer is still usable, but not for the parameter value written on it.

Parameters:

• gene_library – Result of simulate_library(); needs at least gene and gene_abundance, and all its columns are carried through.

• n_wells_per_screen – Number of wells, must be > 0.

• well_abundance_factor_mu – Mean per-well abundance factor, must be > 0.

• well_abundance_factor_var – Variance of it, must be > 0.

• rng – Generator to draw from; mutually exclusive with seed.

• seed – Seed for a fresh generator; mutually exclusive with rng.

Returns:

DataFrame with one row per (gene, well) pair — gene-major, well varying fastest, matching tidyr::expand_grid — and columns [gene, well, <library columns>, well_abundance_factor_mu, well_abundance_factor_var, well_abundance, gene_in_well]. DataFrame.attrs['n_prob_clipped'] records the clip count; note pandas drops attrs through most merges, so the warning is the durable record.

Raises:

• ScreenDesignError – If there are no wells, the abundance moments are not positive, or neither/both of rng and seed given. Zero wells is refused because R’s 1:0 idiom silently yields two wells.

• MalformedPlateError – If gene_library is empty or is missing gene / gene_abundance.

Example:

>>> lib = simulate_library(50, 10.0, 0.1, seed=6)
>>> spot = simulate_spot_plate(lib, 8, 2.0, 0.1, seed=7)
>>> len(spot), sorted(spot['well'].unique().tolist())[:3]
(400, [1, 2, 3])