{ "cells": [ { "cell_type": "markdown", "id": "db5ef5e6", "metadata": {}, "source": [ "# Hit callers on well profiles\n", "\n", "{func}`~mantispy.tl.hit_calling` scores a treatment against the negative controls. On well profiles, where a\n", "treatment has only its replicate wells, two of its methods apply, and this page runs both on two compound screens:\n", "\n", "- `\"mahalanobis\"`, the default: the median distance of the treatment's wells from the control centroid, in the\n", " controls' covariance.\n", "- `\"mp_value\"`: the squared Mahalanobis distance between the treatment's centroid and the controls', in the average\n", " of the two covariances. This is the multidimensional perturbation value (Hutz et al., *J Biomol Screen* 2013,\n", " [doi:10.1177/1087057112469257](https://doi.org/10.1177/1087057112469257)), and unlike the default it counts the\n", " treatment's own spread: a treatment whose replicate wells disagree is measured in a wider covariance.\n", "\n", "Three wells are too few to estimate a treatment's covariance on its own, but enough to add to the controls': the\n", "average is positive definite whenever the controls' covariance is, so `\"mp_value\"` needs only two wells per\n", "treatment. `\"hellinger\"` compares the two covariances themselves and needs at least twice as many rows as features\n", "on each side, so it does not apply here; the [cell and guide comparison](../genetics/hit_calling_methods.ipynb)\n", "covers it." ] }, { "cell_type": "code", "execution_count": 1, "id": "5fe51262", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T15:23:58.196824Z", "iopub.status.busy": "2026-11-11T15:23:58.196454Z", "iopub.status.idle": "2026-11-11T15:26:07.718899Z", "shell.execute_reply": "2026-11-11T15:26:07.717283Z" } }, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import scanpy as sc\n", "from scipy.stats import binom\n", "\n", "import mantispy as mt\n", "\n", "METHODS = {\"mahalanobis\": {}, \"mp_value\": {\"method\": \"mp_value\"}}\n", "\n", "\n", "def call(adata):\n", " \"\"\"Run both methods on the PCA representation and return each one's calls at q < 0.05.\"\"\"\n", " calls = {}\n", " for name, keywords in METHODS.items():\n", " mt.tl.hit_calling(adata, use_rep=\"X_pca\", n_permutations=1000, key_added=name, **keywords)\n", " calls[name] = adata.uns[\"mantispy\"][name].set_index(\"group\")[\"is_hit\"]\n", " return pd.DataFrame(calls)\n", "\n", "\n", "def pseudo_rate(adata, size, n_draws=100, seed=0):\n", " \"\"\"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.\"\"\"\n", " controls = adata[adata.obs[\"Metadata_Control_Type\"].astype(str).to_numpy() == \"negcon\"].copy()\n", " rng = np.random.default_rng(seed)\n", " called = dict.fromkeys(METHODS, 0)\n", " for draw in range(n_draws):\n", " picked = np.zeros(controls.n_obs, dtype=bool)\n", " picked[rng.choice(controls.n_obs, size, replace=False)] = True\n", " controls.obs[\"Metadata_Perturbation\"] = np.where(picked, \"pseudo\", \"DMSO\")\n", " # The reference is selected through Metadata_Control_Type, so flip it along with Metadata_Control.\n", " controls.obs[\"Metadata_Control_Type\"] = np.where(picked, \"treatment\", \"negcon\")\n", " controls.obs[\"Metadata_Control\"] = ~picked\n", " for name, keywords in METHODS.items():\n", " mt.tl.hit_calling(controls, use_rep=\"X_pca\", n_permutations=1000, seed=draw, key_added=name, **keywords)\n", " called[name] += int(controls.uns[\"mantispy\"][name].set_index(\"group\").loc[\"pseudo\", \"pvalue\"] < 0.05)\n", " # A calibrated test stays at or below this rate in 95% of runs of n_draws draws.\n", " chance = binom.ppf(0.95, n_draws, 0.05) / n_draws\n", " return pd.DataFrame(\n", " {\"rate at p < 0.05\": {name: count / n_draws for name, count in called.items()}, \"95% chance limit\": chance}\n", " ).rename_axis(f\"{n_draws} pseudo-treatments of {size} wells\")" ] }, { "cell_type": "markdown", "id": "4aa9896c", "metadata": {}, "source": [ "## BBBC021\n", "\n", "The overview's screen and preprocessing: MCF-7 cells, a compound at one concentration is a treatment, and most\n", "treatments have three replicate wells." ] }, { "cell_type": "code", "execution_count": 2, "id": "44475bf9", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T15:26:07.723835Z", "iopub.status.busy": "2026-11-11T15:26:07.723399Z", "iopub.status.idle": "2026-11-11T15:26:48.958056Z", "shell.execute_reply": "2026-11-11T15:26:48.956536Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":2: UserWarning: 2 of 467 features have no spread in at least one group of 'Metadata_Plate' among the rows selected by reference='negcon' (Cells_AreaShape_EulerNumber, Nuclei_AreaShape_EulerNumber). They cannot be scaled there, so their values in those groups are set to NaN. They are flagged in var['degenerate_scale'], and pp.feature_select drops them by default; to drop them now, adata = adata[:, ~adata.var['degenerate_scale']].copy(). A variance or outlier cut does not catch a feature that varies across a plate but is constant among its control wells; the flag does.\n", " mt.pp.normalize(bbbc, method=\"mad_robustize\", by=\"Metadata_Plate\", reference=\"negcon\", keep_raw=True)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ ":96: UserWarning: 1 selected feature(s) hold values more than outlier_cutoff=500 robust standard deviations from their own median (Cytoplasm_Texture_Gabor_CorrActin_3). They can dominate distances and PCA. Add 'drop_outliers' to operations to drop them, or check where the extreme values come from; on single cells, before aggregation, since a median profile hides a few extreme cells.\n", " func(target, **params, **extra)\n" ] }, { "data": { "text/plain": [ "wells per treatment\n", "3 96\n", "2 7\n", "Name: treatments, dtype: int64" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bbbc = mt.ds.bbbc021()\n", "mt.pp.normalize(bbbc, method=\"mad_robustize\", by=\"Metadata_Plate\", reference=\"negcon\", keep_raw=True)\n", "mt.pp.feature_select(bbbc, na_cutoff=0.0)\n", "bbbc = mt.pp.subset_features(bbbc)\n", "sc.pp.pca(bbbc, n_comps=30)\n", "\n", "treated = bbbc.obs[\"Metadata_Control_Type\"].astype(str).to_numpy() != \"negcon\"\n", "wells = bbbc.obs.loc[treated, \"Metadata_Perturbation\"].astype(str).value_counts()\n", "wells.value_counts().rename(\"treatments\").rename_axis(\"wells per treatment\")" ] }, { "cell_type": "markdown", "id": "5619def2", "metadata": {}, "source": [ "First the calibration. Control wells relabelled as a treatment of the typical size should be called at about\n", "the nominal rate." ] }, { "cell_type": "code", "execution_count": 3, "id": "782e7173", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T15:26:48.961138Z", "iopub.status.busy": "2026-11-11T15:26:48.960896Z", "iopub.status.idle": "2026-11-11T15:27:55.785902Z", "shell.execute_reply": "2026-11-11T15:27:55.783812Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
rate at p < 0.0595% chance limit
100 pseudo-treatments of 3 wells
mahalanobis0.070.09
mp_value0.020.09
\n", "
" ], "text/plain": [ " rate at p < 0.05 95% chance limit\n", "100 pseudo-treatments of 3 wells \n", "mahalanobis 0.07 0.09\n", "mp_value 0.02 0.09" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pseudo_rate(bbbc, size=int(wells.median()))" ] }, { "cell_type": "code", "execution_count": 4, "id": "3735e387", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T15:27:55.789637Z", "iopub.status.busy": "2026-11-11T15:27:55.789204Z", "iopub.status.idle": "2026-11-11T15:27:57.942394Z", "shell.execute_reply": "2026-11-11T15:27:57.940666Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
mp_valueTrue
mahalanobis
False5
True98
\n", "
" ], "text/plain": [ "mp_value True\n", "mahalanobis \n", "False 5\n", "True 98" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bbbc_calls = call(bbbc).loc[wells.index]\n", "pd.crosstab(bbbc_calls[\"mahalanobis\"], bbbc_calls[\"mp_value\"])" ] }, { "cell_type": "code", "execution_count": 5, "id": "8cae67ee", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T15:27:57.945571Z", "iopub.status.busy": "2026-11-11T15:27:57.945150Z", "iopub.status.idle": "2026-11-11T15:27:58.074153Z", "shell.execute_reply": "2026-11-11T15:27:58.072536Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
wellscells per wellcontrol_cells_per_well
bryostatin@0.321576.0878
PD-169316@10.021102.0878
PD-169316@3.021095.0878
alsterpaullone@3.02990.0878
alsterpaullone@1.021302.0878
\n", "
" ], "text/plain": [ " wells cells per well control_cells_per_well\n", "bryostatin@0.3 2 1576.0 878\n", "PD-169316@10.0 2 1102.0 878\n", "PD-169316@3.0 2 1095.0 878\n", "alsterpaullone@3.0 2 990.0 878\n", "alsterpaullone@1.0 2 1302.0 878" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cells = bbbc.obs.groupby(\"Metadata_Perturbation\", observed=True)[\"Metadata_CellCount\"].mean()\n", "only = bbbc_calls.index[bbbc_calls[\"mp_value\"] & ~bbbc_calls[\"mahalanobis\"]]\n", "pd.DataFrame({\"wells\": wells[only], \"cells per well\": cells.reindex(only).round()}).assign(\n", " control_cells_per_well=round(bbbc.obs.loc[~treated, \"Metadata_CellCount\"].mean())\n", ")" ] }, { "cell_type": "markdown", "id": "34cc51cd", "metadata": {}, "source": [ "Of 100 control-only pseudo-treatments, the default called 7 at p < 0.05 and `\"mp_value\"` called 2. Both counts are\n", "within chance: a calibrated test stays at or below 9 of 100 in at least 95% of runs.\n", "\n", "The two methods agree on 98 of the 103 treatments, and `\"mp_value\"` calls the other five as well, so it calls every\n", "treatment. All five have two wells rather than three, which leaves the default calling two of the seven two-well\n", "treatments. None of the five has lost cells: each keeps more cells per well than the controls' 878." ] }, { "cell_type": "markdown", "id": "51cba4be", "metadata": {}, "source": [ "## JUMP-Target-2\n", "\n", "{func}`~mantispy.ds.jump_target2` runs one plate map in eleven laboratories, and its compounds were chosen to cover\n", "many targets rather than to be active, so some of them should not be called. There is no list of which. As a\n", "proxy, {func}`~mantispy.tl.map` with `mode=\"activity\"` asks whether a compound's replicates retrieve each other\n", "against the controls on their own plate, and a compound it misses is not known to be inactive." ] }, { "cell_type": "code", "execution_count": 6, "id": "8cd52b5a", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T15:27:58.077061Z", "iopub.status.busy": "2026-11-11T15:27:58.076764Z", "iopub.status.idle": "2026-11-11T15:28:08.864706Z", "shell.execute_reply": "2026-11-11T15:28:08.863831Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":2: UserWarning: 111 of 3616 features have no spread in at least one group of 'Metadata_Plate' among the rows selected by reference='negcon' (Cells_AreaShape_EulerNumber, Cells_AreaShape_FormFactor, Cells_Correlation_Manders_DNA_AGP, Cells_Correlation_Manders_ER_AGP, Cells_Correlation_Manders_Mito_AGP and 106 more). They cannot be scaled there, so their values in those groups are set to NaN. They are flagged in var['degenerate_scale'], and pp.feature_select drops them by default; to drop them now, adata = adata[:, ~adata.var['degenerate_scale']].copy(). A variance or outlier cut does not catch a feature that varies across a plate but is constant among its control wells; the flag does.\n", " mt.pp.normalize(jump, method=\"mad_robustize\", by=\"Metadata_Plate\", reference=\"negcon\")\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ ":96: UserWarning: 28 selected feature(s) hold values more than outlier_cutoff=500 robust standard deviations from their own median (Cells_AreaShape_Compactness, Cells_Correlation_Manders_DNA_ER, Cells_Correlation_Manders_ER_DNA, Cells_Intensity_IntegratedIntensityEdge_Mito, Cells_RadialDistribution_FracAtD_mito_tubeness_1of20 and 23 more). They can dominate distances and PCA. Add 'drop_outliers' to operations to drop them, or check where the extreme values come from; on single cells, before aggregation, since a median profile hides a few extreme cells.\n", " func(target, **params, **extra)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "301 compounds, 14 wells each (median), 898 control wells\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
rate at p < 0.0595% chance limit
100 pseudo-treatments of 14 wells
mahalanobis0.050.09
mp_value0.090.09
\n", "
" ], "text/plain": [ " rate at p < 0.05 95% chance limit\n", "100 pseudo-treatments of 14 wells \n", "mahalanobis 0.05 0.09\n", "mp_value 0.09 0.09" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "jump = mt.ds.jump_target2()\n", "mt.pp.normalize(jump, method=\"mad_robustize\", by=\"Metadata_Plate\", reference=\"negcon\")\n", "mt.pp.feature_select(jump, na_cutoff=0.0)\n", "jump = mt.pp.subset_features(jump)\n", "sc.pp.pca(jump, n_comps=30)\n", "\n", "treated = jump.obs[\"Metadata_Control_Type\"].astype(str).to_numpy() != \"negcon\"\n", "wells = jump.obs.loc[treated, \"Metadata_Perturbation\"].astype(str).value_counts()\n", "print(f\"{wells.size} compounds, {int(wells.median())} wells each (median), {int((~treated).sum())} control wells\")\n", "pseudo_rate(jump, size=int(wells.median()))" ] }, { "cell_type": "code", "execution_count": 7, "id": "2c3b0fec", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T15:28:08.866592Z", "iopub.status.busy": "2026-11-11T15:28:08.866426Z", "iopub.status.idle": "2026-11-11T15:28:27.932238Z", "shell.execute_reply": "2026-11-11T15:28:27.930934Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
calledcalled, mAP activecalled, mAP inactivemAP active, not called
mahalanobis15968911
mp_value172691030
\n", "
" ], "text/plain": [ " called called, mAP active called, mAP inactive \\\n", "mahalanobis 159 68 91 \n", "mp_value 172 69 103 \n", "\n", " mAP active, not called \n", "mahalanobis 1 \n", "mp_value 0 " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "jump_calls = call(jump).loc[wells.index]\n", "mt.tl.map(jump, mode=\"activity\", use_rep=\"X_pca\", key_added=\"activity\")\n", "active = jump.uns[\"mantispy\"][\"activity\"].set_index(\"Metadata_Perturbation\")[\"qvalue\"].reindex(wells.index) < 0.05\n", "pd.DataFrame(\n", " {\n", " name: {\n", " \"called\": int(called.sum()),\n", " \"called, mAP active\": int((called & active).sum()),\n", " \"called, mAP inactive\": int((called & ~active).sum()),\n", " \"mAP active, not called\": int((~called & active).sum()),\n", " }\n", " for name, called in jump_calls.items()\n", " }\n", ").T" ] }, { "cell_type": "markdown", "id": "a0781a10", "metadata": {}, "source": [ "Of 100 control-only pseudo-treatments of 14 wells, the default called 5 and `\"mp_value\"` called 9. Both are within\n", "chance, though `\"mp_value\"` sits at the 95% limit. `\"mp_value\"` calls\n", "172 of the 301 compounds and the default 159. Of the 69 compounds mAP calls active, `\"mp_value\"` calls all 69 and the\n", "default 68. Of the 232 that mAP does not call, `\"mp_value\"` calls 103 and the default 91.\n", "\n", "mAP is not ground truth. A compound it misses may still have a small, real effect, so these extra calls are not known\n", "to be false. They do show that `\"mp_value\"` returns the longer list on this screen." ] }, { "cell_type": "markdown", "id": "1a0e257c", "metadata": {}, "source": [ "## Which to use\n", "\n", "`\"mahalanobis\"` stays the default, and `\"mp_value\"` is opt-in. Both methods stayed within chance on control-only\n", "pseudo-treatments of both screens. On BBBC021 the treatments `\"mp_value\"` adds are the two-well ones. On JUMP-Target-2 it\n", "calls more of the compounds mAP does not call than the default does, and nothing on this page says whether those\n", "extra calls are right.\n", "\n", "On single cells the picture changes. In the [cell and guide comparison](../genetics/hit_calling_methods.ipynb),\n", "`\"mp_value\"` called nine of the twelve olfactory-receptor negatives of a CRISPR screen, where the default called one.\n", "Whichever method you pick, relabel control wells as pseudo-treatments of your group size, as `pseudo_rate` does\n", "above, to see how often it calls a group that has no effect." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.15" } }, "nbformat": 4, "nbformat_minor": 5 }