{ "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 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." ] }, { "cell_type": "code", "execution_count": 1, "id": "c001", "metadata": { "execution": { "iopub.execute_input": "2026-10-27T22:17:01.858606Z", "iopub.status.busy": "2026-10-27T22:17:01.858205Z", "iopub.status.idle": "2026-10-27T22:18:17.047914Z", "shell.execute_reply": "2026-10-27T22:18:17.045970Z" } }, "outputs": [], "source": [ "import mantispy as mt" ] }, { "cell_type": "markdown", "id": "c002", "metadata": {}, "source": [ "## A real screen\n", "\n", "`pki` is the JUMP pilot's kinase inhibitor set {cite:p}`Chandrasekaran_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](../compounds/dose_response.ipynb)." ] }, { "cell_type": "code", "execution_count": 2, "id": "c003", "metadata": { "execution": { "iopub.execute_input": "2026-10-27T22:18:17.052375Z", "iopub.status.busy": "2026-10-27T22:18:17.051882Z", "iopub.status.idle": "2026-10-27T22:18:38.402388Z", "shell.execute_reply": "2026-10-27T22:18:38.401122Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":2: UserWarning: 133 of 5839 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": [ "(2112, 850)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata = mt.ds.pki()\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": "892fe9287a2d", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "id": "f394ae3e3e14", "metadata": {}, "source": [ "## Did the cells survive?\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 3, "id": "07f51a189343", "metadata": { "execution": { "iopub.execute_input": "2026-10-27T22:18:38.405848Z", "iopub.status.busy": "2026-10-27T22:18:38.405395Z", "iopub.status.idle": "2026-10-27T22:18:38.467272Z", "shell.execute_reply": "2026-10-27T22:18:38.465903Z" } }, "outputs": [ { "data": { "text/plain": [ "{'control median cells/well': 1679,\n", " 'treated median cells/well': 1664,\n", " 'treated wells below half the control median': 578,\n", " 'of treated wells': 2112}" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "\n", "counts = adata.obs[\"Metadata_CellCount\"].to_numpy(float)\n", "is_control = adata.obs[\"Metadata_Control\"].to_numpy(bool)\n", "control_median = np.median(counts[is_control])\n", "{\n", " \"control median cells/well\": int(control_median),\n", " \"treated median cells/well\": int(np.median(counts[~is_control])),\n", " \"treated wells below half the control median\": int((counts[~is_control] < 0.5 * control_median).sum()),\n", " \"of treated wells\": int((~is_control).sum()),\n", "}" ] }, { "cell_type": "code", "execution_count": 4, "id": "424d5b9cdd73", "metadata": { "execution": { "iopub.execute_input": "2026-10-27T22:18:38.470664Z", "iopub.status.busy": "2026-10-27T22:18:38.470060Z", "iopub.status.idle": "2026-10-27T22:19:00.508464Z", "shell.execute_reply": "2026-10-27T22:19:00.507240Z" } }, "outputs": [ { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mt.pl.cell_counts(adata, groupby=\"Metadata_Control\", labels={True: \"control\", False: \"treated\"});" ] }, { "cell_type": "markdown", "id": "af5ea91e4d06", "metadata": {}, "source": [ "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](../compounds/hits.ipynb) follows the flagged compounds down the dose ladder with `tl.cytotoxicity`, and [Plate artifacts](../profiles/artifacts.ipynb) regresses cell count out of the features. The rule: if the strongest hits are also the emptiest wells, suspect cytotoxicity, not phenotype." ] }, { "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 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.\n", "\n", "`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](https://github.com/cytomining/copairs) {cite:p}`Kalinin_2025`, the Broad's implementation, 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 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." ] }, { "cell_type": "code", "execution_count": 5, "id": "c005", "metadata": { "execution": { "iopub.execute_input": "2026-10-27T22:19:00.514955Z", "iopub.status.busy": "2026-10-27T22:19:00.514005Z", "iopub.status.idle": "2026-10-27T22:19:17.485353Z", "shell.execute_reply": "2026-10-27T22:19:17.483892Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "36 of 37 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", " \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", " \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
31BRD-U00086681-001-01-9@0.040.4559780.071393
25BRD-U00086679-001-01-9@0.20.4726220.035866
21BRD-U00086677-001-01-9@0.040.4809760.027800
19BRD-U00086677-001-01-9@0.0040.4828190.027312
7BRD-U00086673-001-01-9@0.20.4913240.017713
26BRD-U00086679-001-01-9@0.40.4936170.015955
22BRD-U00086678-001-01-9@0.0040.4936890.015955
20BRD-U00086677-001-01-9@0.010.4992540.013072
32BRD-U00086681-001-01-9@0.20.5000120.012885
28BRD-U00086680-001-01-9@0.040.5049080.009778
5BRD-U00086672-001-01-9@0.40.5090000.007947
33BRD-U00086681-001-01-9@0.40.5153770.006119
4BRD-U00086672-001-01-9@0.20.5233920.004292
23BRD-U00086678-001-01-9@0.010.5255360.004292
8BRD-U00086673-001-01-9@0.40.5258810.004292
2BRD-K40109029-001-06-9@2.00.5463320.000185
24BRD-U00086678-001-01-9@0.040.5539900.001009
27BRD-U00086679-001-01-9@1.00.5941870.000352
6BRD-U00086672-001-01-9@1.00.6048750.000185
9BRD-U00086673-001-01-9@1.00.6106270.000185
10BRD-U00086674-001-01-9@0.040.6425090.000185
34BRD-U00086682-001-01-9@0.040.6436230.000185
29BRD-U00086680-001-01-9@0.20.6517590.000185
13BRD-U00086675-001-01-9@0.040.6940400.000185
16BRD-U00086676-001-01-9@0.0040.6999550.000185
30BRD-U00086680-001-01-9@0.40.7292620.000185
17BRD-U00086676-001-01-9@0.010.7529510.000185
18BRD-U00086676-001-01-9@0.040.7557860.000185
1BRD-K15819326-001-01-6@2.00.7558160.000185
14BRD-U00086675-001-01-9@0.20.7974460.000185
11BRD-U00086674-001-01-9@0.20.8030300.000185
15BRD-U00086675-001-01-9@0.40.8257110.000185
12BRD-U00086674-001-01-9@0.40.8598700.000185
35BRD-U00086682-001-01-9@0.20.8687330.000185
0BRD-K15179513-001-03-4@2.00.9152010.000185
3BRD-K95785537-001-26-9@2.00.9430410.000185
36BRD-U00086682-001-01-9@0.40.9725190.000185
\n", "
" ], "text/plain": [ " Metadata_Perturbation score qvalue\n", "31 BRD-U00086681-001-01-9@0.04 0.455978 0.071393\n", "25 BRD-U00086679-001-01-9@0.2 0.472622 0.035866\n", "21 BRD-U00086677-001-01-9@0.04 0.480976 0.027800\n", "19 BRD-U00086677-001-01-9@0.004 0.482819 0.027312\n", "7 BRD-U00086673-001-01-9@0.2 0.491324 0.017713\n", "26 BRD-U00086679-001-01-9@0.4 0.493617 0.015955\n", "22 BRD-U00086678-001-01-9@0.004 0.493689 0.015955\n", "20 BRD-U00086677-001-01-9@0.01 0.499254 0.013072\n", "32 BRD-U00086681-001-01-9@0.2 0.500012 0.012885\n", "28 BRD-U00086680-001-01-9@0.04 0.504908 0.009778\n", "5 BRD-U00086672-001-01-9@0.4 0.509000 0.007947\n", "33 BRD-U00086681-001-01-9@0.4 0.515377 0.006119\n", "4 BRD-U00086672-001-01-9@0.2 0.523392 0.004292\n", "23 BRD-U00086678-001-01-9@0.01 0.525536 0.004292\n", "8 BRD-U00086673-001-01-9@0.4 0.525881 0.004292\n", "2 BRD-K40109029-001-06-9@2.0 0.546332 0.000185\n", "24 BRD-U00086678-001-01-9@0.04 0.553990 0.001009\n", "27 BRD-U00086679-001-01-9@1.0 0.594187 0.000352\n", "6 BRD-U00086672-001-01-9@1.0 0.604875 0.000185\n", "9 BRD-U00086673-001-01-9@1.0 0.610627 0.000185\n", "10 BRD-U00086674-001-01-9@0.04 0.642509 0.000185\n", "34 BRD-U00086682-001-01-9@0.04 0.643623 0.000185\n", "29 BRD-U00086680-001-01-9@0.2 0.651759 0.000185\n", "13 BRD-U00086675-001-01-9@0.04 0.694040 0.000185\n", "16 BRD-U00086676-001-01-9@0.004 0.699955 0.000185\n", "30 BRD-U00086680-001-01-9@0.4 0.729262 0.000185\n", "17 BRD-U00086676-001-01-9@0.01 0.752951 0.000185\n", "18 BRD-U00086676-001-01-9@0.04 0.755786 0.000185\n", "1 BRD-K15819326-001-01-6@2.0 0.755816 0.000185\n", "14 BRD-U00086675-001-01-9@0.2 0.797446 0.000185\n", "11 BRD-U00086674-001-01-9@0.2 0.803030 0.000185\n", "15 BRD-U00086675-001-01-9@0.4 0.825711 0.000185\n", "12 BRD-U00086674-001-01-9@0.4 0.859870 0.000185\n", "35 BRD-U00086682-001-01-9@0.2 0.868733 0.000185\n", "0 BRD-K15179513-001-03-4@2.0 0.915201 0.000185\n", "3 BRD-K95785537-001-26-9@2.0 0.943041 0.000185\n", "36 BRD-U00086682-001-01-9@0.4 0.972519 0.000185" ] }, "execution_count": 5, "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": "code", "execution_count": 6, "id": "09015b24fa6d", "metadata": { "execution": { "iopub.execute_input": "2026-10-27T22:19:17.488867Z", "iopub.status.busy": "2026-10-27T22:19:17.488219Z", "iopub.status.idle": "2026-10-27T22:19:17.928590Z", "shell.execute_reply": "2026-10-27T22:19:17.927484Z" } }, "outputs": [ { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mt.pl.map(adata);" ] }, { "cell_type": "markdown", "id": "c006", "metadata": {}, "source": [ "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 `@`, 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.\n", "\n", "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." ] }, { "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": 7, "id": "c008", "metadata": { "execution": { "iopub.execute_input": "2026-10-27T22:19:17.932037Z", "iopub.status.busy": "2026-10-27T22:19:17.931721Z", "iopub.status.idle": "2026-10-27T22:19:18.128003Z", "shell.execute_reply": "2026-10-27T22:19:18.127032Z" } }, "outputs": [ { "data": { "text/plain": [ "{'median': 0.233,\n", " 'lower quartile': 0.136,\n", " 'upper quartile': 0.34,\n", " 'features above 0.2': 495,\n", " 'of': 850}" ] }, "execution_count": 7, "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": 8, "id": "c010", "metadata": { "execution": { "iopub.execute_input": "2026-10-27T22:19:18.131191Z", "iopub.status.busy": "2026-10-27T22:19:18.131013Z", "iopub.status.idle": "2026-10-27T22:19:20.553444Z", "shell.execute_reply": "2026-10-27T22:19:20.552194Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":4: UserWarning: with null_size=500 no p-value can fall below 1/501, so the correction over 33 groups calls none of them unless at least 2 reach that floor together. Raise null_size above 660 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: with null_size=500 no p-value can fall below 1/501, so the correction over 33 groups calls none of them unless at least 2 reach that floor together. Raise null_size above 660 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: with null_size=500 no p-value can fall below 1/501, so the correction over 33 groups calls none of them unless at least 2 reach that floor together. Raise null_size above 660 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.238, 'ICC > 0.2': 0.247, 'ICC > 0.4': 0.26}" ] }, "execution_count": 8, "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 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.\n", "\n", "Mechanism retrieval, whether a compound finds a different compound with the same mechanism, is the task from [Mechanism of action](../compounds/mechanism_of_action.ipynb). On this screen it cannot be measured at all." ] }, { "cell_type": "markdown", "id": "eb2a73641c7a", "metadata": {}, "source": [ "The next cell returns 0.00 on `pki` by design, not because anything failed; the section after it explains why." ] }, { "cell_type": "code", "execution_count": 9, "id": "c012", "metadata": { "execution": { "iopub.execute_input": "2026-10-27T22:19:20.556361Z", "iopub.status.busy": "2026-10-27T22:19:20.556165Z", "iopub.status.idle": "2026-10-27T22:19:21.227029Z", "shell.execute_reply": "2026-10-27T22:19:21.225407Z" } }, "outputs": [ { "data": { "text/plain": [ "{'all features': 0.0, 'ICC > 0.2': 0.0, 'ICC > 0.4': 0.0}" ] }, "execution_count": 9, "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": [ "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](../compounds/mechanism_of_action.ipynb), is where ICC can be watched to help or hurt that task.\n", "\n", "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).\n", "\n", "## How much of this is batch effect?\n", "\n", "`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." ] }, { "cell_type": "code", "execution_count": 10, "id": "c014", "metadata": { "execution": { "iopub.execute_input": "2026-10-27T22:19:21.230626Z", "iopub.status.busy": "2026-10-27T22:19:21.230146Z", "iopub.status.idle": "2026-10-27T22:19:23.674219Z", "shell.execute_reply": "2026-10-27T22:19:23.673353Z" } }, "outputs": [ { "data": { "text/plain": [ "{'controls, before normalization': 850,\n", " 'treated, after normalization': 847,\n", " 'of': 850}" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "raw = mt.ds.pki()[:, 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", "treated_normalized = adata[~adata.obs[\"Metadata_Control\"].to_numpy()].copy()\n", "mt.pp.feature_batch_sensitivity(treated_normalized, batch_key=\"Metadata_Plate\")\n", "{\n", " \"controls, before normalization\": int(controls.var[\"batch_sensitive\"].sum()),\n", " \"treated, after normalization\": int(treated_normalized.var[\"batch_sensitive\"].sum()),\n", " \"of\": int(controls.n_vars),\n", "}" ] }, { "cell_type": "markdown", "id": "c015", "metadata": {}, "source": [ "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](../profiles/artifacts.ipynb) does this and corrects what remains.\n", "\n", "## How many replicates do I need?\n", "\n", "`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:\n", "\n", "- `\"signature_stability\"` (the default) correlates two disjoint subsets. It is unbiased and needs `2 * depth` replicates.\n", "- `\"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.\n", "\n", "`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." ] }, { "cell_type": "code", "execution_count": 11, "id": "c016", "metadata": { "execution": { "iopub.execute_input": "2026-10-27T22:19:23.677480Z", "iopub.status.busy": "2026-10-27T22:19:23.677242Z", "iopub.status.idle": "2026-10-27T22:19:36.560426Z", "shell.execute_reply": "2026-10-27T22:19:36.559583Z" } }, "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", " \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", " \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", " \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_replicatesmeanstdn_draws
010.1730.03910
120.3160.03710
230.3440.03810
340.4500.04810
450.4510.02710
560.5250.02110
670.5110.01910
780.5970.02610
890.5860.02710
9100.6360.02310
10110.6330.02010
11120.6660.03110
12130.6670.01610
13140.6970.02110
14150.7000.01910
15160.7120.02410
16170.7150.02310
17180.7320.01610
18190.7400.02410
19200.7530.02210
20210.7460.01310
21220.7670.01910
22230.7750.02310
23240.7810.01910
24250.7890.00910
25260.7990.01910
26270.8080.01210
27280.8110.01410
28290.8070.01110
29300.8150.01410
30310.8250.01610
31320.8190.01210
\n", "
" ], "text/plain": [ " n_replicates mean std n_draws\n", "0 1 0.173 0.039 10\n", "1 2 0.316 0.037 10\n", "2 3 0.344 0.038 10\n", "3 4 0.450 0.048 10\n", "4 5 0.451 0.027 10\n", "5 6 0.525 0.021 10\n", "6 7 0.511 0.019 10\n", "7 8 0.597 0.026 10\n", "8 9 0.586 0.027 10\n", "9 10 0.636 0.023 10\n", "10 11 0.633 0.020 10\n", "11 12 0.666 0.031 10\n", "12 13 0.667 0.016 10\n", "13 14 0.697 0.021 10\n", "14 15 0.700 0.019 10\n", "15 16 0.712 0.024 10\n", "16 17 0.715 0.023 10\n", "17 18 0.732 0.016 10\n", "18 19 0.740 0.024 10\n", "19 20 0.753 0.022 10\n", "20 21 0.746 0.013 10\n", "21 22 0.767 0.019 10\n", "22 23 0.775 0.023 10\n", "23 24 0.781 0.019 10\n", "24 25 0.789 0.009 10\n", "25 26 0.799 0.019 10\n", "26 27 0.808 0.012 10\n", "27 28 0.811 0.014 10\n", "28 29 0.807 0.011 10\n", "29 30 0.815 0.014 10\n", "30 31 0.825 0.016 10\n", "31 32 0.819 0.012 10" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mt.tl.replicate_saturation(treated, n_draws=10)\n", "treated.uns[\"mantispy\"][\"replicate_saturation\"].round(3)" ] }, { "cell_type": "code", "execution_count": 12, "id": "c017", "metadata": { "execution": { "iopub.execute_input": "2026-10-27T22:19:36.563238Z", "iopub.status.busy": "2026-10-27T22:19:36.563058Z", "iopub.status.idle": "2026-10-27T22:19:36.693095Z", "shell.execute_reply": "2026-10-27T22:19:36.692108Z" } }, "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, 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." ] }, { "cell_type": "markdown", "id": "c019", "metadata": {}, "source": [ "## Summary\n", "\n", "A quick read on whether a screen is worth analysing:\n", "\n", "- **Cells survived?** If the strongest hits are also the emptiest wells, suspect cytotoxicity, not phenotype.\n", "- **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.\n", "- **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.\n", "- **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.\n", "- **Enough replicates?** A saturation curve still rising at your maximum depth says the next screen needs more wells per treatment." ] } ], "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 }