Hit callers on well profiles#

hit_calling() scores a treatment against the negative controls. On well profiles, where a treatment has only its replicate wells, two of its methods apply, and this page runs both on two compound screens:

"mahalanobis", the default: the median distance of the treatment’s wells from the control centroid, in the controls’ covariance.
"mp_value": the squared Mahalanobis distance between the treatment’s centroid and the controls’, in the average of the two covariances. This is the multidimensional perturbation value (Hutz et al., J Biomol Screen 2013, doi:10.1177/1087057112469257), and unlike the default it counts the treatment’s own spread: a treatment whose replicate wells disagree is measured in a wider covariance.

Three wells are too few to estimate a treatment’s covariance on its own, but enough to add to the controls’: the average is positive definite whenever the controls’ covariance is, so "mp_value" needs only two wells per treatment. "hellinger" compares the two covariances themselves and needs at least twice as many rows as features on each side, so it does not apply here; the cell and guide comparison covers it.

import numpy as np
import pandas as pd
import scanpy as sc
from scipy.stats import binom

import mantispy as mt

METHODS = {"mahalanobis": {}, "mp_value": {"method": "mp_value"}}


def call(adata):
    """Run both methods on the PCA representation and return each one's calls at q < 0.05."""
    calls = {}
    for name, keywords in METHODS.items():
        mt.tl.hit_calling(adata, use_rep="X_pca", n_permutations=1000, key_added=name, **keywords)
        calls[name] = adata.uns["mantispy"][name].set_index("group")["is_hit"]
    return pd.DataFrame(calls)


def pseudo_rate(adata, size, n_draws=100, seed=0):
    """Share of control-only pseudo-treatments of ``size`` wells that each method calls at p < 0.05, beside the share a calibrated test stays under."""
    controls = adata[adata.obs["Metadata_Control_Type"].astype(str).to_numpy() == "negcon"].copy()
    rng = np.random.default_rng(seed)
    called = dict.fromkeys(METHODS, 0)
    for draw in range(n_draws):
        picked = np.zeros(controls.n_obs, dtype=bool)
        picked[rng.choice(controls.n_obs, size, replace=False)] = True
        controls.obs["Metadata_Perturbation"] = np.where(picked, "pseudo", "DMSO")
        # The reference is selected through Metadata_Control_Type, so flip it along with Metadata_Control.
        controls.obs["Metadata_Control_Type"] = np.where(picked, "treatment", "negcon")
        controls.obs["Metadata_Control"] = ~picked
        for name, keywords in METHODS.items():
            mt.tl.hit_calling(controls, use_rep="X_pca", n_permutations=1000, seed=draw, key_added=name, **keywords)
            called[name] += int(controls.uns["mantispy"][name].set_index("group").loc["pseudo", "pvalue"] < 0.05)
    # A calibrated test stays at or below this rate in 95% of runs of n_draws draws.
    chance = binom.ppf(0.95, n_draws, 0.05) / n_draws
    return pd.DataFrame(
        {"rate at p < 0.05": {name: count / n_draws for name, count in called.items()}, "95% chance limit": chance}
    ).rename_axis(f"{n_draws} pseudo-treatments of {size} wells")

BBBC021#

The overview’s screen and preprocessing: MCF-7 cells, a compound at one concentration is a treatment, and most treatments have three replicate wells.

bbbc = mt.ds.bbbc021()
mt.pp.normalize(bbbc, method="mad_robustize", by="Metadata_Plate", reference="negcon", keep_raw=True)
mt.pp.feature_select(bbbc, na_cutoff=0.0)
bbbc = mt.pp.subset_features(bbbc)
sc.pp.pca(bbbc, n_comps=30)

treated = bbbc.obs["Metadata_Control_Type"].astype(str).to_numpy() != "negcon"
wells = bbbc.obs.loc[treated, "Metadata_Perturbation"].astype(str).value_counts()
wells.value_counts().rename("treatments").rename_axis("wells per treatment")

wells per treatment
3    96
2     7
Name: treatments, dtype: int64

First the calibration. Control wells relabelled as a treatment of the typical size should be called at about the nominal rate.

pseudo_rate(bbbc, size=int(wells.median()))

	rate at p < 0.05	95% chance limit
100 pseudo-treatments of 3 wells
mahalanobis	0.07	0.09
mp_value	0.02	0.09

bbbc_calls = call(bbbc).loc[wells.index]
pd.crosstab(bbbc_calls["mahalanobis"], bbbc_calls["mp_value"])

mp_value	True
mahalanobis
False	5
True	98

cells = bbbc.obs.groupby("Metadata_Perturbation", observed=True)["Metadata_CellCount"].mean()
only = bbbc_calls.index[bbbc_calls["mp_value"] & ~bbbc_calls["mahalanobis"]]
pd.DataFrame({"wells": wells[only], "cells per well": cells.reindex(only).round()}).assign(
    control_cells_per_well=round(bbbc.obs.loc[~treated, "Metadata_CellCount"].mean())
)

	wells	cells per well	control_cells_per_well
bryostatin@0.3	2	1576.0	878
PD-169316@10.0	2	1102.0	878
PD-169316@3.0	2	1095.0	878
alsterpaullone@3.0	2	990.0	878
alsterpaullone@1.0	2	1302.0	878

Of 100 control-only pseudo-treatments, the default called 7 at p < 0.05 and "mp_value" called 2. Both counts are within chance: a calibrated test stays at or below 9 of 100 in at least 95% of runs.

The two methods agree on 98 of the 103 treatments, and "mp_value" calls the other five as well, so it calls every treatment. All five have two wells rather than three, which leaves the default calling two of the seven two-well treatments. None of the five has lost cells: each keeps more cells per well than the controls’ 878.

JUMP-Target-2#

jump_target2() runs one plate map in eleven laboratories, and its compounds were chosen to cover many targets rather than to be active, so some of them should not be called. There is no list of which. As a proxy, map() with mode="activity" asks whether a compound’s replicates retrieve each other against the controls on their own plate, and a compound it misses is not known to be inactive.

jump = mt.ds.jump_target2()
mt.pp.normalize(jump, method="mad_robustize", by="Metadata_Plate", reference="negcon")
mt.pp.feature_select(jump, na_cutoff=0.0)
jump = mt.pp.subset_features(jump)
sc.pp.pca(jump, n_comps=30)

treated = jump.obs["Metadata_Control_Type"].astype(str).to_numpy() != "negcon"
wells = jump.obs.loc[treated, "Metadata_Perturbation"].astype(str).value_counts()
print(f"{wells.size} compounds, {int(wells.median())} wells each (median), {int((~treated).sum())} control wells")
pseudo_rate(jump, size=int(wells.median()))

301 compounds, 14 wells each (median), 898 control wells

	rate at p < 0.05	95% chance limit
100 pseudo-treatments of 14 wells
mahalanobis	0.05	0.09
mp_value	0.09	0.09

jump_calls = call(jump).loc[wells.index]
mt.tl.map(jump, mode="activity", use_rep="X_pca", key_added="activity")
active = jump.uns["mantispy"]["activity"].set_index("Metadata_Perturbation")["qvalue"].reindex(wells.index) < 0.05
pd.DataFrame(
    {
        name: {
            "called": int(called.sum()),
            "called, mAP active": int((called & active).sum()),
            "called, mAP inactive": int((called & ~active).sum()),
            "mAP active, not called": int((~called & active).sum()),
        }
        for name, called in jump_calls.items()
    }
).T

	called	called, mAP active	called, mAP inactive	mAP active, not called
mahalanobis	159	68	91	1
mp_value	172	69	103	0

Of 100 control-only pseudo-treatments of 14 wells, the default called 5 and "mp_value" called 9. Both are within chance, though "mp_value" sits at the 95% limit. "mp_value" calls 172 of the 301 compounds and the default 159. Of the 69 compounds mAP calls active, "mp_value" calls all 69 and the default 68. Of the 232 that mAP does not call, "mp_value" calls 103 and the default 91.

mAP is not ground truth. A compound it misses may still have a small, real effect, so these extra calls are not known to be false. They do show that "mp_value" returns the longer list on this screen.

Which to use#

"mahalanobis" stays the default, and "mp_value" is opt-in. Both methods stayed within chance on control-only pseudo-treatments of both screens. On BBBC021 the treatments "mp_value" adds are the two-well ones. On JUMP-Target-2 it calls more of the compounds mAP does not call than the default does, and nothing on this page says whether those extra calls are right.

On single cells the picture changes. In the cell and guide comparison, "mp_value" called nine of the twelve olfactory-receptor negatives of a CRISPR screen, where the default called one. Whichever method you pick, relabel control wells as pseudo-treatments of your group size, as pseudo_rate does above, to see how often it calls a group that has no effect.

Core

Ecosystem

Hit callers on well profiles#

BBBC021#

JUMP-Target-2#

Which to use#