High and low MOI, both control groups. The complement control group in
either MOI and the NT-cells control group in low MOI, for all three analyses.
A low-MOI analysis needs each cell to carry at most one gRNA, which sceptre's
QC enforces and pysceptre can only partly check -- see
Low MOI.
Both resampling mechanisms are validated against R, in different ways.
Permutation p-values are checked value for value on R's own draws, which a
test-only replica of R's two permutation samplers replays, under both
control groups. R's CRT draws are not replicated, so the CRT is checked
exactly on everything deterministic and on its p-values in distribution.
Run normally, pysceptre draws with numpy rather than R's generator, so a
result from either mechanism agrees with R in distribution, not draw for
draw -- see
Validation,
items 5 and 6.
assign_grnas(): all three of sceptre's methods. No run_qc(), with
one carve-out: the calibration and power checks apply the pairwise
nonzero-count thresholds, because they build or receive their own pairs and
cannot select them otherwise. The mixture, thresholding and maximum
assignment methods are ported; cell-level and gRNA-level QC are still out of
scope, including low MOI's removal of cells with zero or two or more gRNAs.
pysceptre reports which cells those are (cells_w_zero_or_twoplus_grnas,
and a maximum assignment's field of that name) but never removes them.
Feed run_discovery_analysis pairs that have already passed QC (e.g. from a
real sceptre_object's @discovery_pairs_with_info, filtered to
pass_qc == TRUE).
assign_grnas_fishash is not sceptre's. It ports fishash 0.99.5, the
version after the preprint's 0.3.0; the two differ only where 0.3.0 would
divide 0 by 0. All four assignment methods take raw integer counts, which
the .h5mu exports this repository builds carry as the optional
grna_counts assay.
The dose test is not sceptre's. With per-cell weights
(grna_target_weights, built by dose_weights from gRNA counts) the three
analyses run a score test in which each cell carries a target with a weight
rather than yes or no. It runs with the complement control group, the CRT
and the union strategy only; without weights, or with weights all 1, every
result is sceptre's. See
Design decisions.
No formula DSL.covariate_matrix must already be a plain numeric
design matrix; there's no model.matrix()-equivalent formula parser here.
RNG is not bit-for-bit reproducible against R. sceptre seeds
boost::mt19937 with a fixed literal seed on every resampling call;
pysceptre uses numpy.random.Generator with a different algorithm and
seeding scheme. Validated against R by matching distributions (the
resulting null-statistic and p-value distributions), not exact draws --
see Validation.
B1/B2/B3 are not exposed at the top-level API, but they are no
longer fixed: run_discovery_analysis derives them from
resampling_approximation exactly as R does -- B1=499 always, then
(B2, B3) = (4999, 0) for skew_normal and (0, ceil(mult * n_pairs /
multiple_testing_alpha)) for no_approximation. B3=0 on the
skew_normal path is parity with R, which only uses B3=24999 for the
permutations mechanism -- which this package does implement, and sizes the
same way.
run_discovery_ntcells_complement, run_discovery_nt_cells and
run_low_level_test_full accept all three directly if you need to override
them.
no_approximation is expensive, and can be coarser at small scale. Its
B3 grows linearly in pair count: a 33,066-pair one-sided run needs
1,653,300 draws per target (5.2 GB), so the chunk collapses to a
single target and the run warns that it may still exhaust memory.
Conversely, at 4 pairs one-sided B3 = 200, belowB1 = 499. Both are
R's behavior, reproduced rather than corrected.
All three gRNA integration strategies."union" (the default),
"singleton" and "bonferroni", matching sceptre's own option. Under the
latter two, grna_target_cells is keyed by guide rather than by target and
grna_target_data_frame supplies the design; the pair expansion is
validated against sceptre's own on a real screen, 36,450 pairs to 515,972.
Nothing statistical differs between them -- only which cells count as
treated, and what happens to the results afterwards.
compute_power is not a sceptre path, and its limits are its
own. It estimates in closed form what a screen could have detected,
which is a different question from the three analyses above, and it is a
port of PerturbPlan (MIT)
rather than of sceptre. It is validated against that package's own output,
not against sceptre's. Within it: complement control group only, an explicit
cutoff only, per-gRNA cell counts required rather than the per-target
union, and no minimum-detectable-effect-size path. It is accurate enough to
plan a screen and triage its negatives, and not accurate enough to
settle a question about one element-gene pair -- see
Design decisions.
Two ways to build its baseline statistics, and they are not
interchangeable.baseline_expression_stats_from_fits takes the mean and
theta from the same negative-binomial fit the discovery test uses, so the
estimate and the test it predicts are on one expression scale. Prefer it.baseline_expression_stats instead computes a size-factor-normalised mean,
which is what the published comparison used; measured on one screen that
mean sits about 16 % below the scale sceptre's own model works on, which
makes the estimate conservative rather than wrong-shaped. Use it only to
reproduce those numbers.
A per-pair alternative that carries the covariates.matched_expression_stats gives each (target, gene) pair the mean that
matches the information sceptre's per-cell fit gives that target's cells,
which corrects batch-driven genes and the larger libraries of perturbed
cells. It improved agreement with simulation on two DC-TAP screens but has
not been checked on day0, so it is an option beside
baseline_expression_stats_from_fits, not a replacement -- see
Design decisions.
run_specificity_check is also not sceptre's, and has no R reference.
It counts how many discovered links exceed a background measured on tests
across chromosomes. Its reference is the analysis it was ported from, which it
reproduces on three screens; nothing external has checked it. The
fingerprint part of that analysis is not included -- see
Design decisions.
run_trans_analysis, which computes that check's background, is not
sceptre's either, and is checked against pysceptre only. It tests
background pairs with a reference discovery run's own test, and a pair
tested in full gets that run's result bit for bit, given the same CPU
class, BLAS and covariate layout; nothing compares it with R. Complement
control group, CRT, skew_normal and the union strategy only, and it needs
numba. A screened row has no p-value, so a trans result serves the
specificity check and not a p-value distribution -- see
Design decisions.
glm_solver="structured" is not sceptre's route to its GLMs. It fits
the same models by per-level sums over the covariate matrix's indicator
columns, and its fits are checked against R's and against the default
route's, but not to the same bits, so the default stays "outer", which
reproduces 0.5.0's results. Complement control group only -- see
Design decisions.