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 a kinase inhibitor set from the Broad Institute and the Gray lab: fifteen compounds, 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. The plates come in pairs, one pair for each of four cell lines: PATU-8988T, MDA-MB-231, U2OS and A549. Every treatment therefore sits on all eight plates with 32 to 64 replicate wells, 8 to 16 in each cell line. 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.

This page pools the four cell lines. It treats a compound’s wells in all four lines as one replicate set, which they are not: the lines differ biologically. Read the activity, batch and saturation numbers below as numbers for the pooled wells. Metadata_CellLine records the line, so subset by it to ask these questions of one line.

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_Type"].to_numpy() != "negcon"].copy()
treated.shape

(2240, 854)

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 2240 treated wells (the DMSO controls are set aside through Metadata_Control_Type) by 854 features, and every number below is counted over those 854.

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_Type"].to_numpy() == "negcon"
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': 1726,
 'treated median cells/well': 1649,
 'treated wells below half the control median': 615,
 'of treated wells': 2240}

mt.pl.cell_counts(
    adata,
    groupby="Metadata_Control_Type",
    labels={"negcon": "DMSO", "poscon": "positive control", "treatment": "treated"},
);

The typical treated well holds about as many cells as a control (1649 versus 1726), but 615 of the 2240 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.455314 0.070293
25 BRD-U00086679-001-01-9@0.2 0.472155 0.036894
21 BRD-U00086677-001-01-9@0.04 0.480275 0.025897
19 BRD-U00086677-001-01-9@0.004 0.483003 0.023939
7 BRD-U00086673-001-01-9@0.2 0.492079 0.016592
22 BRD-U00086678-001-01-9@0.004 0.492994 0.016301
26 BRD-U00086679-001-01-9@0.4 0.494026 0.015634
20 BRD-U00086677-001-01-9@0.01 0.497586 0.013442
32 BRD-U00086681-001-01-9@0.2 0.499325 0.012630
28 BRD-U00086680-001-01-9@0.04 0.505498 0.008985
5 BRD-U00086672-001-01-9@0.4 0.508142 0.007810
33 BRD-U00086681-001-01-9@0.4 0.515569 0.005834
4 BRD-U00086672-001-01-9@0.2 0.523089 0.004440
23 BRD-U00086678-001-01-9@0.01 0.524732 0.004316
8 BRD-U00086673-001-01-9@0.4 0.525179 0.004316
2 BRD-K40109029-001-06-9@2.0 0.547813 0.000185
24 BRD-U00086678-001-01-9@0.04 0.553644 0.001177
27 BRD-U00086679-001-01-9@1.0 0.594475 0.000352
6 BRD-U00086672-001-01-9@1.0 0.603985 0.000185
9 BRD-U00086673-001-01-9@1.0 0.609959 0.000185
10 BRD-U00086674-001-01-9@0.04 0.642027 0.000185
34 BRD-U00086682-001-01-9@0.04 0.644617 0.000185
29 BRD-U00086680-001-01-9@0.2 0.651402 0.000185
13 BRD-U00086675-001-01-9@0.04 0.693561 0.000185
16 BRD-U00086676-001-01-9@0.004 0.700006 0.000185
30 BRD-U00086680-001-01-9@0.4 0.729580 0.000185
17 BRD-U00086676-001-01-9@0.01 0.752346 0.000185
1 BRD-K15819326-001-01-6@2.0 0.755036 0.000185
18 BRD-U00086676-001-01-9@0.04 0.756370 0.000185
14 BRD-U00086675-001-01-9@0.2 0.795772 0.000185
11 BRD-U00086674-001-01-9@0.2 0.802644 0.000185
15 BRD-U00086675-001-01-9@0.4 0.824875 0.000185
12 BRD-U00086674-001-01-9@0.4 0.858894 0.000185
35 BRD-U00086682-001-01-9@0.2 0.868337 0.000185
0 BRD-K15179513-001-03-4@2.0 0.915037 0.000185
3 BRD-K95785537-001-26-9@2.0 0.944720 0.000185
36 BRD-U00086682-001-01-9@0.4 0.972366 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 the platemap’s mmoles_per_liter of 0.04 (millimolar, and most likely a source concentration since the final well dose is not recorded) 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.

Scored within one line, the count changes; quest #236 computes it per line.

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.

Fewer control wells, fewer calls#

mode="activity" scores each treatment against the DMSO wells, so the number of DMSO wells on a plate sets how much reference there is to compare against. The next cell keeps a few DMSO wells per plate, drops the rest, and counts the calls again.

import pandas as pd

rng = np.random.default_rng(0)
plates = adata.obs["Metadata_Plate"].to_numpy()
per_plate_total = int(min((is_control & (plates == plate)).sum() for plate in np.unique(plates)))
rows = []
for per_plate in (2, 4, 8, 16, per_plate_total):
    keep = ~is_control
    for plate in np.unique(plates):
        dmso_rows = np.flatnonzero(is_control & (plates == plate))
        keep[rng.choice(dmso_rows, min(per_plate, len(dmso_rows)), replace=False)] = True
    subset = adata[keep].copy()
    mt.tl.map(subset, mode="activity", seed=0)
    result = subset.uns["mantispy"]["map"]
    rows.append(
        {
            "DMSO wells per plate": per_plate,
            "treatments called": int(result["is_significant"].sum()),
            "median mAP": round(float(result["score"].median()), 3),
        }
    )
pd.DataFrame(rows)
DMSO wells per plate treatments called median mAP
0 2 2 0.985
1 4 3 0.969
2 8 14 0.946
3 16 21 0.894
4 104 36 0.604

With two DMSO wells per plate, 2 of the 37 treatments are called. With 8 it is 14, with 16 it is 21, and with all 104 per plate it is the 36 from above. The treatments did not change, only the reference did. The median mAP goes the other way, from 0.604 with every DMSO well to 0.985 with two, because ranking a treatment’s replicates above a handful of DMSO wells is easy. The permutation null is drawn from those same few wells, so it rises with the score and fewer treatments clear it. A plate with few control wells can only call strong effects, which is worth knowing when the plate map is designed.

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.266,
 'lower quartile': 0.166,
 'upper quartile': 0.374,
 'features above 0.2': 571,
 'of': 854}

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.252, 'ICC > 0.2': 0.261, 'ICC > 0.4': 0.275}

Replicate retrieval improves with ICC, from 0.252 with every feature to 0.275 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 on the DMSO controls and on the treated wells, which carry the same plate map on every plate. On pki a plate-to-plate difference is not only technical: the plates come in pairs from four cell lines, so part of what the test finds is the cell line.

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

treated_normalized = adata[adata.obs["Metadata_Control_Type"].to_numpy() != "negcon"].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': 854,
 'treated, after normalization': 852,
 'of': 854}

Before normalization every one of the 854 features separates the plates even among the DMSO wells. Those wells are not biologically identical across pki: each pair of plates is a different cell line, so part of this is the cell line rather than the plate. With about a hundred control wells per plate a Kruskal-Wallis test calls almost anything, so read 854 of 854 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 next check is the treated wells, where the same plate map runs on every plate, though here too the cell line differs between plate pairs. There the count barely falls, from 854 to 852, 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 pooled over the four cell lines, deep enough for the unbiased default to trace a real curve rather than the single point a three-replicate screen like BBBC021 allows. Within one line it has 8 to 16.

mt.tl.replicate_saturation(treated, n_draws=10)
treated.uns["mantispy"]["replicate_saturation"].round(3)
n_replicates mean std n_draws
0 1 0.188 0.036 10
1 2 0.350 0.051 10
2 3 0.361 0.043 10
3 4 0.462 0.035 10
4 5 0.481 0.039 10
5 6 0.554 0.044 10
6 7 0.552 0.037 10
7 8 0.605 0.035 10
8 9 0.605 0.022 10
9 10 0.666 0.024 10
10 11 0.670 0.033 10
11 12 0.683 0.027 10
12 13 0.693 0.025 10
13 14 0.707 0.017 10
14 15 0.723 0.021 10
15 16 0.744 0.019 10
16 17 0.715 0.023 10
17 18 0.732 0.016 10
18 19 0.739 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.774 0.023 10
23 24 0.781 0.020 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.015 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.19 to 0.35, and it keeps climbing with depth: about 0.75 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.

The curve above pools the four cell lines, so its deeper draws mix wells from different lines and it is not a plan for a single-line screen. Within one line the 33 non-control treatments have 16 wells (the four positive controls have 8), so the curve stops at eight. Scored within one line, the curve changes; quest #236 computes it per line.

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.