{ "cells": [ { "cell_type": "markdown", "id": "c000", "metadata": {}, "source": [ "# Is the screen any good?\n", "\n", "Before calling hits, check that the screen measured something: that perturbations can be told from the\n", "controls, which features are reproducible, how much plate structure is left, and whether more replicates\n", "would help." ] }, { "cell_type": "code", "execution_count": 1, "id": "c001", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:14:51.598245Z", "iopub.status.busy": "2026-10-21T01:14:51.597849Z", "iopub.status.idle": "2026-10-21T01:14:55.166711Z", "shell.execute_reply": "2026-10-21T01:14:55.165375Z" } }, "outputs": [], "source": [ "import mantispy as mt" ] }, { "cell_type": "markdown", "id": "c002", "metadata": {}, "source": [ "## A real screen\n", "\n", "BBBC021 {cite:p}`Ljosa_2013`: MCF-7 cells, 38 compounds at one to seven concentrations with 12 mechanisms of\n", "action, two or three replicate wells per treatment. This downloads about 10 MB." ] }, { "cell_type": "code", "execution_count": 2, "id": "c003", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:14:55.170299Z", "iopub.status.busy": "2026-10-21T01:14:55.170063Z", "iopub.status.idle": "2026-10-21T01:14:55.699855Z", "shell.execute_reply": "2026-10-21T01:14:55.699078Z" } }, "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'. Adding epsilon=1e-18 to their scale multiplies them by up to 1e18, so they dominate every distance downstream. 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(adata, method=\"mad_robustize\", by=\"Metadata_Plate\", reference=\"negcon\")\n" ] }, { "data": { "text/plain": [ "(302, 343)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata = mt.ds.bbbc021()\n", "mt.pp.normalize(adata, method=\"mad_robustize\", by=\"Metadata_Plate\", reference=\"negcon\")\n", "mt.pp.feature_select(adata, na_cutoff=0.0)\n", "adata = mt.pp.subset_features(adata)\n", "\n", "treated = adata[~adata.obs[\"Metadata_Control\"].to_numpy()].copy()\n", "treated.shape" ] }, { "cell_type": "markdown", "id": "c004", "metadata": {}, "source": [ "## Can a treatment be told from the controls?\n", "\n", "Mean average precision (mAP) asks a retrieval question: when profiles are ranked by\n", "similarity to one profile, do its replicates rank above the profiles that should not match? It\n", "is rank-based, so it needs no correlation threshold, and its permutation null gives every\n", "group a p-value.\n", "\n", "`tl.map` wraps [copairs](https://github.com/cytomining/copairs) {cite:p}`Kalinin_2025`, the Broad's implementation,\n", "so the numbers are comparable with published ones.\n", "\n", "Compare mAP with its null, not with an intuition about scores between zero and one. Unlike\n", "accuracy, it has no fixed chance level: a random ranking's score depends on how many positives\n", "each query has, so 0.14 can be a strong result on one screen and noise on another. `null_size`\n", "and the `qvalue` column account for this." ] }, { "cell_type": "code", "execution_count": 3, "id": "c005", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:14:55.704050Z", "iopub.status.busy": "2026-10-21T01:14:55.703851Z", "iopub.status.idle": "2026-10-21T01:14:55.991375Z", "shell.execute_reply": "2026-10-21T01:14:55.989948Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "80 of 103 treatments told apart from DMSO at q < 0.05\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", " \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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Metadata_Perturbationscoreqvalue
48colchicine@0.030.3142860.606939
72floxuridine@30.00.5587300.281404
73lactacystin@10.00.6111110.247481
88proteasome inhibitor I@0.10.6333330.210408
42bryostatin@0.30.6666670.148555
............
98vincristine@0.031.0000000.000156
99vincristine@0.11.0000000.000156
100vincristine@0.31.0000000.000156
101vincristine@1.01.0000000.000156
102vincristine@3.01.0000000.000156
\n", "

103 rows × 3 columns

\n", "
" ], "text/plain": [ " Metadata_Perturbation score qvalue\n", "48 colchicine@0.03 0.314286 0.606939\n", "72 floxuridine@30.0 0.558730 0.281404\n", "73 lactacystin@10.0 0.611111 0.247481\n", "88 proteasome inhibitor I@0.1 0.633333 0.210408\n", "42 bryostatin@0.3 0.666667 0.148555\n", ".. ... ... ...\n", "98 vincristine@0.03 1.000000 0.000156\n", "99 vincristine@0.1 1.000000 0.000156\n", "100 vincristine@0.3 1.000000 0.000156\n", "101 vincristine@1.0 1.000000 0.000156\n", "102 vincristine@3.0 1.000000 0.000156\n", "\n", "[103 rows x 3 columns]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mt.tl.map(adata, mode=\"activity\")\n", "activity = adata.uns[\"mantispy\"][\"map\"]\n", "print(f\"{int(activity['is_significant'].sum())} of {len(activity)} treatments told apart from DMSO at q < 0.05\")\n", "activity.sort_values(\"score\")[[\"Metadata_Perturbation\", \"score\", \"qvalue\"]]" ] }, { "cell_type": "markdown", "id": "c006", "metadata": {}, "source": [ "80 of the 103 treatments are told apart from DMSO at q < 0.05, and several retrieve their replicates perfectly. At\n", "the bottom is colchicine, tested at a single dose of 0.03 uM. There is no DMSO row: in `mode=\"activity\"` the negative\n", "controls are the reference set every treatment is scored against, so they are not scored themselves.\n", "\n", "`mode=\"cross_plate\"` asks a stricter question: replicates must retrieve each other across plates, which separates\n", "reproducible biology from a plate that happened to look uniform. `percent_replicating` and `grit` are older\n", "readouts, kept because many published numbers use them." ] }, { "cell_type": "markdown", "id": "c007", "metadata": {}, "source": [ "## Which features can I trust?\n", "\n", "A feature with high variance can still be useless. If that variance lies within replicates of\n", "the same treatment rather than between treatments, it is noise. A variance filter cannot tell\n", "the difference because it ignores the replicate structure.\n", "\n", "`pp.feature_reproducibility` computes the intraclass correlation (ICC), the share of a\n", "feature's variance that lies between replicate groups." ] }, { "cell_type": "code", "execution_count": 4, "id": "c008", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:14:55.995251Z", "iopub.status.busy": "2026-10-21T01:14:55.994853Z", "iopub.status.idle": "2026-10-21T01:14:56.023332Z", "shell.execute_reply": "2026-10-21T01:14:56.022302Z" } }, "outputs": [ { "data": { "text/plain": [ "{'median': 0.387,\n", " 'lower quartile': 0.207,\n", " 'upper quartile': 0.537,\n", " 'features above 0.2': 262,\n", " 'of': 343}" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mt.pp.feature_reproducibility(treated, by=\"Metadata_Perturbation\")\n", "icc = treated.var[\"icc\"]\n", "{\n", " \"median\": round(float(icc.median()), 3),\n", " \"lower quartile\": round(float(icc.quantile(0.25)), 3),\n", " \"upper quartile\": round(float(icc.quantile(0.75)), 3),\n", " \"features above 0.2\": int((icc > 0.2).sum()),\n", " \"of\": int(treated.n_vars),\n", "}" ] }, { "cell_type": "markdown", "id": "c009", "metadata": {}, "source": [ "Does selecting on ICC help? The answer differs between two retrieval tasks.\n", "\n", "Replicate retrieval, whether a treatment's replicates find each other, is what mAP measures." ] }, { "cell_type": "code", "execution_count": 5, "id": "c010", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:14:56.026477Z", "iopub.status.busy": "2026-10-21T01:14:56.026328Z", "iopub.status.idle": "2026-10-21T01:14:56.666843Z", "shell.execute_reply": "2026-10-21T01:14:56.665689Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":4: UserWarning: 3 profile(s) have replicates but no negative pair to rank them against, and are left out. They are in [{'Metadata_Plate': 'Week6_31641'}, {'Metadata_Plate': 'Week6_31661'}, {'Metadata_Plate': 'Week6_31681'}].\n", " mt.tl.map(scratch, mode=\"replicability\", null_size=500, seed=0)\n", ":4: UserWarning: with null_size=500 no p-value can fall below 1/501, so the correction over 102 groups calls none of them unless at least 5 reach that floor together. Raise null_size above 2040 for one group to be callable on its own.\n", " mt.tl.map(scratch, mode=\"replicability\", null_size=500, seed=0)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ ":4: UserWarning: 3 profile(s) have replicates but no negative pair to rank them against, and are left out. They are in [{'Metadata_Plate': 'Week6_31641'}, {'Metadata_Plate': 'Week6_31661'}, {'Metadata_Plate': 'Week6_31681'}].\n", " mt.tl.map(scratch, mode=\"replicability\", null_size=500, seed=0)\n", ":4: UserWarning: with null_size=500 no p-value can fall below 1/501, so the correction over 102 groups calls none of them unless at least 5 reach that floor together. Raise null_size above 2040 for one group to be callable on its own.\n", " mt.tl.map(scratch, mode=\"replicability\", null_size=500, seed=0)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ ":4: UserWarning: 3 profile(s) have replicates but no negative pair to rank them against, and are left out. They are in [{'Metadata_Plate': 'Week6_31641'}, {'Metadata_Plate': 'Week6_31661'}, {'Metadata_Plate': 'Week6_31681'}].\n", " mt.tl.map(scratch, mode=\"replicability\", null_size=500, seed=0)\n", ":4: UserWarning: with null_size=500 no p-value can fall below 1/501, so the correction over 102 groups calls none of them unless at least 5 reach that floor together. Raise null_size above 2040 for one group to be callable on its own.\n", " mt.tl.map(scratch, mode=\"replicability\", null_size=500, seed=0)\n" ] }, { "data": { "text/plain": [ "{'all features': 0.317, 'ICC > 0.2': 0.343, 'ICC > 0.4': 0.369}" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def replicate_map(profiles):\n", " \"\"\"Mean average precision of replicate retrieval.\"\"\"\n", " scratch = profiles.copy()\n", " mt.tl.map(scratch, mode=\"replicability\", null_size=500, seed=0)\n", " return float(scratch.uns[\"mantispy\"][\"map\"][\"score\"].mean())\n", "\n", "\n", "{\n", " \"all features\": round(replicate_map(treated), 3),\n", " \"ICC > 0.2\": round(replicate_map(treated[:, (icc > 0.2).to_numpy()].copy()), 3),\n", " \"ICC > 0.4\": round(replicate_map(treated[:, (icc > 0.4).to_numpy()].copy()), 3),\n", "}" ] }, { "cell_type": "markdown", "id": "c011", "metadata": {}, "source": [ "Replicate retrieval improves, from 0.317 with every feature to 0.369 above an ICC of 0.4. This gain is partly\n", "circular: ICC is computed on the same replicate groups that mAP then scores, so it favors features that agree within them.\n", "\n", "Mechanism retrieval, whether a compound finds a different compound with the same mechanism,\n", "is the task from [Mechanism of action](../compounds/mechanism_of_action.ipynb). Here the same filter does not help." ] }, { "cell_type": "code", "execution_count": 6, "id": "c012", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:14:56.670195Z", "iopub.status.busy": "2026-10-21T01:14:56.669999Z", "iopub.status.idle": "2026-10-21T01:14:56.815856Z", "shell.execute_reply": "2026-10-21T01:14:56.814747Z" } }, "outputs": [ { "data": { "text/plain": [ "{'all features': 0.777, 'ICC > 0.2': 0.767, 'ICC > 0.4': 0.757}" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def not_same_compound(profiles):\n", " \"\"\"Not-same-compound mechanism retrieval.\"\"\"\n", " signatures = mt.tl.consensus(profiles.copy(), method=\"median\", min_replicates=1)\n", " signatures = signatures[signatures.obs[\"Metadata_MOA\"].notna().to_numpy()].copy()\n", " mt.tl.nn_moa_classify(signatures, scheme=\"nsc\")\n", " return float(signatures.uns[\"mantispy\"][\"moa\"][\"accuracy\"])\n", "\n", "\n", "{\n", " \"all features\": round(not_same_compound(treated), 3),\n", " \"ICC > 0.2\": round(not_same_compound(treated[:, (icc > 0.2).to_numpy()].copy()), 3),\n", " \"ICC > 0.4\": round(not_same_compound(treated[:, (icc > 0.4).to_numpy()].copy()), 3),\n", "}" ] }, { "cell_type": "markdown", "id": "c013", "metadata": {}, "source": [ "Reproducible and generalizable are different properties. ICC favors features whose\n", "replicates agree, which is what replicate retrieval rewards. Transferring from one molecule to\n", "another with the same mechanism needs something else, and filtering on reproducibility\n", "slightly hurts it here.\n", "\n", "Use ICC when the question is reproducibility. For other tasks, measure its effect first, and\n", "choose the cutoff from the ICC distribution instead of using a default.\n", "\n", "## How much of this is batch effect?\n", "\n", "`pp.feature_batch_sensitivity` tests each feature for dependence on batch with a\n", "Kruskal-Wallis test and corrects for multiple testing across features. Run it on wells that\n", "received the same treatment on every plate, such as the controls: BBBC021's plates hold different compounds,\n", "so a test on its treated wells would mix batch with treatment." ] }, { "cell_type": "code", "execution_count": 7, "id": "c014", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:14:56.818890Z", "iopub.status.busy": "2026-10-21T01:14:56.818711Z", "iopub.status.idle": "2026-10-21T01:14:57.734815Z", "shell.execute_reply": "2026-10-21T01:14:57.733889Z" } }, "outputs": [ { "data": { "text/plain": [ "{'plates': 55,\n", " 'features plate-dependent among the DMSO wells, before normalization': 335,\n", " 'of': 343}" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "raw = mt.ds.bbbc021()[:, adata.var_names].copy()\n", "controls = raw[raw.obs[\"Metadata_Control\"].to_numpy()].copy()\n", "mt.pp.feature_batch_sensitivity(controls, batch_key=\"Metadata_Plate\")\n", "{\n", " \"plates\": int(controls.obs[\"Metadata_Plate\"].nunique()),\n", " \"features plate-dependent among the DMSO wells, before normalization\": int(controls.var[\"batch_sensitive\"].sum()),\n", " \"of\": int(controls.n_vars),\n", "}" ] }, { "cell_type": "markdown", "id": "c015", "metadata": {}, "source": [ "Before normalization, 335 of the 343 features differ between the plates' DMSO wells. Normalizing each plate\n", "against its own controls centres them, so after it the controls cannot differ; what it leaves is how features\n", "respond within a plate, which only plates holding the same treatments can show.\n", "\n", "Treat this as a warning against interpreting small effects across plates. Compare a corrected\n", "run with an uncorrected one to see whether the count drops, and evaluate corrections with\n", "`mt.metrics.evaluate_integration`.\n", "\n", "## How many replicates do I need?\n", "\n", "`tl.replicate_saturation` subsamples replicates and reports how much the signature still\n", "changes at each depth. There are two metrics, and the number of replicates in the screen\n", "decides which one applies:\n", "\n", "- `\"signature_stability\"` correlates two disjoint subsets. It is unbiased and needs\n", " `2 * depth` replicates.\n", "- `\"convergence\"` correlates a subset with the group's full signature. It is defined up to\n", " one less than the group size and is optimistic by construction.\n", "\n", "BBBC021 has at most three replicates per treatment, so the unbiased measure gives a single point." ] }, { "cell_type": "code", "execution_count": 8, "id": "c016", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:14:57.738299Z", "iopub.status.busy": "2026-10-21T01:14:57.738115Z", "iopub.status.idle": "2026-10-21T01:14:59.451792Z", "shell.execute_reply": "2026-10-21T01:14:59.451148Z" } }, "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_replicatesmeanstdn_draws
010.7910.00910
120.8730.00510
\n", "
" ], "text/plain": [ " n_replicates mean std n_draws\n", "0 1 0.791 0.009 10\n", "1 2 0.873 0.005 10" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mt.tl.replicate_saturation(treated, metric=\"convergence\", n_draws=10)\n", "treated.uns[\"mantispy\"][\"replicate_saturation\"].round(3)" ] }, { "cell_type": "code", "execution_count": 9, "id": "c017", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:14:59.453692Z", "iopub.status.busy": "2026-10-21T01:14:59.453543Z", "iopub.status.idle": "2026-10-21T01:15:00.243655Z", "shell.execute_reply": "2026-10-21T01:15:00.242493Z" } }, "outputs": [ { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mt.pl.replicate_saturation(treated);" ] }, { "cell_type": "markdown", "id": "c018", "metadata": {}, "source": [ "From one replicate to two, signature agreement rises from 0.79 to 0.87 and is still rising at\n", "the largest depth the data supports. For an experiment that has to distinguish mechanisms,\n", "three replicates is on the low side, which is useful when planning the next screen." ] }, { "cell_type": "markdown", "id": "c019", "metadata": {}, "source": [ "## Summary\n", "\n", "- mAP against the controls says whether a perturbation did anything; compare it with its null, not with an\n", " intuition about scores between zero and one.\n", "- Reproducible and generalizable are different properties, and a filter that improves one can hurt the\n", " other. Measure the task you care about.\n", "- Plates differ even where they received the same treatment. Test for it on the controls.\n", "- A saturation curve that is still rising says the next experiment needs more replicates." ] } ], "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 }