Is the screen any good?#

Before calling hits, check that the screen measured something: that the cells survived, that treatments can be told from the controls, which features are reproducible, how much plate structure is left, and whether more replicates would help. Hit calling, concentration response and mechanism retrieval each have their own tutorial; this page is the gate they sit behind.

import mantispy as mt

A real screen#

pki is the JUMP pilot’s kinase inhibitor set [Chandrasekaran et al., 2023]: fifteen compounds in U2OS cells, with doses drawn from a seven-point range (eleven compounds at three doses, four at one), across eight plates that each run the same plate map. Every treatment therefore sits on all eight plates with 32 to 64 replicate wells. The controls are DMSO (negcon); four kinase inhibitors at the top dose stand in as positive controls (poscon). The deep, repeated design is what the reproducibility and saturation sections rely on. It is a shallow dose set, so concentration-response fitting lives in Concentration response.

adata = mt.ds.pki()
mt.pp.normalize(adata, method="mad_robustize", by="Metadata_Plate", reference="negcon")
mt.pp.feature_select(adata, na_cutoff=0.0)
adata = mt.pp.subset_features(adata)

treated = adata[~adata.obs["Metadata_Control"].to_numpy()].copy()
treated.shape

(2112, 850)

The preprocessing is three steps. normalize centres each plate on its DMSO median and scales by its robust spread (mad_robustize), so plates become comparable; feature_select drops features with no spread or too many missing values; subset_features keeps the survivors. What is left is 2112 treated wells (the DMSO controls are set aside through Metadata_Control) by 850 features, and every number below is counted over those 850.

Did the cells survive?#

A profile summarises the cells in a well, so it only means something if there were cells. High-dose kinase inhibitors are cytotoxic, and a well with few cells gives a noisy median that drifts from the controls for a reason that is not biology. Check the counts before trusting any separation.

import numpy as np

counts = adata.obs["Metadata_CellCount"].to_numpy(float)
is_control = adata.obs["Metadata_Control"].to_numpy(bool)
control_median = np.median(counts[is_control])
{
    "control median cells/well": int(control_median),
    "treated median cells/well": int(np.median(counts[~is_control])),
    "treated wells below half the control median": int((counts[~is_control] < 0.5 * control_median).sum()),
    "of treated wells": int((~is_control).sum()),
}

{'control median cells/well': 1679,
 'treated median cells/well': 1664,
 'treated wells below half the control median': 578,
 'of treated wells': 2112}

mt.pl.cell_counts(adata, groupby="Metadata_Control", labels={True: "control", False: "treated"});

The typical treated well holds about as many cells as a control (1664 versus 1679), but 578 of the 2112 treated wells, more than a quarter, fall below half the control median: these are the high doses of the cytotoxic compounds, where a well’s median is built from too few cells to trust. That does not sink the screen, but a strong “hit” at a high dose could be cell loss rather than phenotype. Hits, effects and cell loss follows the flagged compounds down the dose ladder with tl.cytotoxicity, and Plate artifacts regresses cell count out of the features. The rule: if the strongest hits are also the emptiest wells, suspect cytotoxicity, not phenotype.

Can a treatment be told from the controls?#

Mean average precision (mAP) asks a retrieval question: when profiles are ranked by similarity to one profile, do its replicates rank above the profiles that should not match? It is rank-based, so it needs no correlation threshold, and its permutation null gives every group a p-value.

tl.map runs in three modes, each a different question. mode="activity" ranks every treatment against the pooled negative controls, used here; mode="cross_plate" asks whether a treatment’s replicates find each other across plates; mode="replicability" scores replicate retrieval within the screen, used later. It wraps copairs [Kalinin et al., 2025], the Broad’s implementation, so the numbers are comparable with published ones.

Compare mAP with its null, not with an intuition about scores between zero and one. Unlike accuracy it has no fixed chance level: a random ranking’s score depends on how many positives each query has, so 0.14 can be a strong result on one screen and noise on another. null_size and the qvalue column (the permutation p-value, corrected across the groups tested) account for this.

mt.tl.map(adata, mode="activity")
activity = adata.uns["mantispy"]["map"]
print(f"{int(activity['is_significant'].sum())} of {len(activity)} treatments told apart from DMSO at q < 0.05")
activity.sort_values("score")[["Metadata_Perturbation", "score", "qvalue"]]

36 of 37 treatments told apart from DMSO at q < 0.05
Metadata_Perturbation score qvalue
31 BRD-U00086681-001-01-9@0.04 0.455978 0.071393
25 BRD-U00086679-001-01-9@0.2 0.472622 0.035866
21 BRD-U00086677-001-01-9@0.04 0.480976 0.027800
19 BRD-U00086677-001-01-9@0.004 0.482819 0.027312
7 BRD-U00086673-001-01-9@0.2 0.491324 0.017713
26 BRD-U00086679-001-01-9@0.4 0.493617 0.015955
22 BRD-U00086678-001-01-9@0.004 0.493689 0.015955
20 BRD-U00086677-001-01-9@0.01 0.499254 0.013072
32 BRD-U00086681-001-01-9@0.2 0.500012 0.012885
28 BRD-U00086680-001-01-9@0.04 0.504908 0.009778
5 BRD-U00086672-001-01-9@0.4 0.509000 0.007947
33 BRD-U00086681-001-01-9@0.4 0.515377 0.006119
4 BRD-U00086672-001-01-9@0.2 0.523392 0.004292
23 BRD-U00086678-001-01-9@0.01 0.525536 0.004292
8 BRD-U00086673-001-01-9@0.4 0.525881 0.004292
2 BRD-K40109029-001-06-9@2.0 0.546332 0.000185
24 BRD-U00086678-001-01-9@0.04 0.553990 0.001009
27 BRD-U00086679-001-01-9@1.0 0.594187 0.000352
6 BRD-U00086672-001-01-9@1.0 0.604875 0.000185
9 BRD-U00086673-001-01-9@1.0 0.610627 0.000185
10 BRD-U00086674-001-01-9@0.04 0.642509 0.000185
34 BRD-U00086682-001-01-9@0.04 0.643623 0.000185
29 BRD-U00086680-001-01-9@0.2 0.651759 0.000185
13 BRD-U00086675-001-01-9@0.04 0.694040 0.000185
16 BRD-U00086676-001-01-9@0.004 0.699955 0.000185
30 BRD-U00086680-001-01-9@0.4 0.729262 0.000185
17 BRD-U00086676-001-01-9@0.01 0.752951 0.000185
18 BRD-U00086676-001-01-9@0.04 0.755786 0.000185
1 BRD-K15819326-001-01-6@2.0 0.755816 0.000185
14 BRD-U00086675-001-01-9@0.2 0.797446 0.000185
11 BRD-U00086674-001-01-9@0.2 0.803030 0.000185
15 BRD-U00086675-001-01-9@0.4 0.825711 0.000185
12 BRD-U00086674-001-01-9@0.4 0.859870 0.000185
35 BRD-U00086682-001-01-9@0.2 0.868733 0.000185
0 BRD-K15179513-001-03-4@2.0 0.915201 0.000185
3 BRD-K95785537-001-26-9@2.0 0.943041 0.000185
36 BRD-U00086682-001-01-9@0.4 0.972519 0.000185 
mt.pl.map(adata);

36 of the 37 treatments separate from DMSO at q < 0.05; the one that does not is a low dose, where the phenotype is weakest. Treatment names read <compound>@<dose>, so BRD-U00086677-001-01-9@0.04 is that compound at 0.04 uM and the same compound recurs at each of its doses. There is no DMSO row: in mode="activity" the negative controls are the reference every treatment is scored against, so they are not scored themselves.

As a rule of thumb a working screen separates most of its treatments; if only a handful cleared the threshold you would suspect the assay, not the compounds. The classic single-number assay-quality metrics, Z’-factor and SSMD, ask the same question from a designated positive control and are univariate; pki’s positive controls are four kinase inhibitors of uneven strength rather than a canonical reference, so this screen leans on the multivariate mAP-against-null instead. percent_replicating and grit are older replicate-retrieval readouts, kept because many published numbers use them.

Which features can I trust?#

A feature with high variance can still be useless. If that variance lies within replicates of the same treatment rather than between treatments, it is noise. A variance filter cannot tell the difference because it ignores the replicate structure.

pp.feature_reproducibility computes the intraclass correlation (ICC), the share of a feature’s variance that lies between replicate groups.

mt.pp.feature_reproducibility(treated, by="Metadata_Perturbation")
icc = treated.var["icc"]
{
    "median": round(float(icc.median()), 3),
    "lower quartile": round(float(icc.quantile(0.25)), 3),
    "upper quartile": round(float(icc.quantile(0.75)), 3),
    "features above 0.2": int((icc > 0.2).sum()),
    "of": int(treated.n_vars),
}

{'median': 0.233,
 'lower quartile': 0.136,
 'upper quartile': 0.34,
 'features above 0.2': 495,
 'of': 850}

Does selecting on ICC help? The answer differs between two retrieval tasks.

Replicate retrieval, whether a treatment’s replicates find each other, is what mAP measures.

def replicate_map(profiles):
    """Mean average precision of replicate retrieval."""
    scratch = profiles.copy()
    mt.tl.map(scratch, mode="replicability", null_size=500, seed=0)
    return float(scratch.uns["mantispy"]["map"]["score"].mean())


{
    "all features": round(replicate_map(treated), 3),
    "ICC > 0.2": round(replicate_map(treated[:, (icc > 0.2).to_numpy()].copy()), 3),
    "ICC > 0.4": round(replicate_map(treated[:, (icc > 0.4).to_numpy()].copy()), 3),
}

{'all features': 0.238, 'ICC > 0.2': 0.247, 'ICC > 0.4': 0.26}

Replicate retrieval improves with ICC, from 0.238 with every feature to 0.260 above an ICC of 0.4. This gain is partly circular: ICC is computed on the same replicate groups that mAP then scores, so it favors features that agree within them.

Mechanism retrieval, whether a compound finds a different compound with the same mechanism, is the task from Mechanism of action. On this screen it cannot be measured at all.

The next cell returns 0.00 on pki by design, not because anything failed; the section after it explains why.

def not_same_compound(profiles):
    """Not-same-compound mechanism retrieval."""
    signatures = mt.tl.consensus(profiles.copy(), method="median", min_replicates=1)
    signatures = signatures[signatures.obs["Metadata_MOA"].notna().to_numpy()].copy()
    mt.tl.nn_moa_classify(signatures, scheme="nsc")
    return float(signatures.uns["mantispy"]["moa"]["accuracy"])


{
    "all features": round(not_same_compound(treated), 3),
    "ICC > 0.2": round(not_same_compound(treated[:, (icc > 0.2).to_numpy()].copy()), 3),
    "ICC > 0.4": round(not_same_compound(treated[:, (icc > 0.4).to_numpy()].copy()), 3),
}

{'all features': 0.0, 'ICC > 0.2': 0.0, 'ICC > 0.4': 0.0}

This is a property of pki, not a failed screen. Reproducible and generalizable are different things, and pki shows the limit by example: mechanism retrieval scores zero with or without the ICC filter, because the task needs at least two compounds that share a mechanism and pki has almost none. Thirteen of its fourteen annotated mechanisms are a single compound measured at several doses, so not-same-compound retrieval has nothing to find. A screen built to measure mechanism transfer, such as BBBC021 in Mechanism of action, is where ICC can be watched to help or hurt that task.

Use ICC when the question is reproducibility. For any other task, measure its effect first, and choose the cutoff from the ICC distribution rather than a default (the 0.2 and 0.4 above are illustrative).

How much of this is batch effect?#

pp.feature_batch_sensitivity tests each feature for dependence on the plate with a Kruskal-Wallis test and corrects across features. Run it where a plate-to-plate difference can only be technical: on the DMSO controls, which are biologically identical, and on the treated wells, which carry the same plate map on every plate so nothing forces them to differ between plates.

raw = mt.ds.pki()[:, adata.var_names].copy()
controls = raw[raw.obs["Metadata_Control"].to_numpy()].copy()
mt.pp.feature_batch_sensitivity(controls, batch_key="Metadata_Plate")

treated_normalized = adata[~adata.obs["Metadata_Control"].to_numpy()].copy()
mt.pp.feature_batch_sensitivity(treated_normalized, batch_key="Metadata_Plate")
{
    "controls, before normalization": int(controls.var["batch_sensitive"].sum()),
    "treated, after normalization": int(treated_normalized.var["batch_sensitive"].sum()),
    "of": int(controls.n_vars),
}

{'controls, before normalization': 850,
 'treated, after normalization': 847,
 'of': 850}

Before normalization every one of the 850 features separates the plates even among the DMSO wells, which are biologically identical. With about a hundred control wells per plate a Kruskal-Wallis test calls almost anything, so read 850 of 850 as “the plate signal is everywhere”, not as a count to track. After normalization the controls are centred, so a location test no longer separates them; the honest check is the treated wells, where the same plate map runs on every plate and nothing forces agreement. There the count barely falls, from 850 to 847, not because normalization failed but because with a few hundred wells per plate the test saturates and measures power, not magnitude. For magnitude, and for whether a correction actually helped, score the plates with mt.metrics.evaluate_integration; Plate artifacts does this and corrects what remains.

How many replicates do I need?#

tl.replicate_saturation subsamples replicates and reports how much the signature still changes at each depth. There are two metrics, and the number of replicates in the screen decides which one applies:

  • "signature_stability" (the default) correlates two disjoint subsets. It is unbiased and needs 2 * depth replicates.

  • "convergence" correlates a subset with the group’s full signature. It is defined up to one less than the group size and is optimistic by construction.

pki has 32 to 64 replicates per treatment, deep enough for the unbiased default to trace a real curve rather than the single point a three-replicate screen like BBBC021 allows.

mt.tl.replicate_saturation(treated, n_draws=10)
treated.uns["mantispy"]["replicate_saturation"].round(3)
n_replicates mean std n_draws
0 1 0.173 0.039 10
1 2 0.316 0.037 10
2 3 0.344 0.038 10
3 4 0.450 0.048 10
4 5 0.451 0.027 10
5 6 0.525 0.021 10
6 7 0.511 0.019 10
7 8 0.597 0.026 10
8 9 0.586 0.027 10
9 10 0.636 0.023 10
10 11 0.633 0.020 10
11 12 0.666 0.031 10
12 13 0.667 0.016 10
13 14 0.697 0.021 10
14 15 0.700 0.019 10
15 16 0.712 0.024 10
16 17 0.715 0.023 10
17 18 0.732 0.016 10
18 19 0.740 0.024 10
19 20 0.753 0.022 10
20 21 0.746 0.013 10
21 22 0.767 0.019 10
22 23 0.775 0.023 10
23 24 0.781 0.019 10
24 25 0.789 0.009 10
25 26 0.799 0.019 10
26 27 0.808 0.012 10
27 28 0.811 0.014 10
28 29 0.807 0.011 10
29 30 0.815 0.014 10
30 31 0.825 0.016 10
31 32 0.819 0.012 10 
mt.pl.replicate_saturation(treated);

From one replicate to two, unbiased signature agreement rises from 0.17 to 0.32, and it keeps climbing with depth: about 0.71 at sixteen replicates and 0.82 at thirty-two, where the curve has nearly flattened. The gains shrink after roughly fifteen to twenty replicates, which is where adding more wells stops buying much reproducibility. As a planning rule, read off the depth where the curve flattens: here it is still creeping up at thirty-two, so a screen that has to tell mechanisms apart would not be over-replicated at this depth. These values stay well below the near-perfect ones the optimistic convergence metric reports, because signature_stability scores two independent halves against each other rather than a subset against the full signature.

Summary#

A quick read on whether a screen is worth analysing:

  • Cells survived? If the strongest hits are also the emptiest wells, suspect cytotoxicity, not phenotype.

  • Activity? Most treatments should clear mAP against the controls at q < 0.05; compare the score with its null, not with an intuition about zero to one. Few clearing it points at the assay.

  • Features trustworthy? ICC separates between-treatment signal from within-replicate noise. Reproducible and generalizable are different properties, so measure the task you care about; pki is deeply reproducible yet cannot test mechanism transfer, because it has one compound per mechanism.

  • Batch? Plates differ even where they ran the same treatments. Test it on the controls and the treated wells, and correct if the structure survives normalization.

  • Enough replicates? A saturation curve still rising at your maximum depth says the next screen needs more wells per treatment.