Tutorials¶
One worked example end to end, then five shorter ones for the tasks that are not a single function call. Every code block here has been run before being written down, and each one names the way it is most easily got wrong.
The three sceptre analyses themselves are in the User guide.
A discovery analysis end to end¶
Synthetic inputs, small enough to run in seconds.
This walks through building every input run_discovery_analysis needs from
scratch, running it, and reading the output. It uses small synthetic data
(no real dataset required) so you can copy-paste and run this directly.
1. Install¶
pip install -e ".[dev,fast]"
2. Build the inputs¶
pysceptre needs five things: a response (gene expression) matrix, gene
IDs, a covariate (design) matrix, a mapping from gRNA target to treated
cells, and the list of (gene, target) pairs to test.
import numpy as np
import pandas as pd
rng = np.random.default_rng(0)
n_cells = 5_000
n_genes = 8
n_targets = 15
# --- covariate_matrix: (n_cells, p) numeric design matrix ---
# In a real analysis this comes from whatever produced your sceptre_object's
# @covariate_matrix (e.g. R's model.matrix() on log(n_umis), batch, etc.) --
# pysceptre does not parse a formula DSL, so build/export this matrix
# yourself. Here: intercept + two continuous covariates.
covariate_matrix = np.column_stack([
np.ones(n_cells),
rng.normal(size=n_cells), # e.g. log(response_n_umis)
rng.normal(size=n_cells), # e.g. log(grna_n_umis + 1)
])
# --- response_matrix: (n_genes, n_cells), one row per gene ---
# A genome-wide matrix is fine: only genes appearing in `pairs` are fitted.
# Untested rows cost storage, not compute (see README).
gene_ids = [f"gene_{i}" for i in range(n_genes)]
beta = rng.normal(loc=[1.0, 0.2, -0.1], scale=0.1, size=(n_genes, 3))
theta_true = rng.uniform(5, 50, size=n_genes) # NB dispersion
mu_true = np.exp(covariate_matrix @ beta.T) # (n_cells, n_genes)
response_matrix = rng.negative_binomial(
n=theta_true[None, :], p=theta_true[None, :] / (theta_true[None, :] + mu_true)
).T.astype(float) # (n_genes, n_cells)
# --- grna_target_cells: dict[target_id -> 0-based treated-cell indices] ---
# "union" convention: one entry per gRNA *target*, not per individual gRNA.
grna_target_cells = {
f"target_{t}": rng.choice(n_cells, size=rng.integers(50, 400), replace=False)
for t in range(n_targets)
}
# --- pairs: DataFrame['response_id', 'grna_target'] ---
# The QC-passed pairs to test. In a real analysis this comes from an
# upstream run_qc() step -- here, every gene x every target.
pairs = pd.DataFrame([
{"response_id": g, "grna_target": t}
for g in gene_ids for t in grna_target_cells
])
3. Run discovery analysis¶
from pysceptre.pipeline.api import run_discovery_analysis
result = run_discovery_analysis(
response_matrix=response_matrix,
gene_ids=gene_ids,
covariate_matrix=covariate_matrix,
grna_target_cells=grna_target_cells,
pairs=pairs,
side="both", # "left" for expected-repression screens (e.g. CRISPRi)
seed=0,
)
print(result.head())
Output columns: 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 (see README.md
for what each means, in particular stage, which tells you whether a pair's
p-value came from the initial empirical draws or a skew-normal tail fit).
4. Inject a real effect (sanity check)¶
A good way to build confidence in a new analysis is to plant a known signal and confirm it comes back significant. Make one target's treated cells have systematically lower expression of one gene (simulating, e.g., CRISPRi knockdown of an enhancer):
strong_target = "target_0"
strong_gene = "gene_0"
trt = grna_target_cells[strong_target]
response_matrix[gene_ids.index(strong_gene), trt] = rng.negative_binomial(
n=theta_true[0], p=theta_true[0] / (theta_true[0] + mu_true[trt, 0] * 0.3)
) # ~70% knockdown in treated cells
result = run_discovery_analysis(
response_matrix=response_matrix, gene_ids=gene_ids,
covariate_matrix=covariate_matrix, grna_target_cells=grna_target_cells,
pairs=pairs, seed=0, # side defaults to "both"; "left" is also reasonable here,
) # since the planted effect is a knockdown
hit = result[(result.response_id == strong_gene) & (result.grna_target == strong_target)]
print(hit) # expect a very small p_value, pct_change_es well below 0, stage == 2
5. Tuning for your dataset's scale¶
target_chunk_size(default200) controls how many gRNA targets' logistic fits + CRT draws are held in memory simultaneously. If you're running out of memory at real dataset scale (thousands of targets, 500k+ cells), lower this. If you have memory to spare, raising it gives a modest speed gain from better batching. You should not normally need to set it --chunk_memory_gbbounds it automatically.- Install the
fastextra (pip install -e ".[fast]") to getnumba-accelerated CRT sampling. This is close to a strict improvement with no downside at real dataset scale; without it,pysceptrestill works correctly, just slower. n_jobs(default1) parallelizes the per-pair tests, which are about 80% of the runtime. A negative value uses every core. Results are bit-identical at any worker count -- only the genes inside an already-drawn target chunk are distributed, so the resampling draws are made in the same order regardless. Linux uses processes; other platforms use threads and hit a lower ceiling, becauseforkafter macOS's Accelerate BLAS can deadlock.- You no longer need to pre-filter
response_matrixto the tested genes. The engine fits only genes that appear inpairs, so passing a genome-wide matrix costs memory to hold but not time to fit. Earlier versions fitted every row, which made this the single biggest real-world mistake in early testing: 38,606 genes fitted to test 244. Holding the matrix is still real memory, though -- never densify one (a 38,606-gene x 131k-cell densify needs ~38 GiB and was OOM-killed once), and for a large matrix prefer a backed reader, which serves genes from disk on demand.
Per-gRNA tests instead of per-target¶
Pass grna_integration_strategy="singleton" with the design table, and
grna_target_cells keyed by guide:
result = run_discovery_analysis(
response_matrix, gene_ids, covariate_matrix,
targeting_grna_cells, # keyed by guide, both from the export
pairs, # still (response_id, grna_target)
grna_integration_strategy="singleton",
grna_target_data_frame=design,
seed=0, # side defaults to "both", as in sceptre
)
# -> response_id, grna_id, grna_target, p_value, pct_change_es, ...
Each pair fans out to one row per guide of its target, and the result carries
both the guide tested and the target it belongs to, sorted by p-value with
missing values last. That is R's shape and R's ordering. "bonferroni" takes
the same per-guide tests and collapses them back to one row per target: the
smallest p-value among the guides that passed QC, multiplied by how many
passed and capped at 1, carrying that guide's other numbers.
Nothing statistical changes between the strategies. The same CRT, the same
score statistic, the same escalation. Only the treated cell set differs, which
is exactly how sceptre works: it tests whatever it calls a grna_group, and
the strategy only decides what a group is.
Three things worth knowing.
It is roughly a guides-per-target multiple of the work. On one real screen 36,450 pairs expand to 515,972, so a singleton run is about 14 times a union run, not a flag flip. A guide shared by two overlapping targets is tested once and reported under both, since the gene and the treated cells are identical.
The multiple-testing burden is per guide. A threshold derived from a union run is the wrong one; derive it from this run's own p-values.
QC is per guide too. A guide's cell count is far below its target's, so pairs that passed QC at target resolution can fail at guide resolution. The discovery path is fed pairs already judged testable and does not recheck, so if that distinction matters, recompute the pairwise counts at guide resolution before you trust a null.
One thing you may be warned about. A design table that lists the same guide
twice against the same target has its rows kept, because that is what R does,
and you will see a warning saying so. It is inert under "union", where a
target's cells are a set union over its guides, but under "singleton" the
guide appears twice in the result and under "bonferroni" it is counted twice
in the correction, inflating the corrected p-value. Pass
drop_duplicate_design_rows=True to collapse them, at the cost of no longer
matching R's row count.
A design matrix with interactions¶
covariate_matrix is a plain numeric design matrix and the engine imposes
nothing on its columns, so anything you can write as columns works. There is
no formula DSL, which means no model.matrix() to fight and also no
model.matrix() to catch your mistakes.
timepoint = ... # (n_cells,)
batch = ... # (n_cells,) 0/1
covariate_matrix = np.column_stack([
np.ones(n_cells), # intercept
np.log(response_n_umis), # sceptre's own covariates, if you want them
timepoint,
batch,
timepoint * batch, # the interaction
])
Hold out a reference level. The one real hazard is rank deficiency, and
hand-built designs are where it happens: a factor's full dummy set beside an
intercept, a covariate repeated on two scales, or a column that is constant
within another's levels. Any of those makes the weighted least squares
singular. run_discovery_analysis refuses such a matrix before fitting
anything and names the redundant columns, so the failure is a sentence rather
than a LinAlgError from inside the solver -- but it is still your design to
fix.
Two smaller notes. Every fit scales with the number of columns, so a wide
design costs real time on a large screen. And the columns you pass are the
columns the null model conditions on, which is what makes a covariate a
control rather than a variable of interest: adding timepoint asks whether
the perturbation has an effect beyond what time explains.
Power for pairs the screen could not test¶
The pairs that fail pairwise QC carry a NaN p-value because no test was ever run on them. They are the clearest case for an analytical estimate, since a p-value has nothing to say and the question -- could this screen ever have seen this -- is exactly what the closed form answers.
power = compute_power(
untested_pairs, # the pass_qc == False rows
cells_per_grna, baseline_expression_stats,
fold_change_mean=0.85, fold_change_sd=0.13,
cutoff=cutoff, num_total_cells=n_cells,
n_nonzero_trt_thresh=7, # the thresholds the QC actually used
n_nonzero_cntrl_thresh=7,
)
Set the thresholds here, unlike everywhere else. They default to 0, which
is right for pairs that already passed QC: their power should not be
discounted for a risk that did not materialise. For a pair that was never
tested the discount is the point, and the returned qc_failure_prob column
says how much of the answer it accounts for. A pair with high power and high
qc_failure_prob is one the screen could have detected had it captured the
cells.
Reporting power alongside a discovery result¶
run_discovery_analysis returns no power column, and compute_power is not
wired into it. That is deliberate: the estimator needs inputs the discovery
path has no business knowing about, and a wrapper is where a wrong cell set or
the per-target union instead of the per-gRNA sum would creep in unseen. The
join is five steps, and each one is a place to check you meant it.
from pysceptre import compute_power, run_discovery_analysis
from pysceptre.analytical_power import (
baseline_expression_stats_from_fits,
bh_nominal_cutoff,
cells_per_grna_from_assignments,
)
from pysceptre.pipeline.discovery import fit_all_genes
ALPHA = 0.1
# 1. the analysis
result = run_discovery_analysis(
response_matrix, gene_ids, covariate_matrix, grna_target_cells, pairs,
side="both", multiple_testing_alpha=ALPHA, seed=0, # "both" is the default
)
# 2. the threshold this run actually applied, then halved. side="both" gives a
# two-sided p-value, and compute_power's default side="left" wants the
# one-sided threshold that corresponds to it. Do NOT halve if you ran the
# discovery one-sided: those p-values are already one-sided.
cutoff = bh_nominal_cutoff(result["p_value"], alpha=ALPHA) / 2
# 3. the estimator's inputs, both on sceptre's own scale
fits = fit_all_genes(response_matrix, gene_ids, covariate_matrix)
baseline = baseline_expression_stats_from_fits(covariate_matrix, fits)
cells_per_grna = cells_per_grna_from_assignments(
grna_target_data_frame, targeting_grna_cells, # both from the export
)
# 4. power for the same pairs, at a 15% knockdown
power = compute_power(
pairs, cells_per_grna, baseline,
fold_change_mean=0.85, fold_change_sd=0.13,
cutoff=cutoff, num_total_cells=covariate_matrix.shape[0],
)
# 5. one table
report = result.merge(
power[["response_id", "grna_target", "power"]],
on=["response_id", "grna_target"], how="left",
)
Four things worth knowing before you run it.
Step 3 refits every gene, and the discovery run already did. It discards
its fits rather than returning them, so this repeats the most expensive part
of the analysis. If that matters, call fit_all_genes once yourself and reuse
the result; the fit depends only on the counts and the covariate matrix, not
on the pairs.
Step 2 raises if nothing is significant, rather than returning a threshold
that would hand every pair zero power. On a screen with no discoveries, pass
cutoff yourself from a plain alpha and say so in whatever you report.
Step 2 needs a real pair list to mean anything. The threshold is the
largest p-value BH calls significant, so on a handful of pairs it is simply
whatever the smallest p-value happens to be: five pairs containing one strong
hit gave 1.2e-215, and every power estimate came back 0. That is the rule R
uses, reproduced, not a defect -- but it only behaves like a threshold at
screen scale, where day0's 34,886 pairs give 7.26e-4. On a toy example pass
cutoff explicitly.
The halving in step 2 belongs to the two-sided run, not to
compute_power. side="both" is sceptre's default and pysceptre's, and it
is the only side the discovery path is validated against R on. Its p-values
are two-sided, so the one-sided threshold compute_power(side="left") wants
is half of it. Run the discovery one-sided and its p-values are already
one-sided: use the threshold unhalved. Halving twice is the
easiest way to get a confidently wrong power estimate here, and it is
invisible in the output.
The pairs that failed QC are the interesting ones. They carry a NaN
p-value because they were never tested, and they are exactly where a power
estimate says something a p-value cannot. Pass them to compute_power too,
with n_nonzero_trt_thresh and n_nonzero_cntrl_thresh set to the
thresholds the QC used, so the estimate includes the probability the pair
would have failed QC at all.
A low-MOI screen¶
Synthetic again: 3,000 cells, each carrying exactly one of fourteen gRNAs, as sceptre's low-MOI QC leaves a screen. Five targets have two guides each, four guides are non-targeting, and one target knocks its gene down.
import numpy as np
import pandas as pd
from scipy import sparse
from pysceptre import run_calibration_check, run_discovery_analysis, run_power_check
rng = np.random.default_rng(0)
n_cells, n_genes = 3000, 8
covariate_matrix = np.column_stack([np.ones(n_cells), rng.normal(size=n_cells)])
guides = [f"t{k}_g{j}" for k in range(5) for j in range(2)] + [f"nt{j}" for j in range(4)]
guide_of_cell = rng.integers(0, len(guides), size=n_cells) # one gRNA per cell
cells = {g: np.flatnonzero(guide_of_cell == i) for i, g in enumerate(guides)}
grna_target_cells = {
f"t{k}": np.sort(np.concatenate([cells[f"t{k}_g0"], cells[f"t{k}_g1"]])) for k in range(5)
}
ntc_grna_cells = {f"nt{j}": cells[f"nt{j}"] for j in range(4)} # one entry per NT guide
mu = np.exp(1.0 + 0.3 * covariate_matrix[:, 1])
counts = rng.negative_binomial(5, 5 / (5 + mu), size=(n_genes, n_cells))
counts[0, grna_target_cells["t0"]] //= 3 # t0 knocks down g0
response_matrix = sparse.csr_matrix(counts.astype(float))
gene_ids = [f"g{i}" for i in range(n_genes)]
pairs = pd.DataFrame(
[(g, t) for g in gene_ids for t in grna_target_cells], columns=["response_id", "grna_target"]
)
result = run_discovery_analysis(
response_matrix, gene_ids, covariate_matrix, grna_target_cells, pairs,
moi="low", ntc_grna_cells=ntc_grna_cells, seed=0,
)
result.sort_values("p_value").head(3)[["response_id", "grna_target", "p_value", "pct_change_es"]]
# response_id grna_target p_value pct_change_es
# g0 t0 1.383583e-92 -70.540403
# g6 t4 9.139733e-03 7.983770
# g4 t2 1.334429e-02 7.187633
moi="low" is what changed: every target is tested against the cells carrying
a non-targeting guide, by permutation, which are sceptre's low-MOI defaults.
The calibration and power checks take the same switch:
calibration = run_calibration_check(
response_matrix, gene_ids, covariate_matrix, ntc_grna_cells,
n_calibration_pairs=len(pairs), calibration_group_size=2, moi="low", seed=0,
)
power = run_power_check(
response_matrix, gene_ids, covariate_matrix, grna_target_cells,
positive_control_pairs=pd.DataFrame({"response_id": ["g0"], "grna_target": ["t0"]}),
moi="low", ntc_grna_cells=ntc_grna_cells, seed=0,
)
# calibration: 40 pairs, 2.5% with p < 0.05
# power: g0 / t0, n_nonzero_trt 213, n_nonzero_cntrl 753 (nonzero NT cells), p 1.38e-92
Three things worth knowing.
One gRNA per cell is a precondition, and it is checked where it can be. sceptre's QC removes cells with zero or two or more gRNAs before a low-MOI analysis; pysceptre does not run that QC. A cell listed under two NT guides, or under a tested target and an NT guide, is refused with the offending target named:
ValueError: 1 tested target(s) share cells with the NT cells, e.g. 't1' (2 cells).
With control_group='nt_cells' a cell must carry at most one gRNA, ...
The cost is per pair. Against the NT cells a gene's GLM is refit for every
pair, as in R, because each pair has its own cells. The fits parallelize across
genes with n_jobs like everything else.
Both defaults can be overridden. resampling_mechanism="crt" keeps the
NT cells and resamples from a logistic fit instead, and results then stop
depending on which other targets are in the run. control_group="complement"
runs the high-MOI engine on a low-MOI screen.
Two longer tutorials are planned and not yet written: a real screen from
.h5mu through to a result frame, and the scanpy/MuData interop path. Until
they exist,
examples/scanpy_interop.py
in the repository shows the second as runnable code.