{ "cells": [ { "cell_type": "markdown", "id": "36f1c57f", "metadata": {}, "source": [ "# Hit callers on cells and on guide profiles\n", "\n", "{func}`~mantispy.tl.hit_calling` offers four statistics, and a screen can be tested at two stages: on the single\n", "cells, or on one median profile per guide. This page runs every method that applies at each stage on the same screen,\n", "and compares the calls with a set of genes whose answer is known.\n", "\n", "The screen is {func}`~mantispy.ds.scallops_arv471`, an optical pooled CRISPR screen in MCF7 cells treated with the\n", "estrogen-receptor degrader ARV-471. The drug needs the CRL4-CRBN ligase to degrade the receptor, so knocking out\n", "`CRBN`, `DDB1`, `CUL4A`, `CUL4B` or the receptor itself (`ESR1`) changes the phenotype: these five are the known\n", "positives. Twelve olfactory-receptor genes are not expressed in these cells and are the known negatives. The other\n", "genes of the library have no known answer, so they are counted but not scored.\n", "[Hit callers on well profiles](../compounds/hit_calling_methods.ipynb) compares the default with `\"mp_value\"` on\n", "compound screens." ] }, { "cell_type": "code", "execution_count": 1, "id": "058bb15b", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T03:54:46.883767Z", "iopub.status.busy": "2026-11-11T03:54:46.883370Z", "iopub.status.idle": "2026-11-11T03:56:30.965880Z", "shell.execute_reply": "2026-11-11T03:56:30.964399Z" } }, "outputs": [ { "data": { "text/plain": [ "AnnData object with n_obs × n_vars = 135549 × 9\n", " obs: 'Metadata_Gene', 'Metadata_sgRNA', 'Metadata_Control_Type', 'Metadata_Control', 'Metadata_Perturbation', 'Metadata_Plate', 'Metadata_Well', 'Metadata_Row', 'Metadata_Col', 'Metadata_ImageID', 'Metadata_ObjectType', 'Metadata_ObjectNumber', 'Metadata_Perturbation_Type'\n", " var: 'object', 'feature_group', 'feature', 'channel', 'scale', 'angle', 'gray_levels', 'radial_bin', 'params', 'is_feature', 'feature_kind', 'degenerate_scale'\n", " uns: 'mantispy'\n", " layers: None (.X)" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import time\n", "\n", "import pandas as pd\n", "\n", "import mantispy as mt\n", "\n", "POSITIVES = [\"CRBN\", \"DDB1\", \"CUL4A\", \"CUL4B\", \"ESR1\"]\n", "\n", "cells = mt.ds.scallops_arv471()\n", "mt.pp.normalize(cells, by=\"Metadata_Plate\", reference=\"negcon\")\n", "genes = cells.obs[\"Metadata_Gene\"].astype(str)\n", "NEGATIVES = sorted(gene for gene in genes.unique() if gene.startswith(\"OR\"))\n", "# Two million cells is more than the comparison needs; keep 200 per guide.\n", "cells = mt.pp.downsample(cells, n_per_group=200, by=\"Metadata_sgRNA\")\n", "cells" ] }, { "cell_type": "code", "execution_count": 2, "id": "00eb44a7", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T03:56:30.969643Z", "iopub.status.busy": "2026-11-11T03:56:30.969440Z", "iopub.status.idle": "2026-11-11T03:56:30.975957Z", "shell.execute_reply": "2026-11-11T03:56:30.974887Z" } }, "outputs": [], "source": [ "def score(calls):\n", " \"\"\"Confusion counts of each method's calls against the known positives and negatives.\"\"\"\n", " rows = {}\n", " for method, called in calls.items():\n", " rows[method] = {\n", " \"positives called\": int(called.reindex(POSITIVES).sum()),\n", " \"positives missed\": int((~called.reindex(POSITIVES)).sum()),\n", " \"negatives called\": int(called.reindex(NEGATIVES).sum()),\n", " \"negatives not called\": int((~called.reindex(NEGATIVES)).sum()),\n", " \"other genes called\": int(called.drop(POSITIVES + NEGATIVES + [\"nontargeting\"], errors=\"ignore\").sum()),\n", " }\n", " return pd.DataFrame(rows).T\n", "\n", "\n", "def agreement(calls):\n", " \"\"\"Share of genes on which two methods make the same call.\"\"\"\n", " frame = pd.DataFrame(calls)\n", " return pd.DataFrame({a: {b: (frame[a] == frame[b]).mean() for b in frame} for a in frame}).round(2)\n", "\n", "\n", "def call(adata, methods):\n", " \"\"\"Run each method on the genes and return its calls at q < 0.05 and its runtime.\"\"\"\n", " calls, seconds = {}, {}\n", " for name, keywords in methods.items():\n", " start = time.perf_counter()\n", " mt.tl.hit_calling(adata, groupby=\"Metadata_Gene\", n_permutations=1000, key_added=name, **keywords)\n", " seconds[name] = round(time.perf_counter() - start, 1)\n", " calls[name] = adata.uns[\"mantispy\"][name].set_index(\"group\")[\"qvalue\"] < 0.05\n", " return calls, pd.Series(seconds, name=\"seconds\")" ] }, { "cell_type": "markdown", "id": "97b94de1", "metadata": {}, "source": [ "## Single cells\n", "\n", "Every well of a pooled screen holds every guide, so the unit the null draws whole is the guide (`replicate=`).\n", "The four statistics ask different questions of a gene's cells:\n", "\n", "- `\"mahalanobis\"`: how far the gene's typical cell sits from the non-targeting cells, in the controls' covariance.\n", "- `\"ks\"`: whether the distribution of the cells' distances differs, which also catches a subpopulation that moved.\n", "- `\"hellinger\"`: how little a Gaussian fitted to the gene's cells overlaps one fitted to the controls, so a change in\n", " spread or shape counts as well as a shift. It follows the robust morphological perturbation value of BioProfiling.jl\n", " (Vulliard et al., *Bioinformatics* 2022, [doi:10.1093/bioinformatics/btab853](https://doi.org/10.1093/bioinformatics/btab853)),\n", " which builds on the multidimensional perturbation value (Hutz et al., *J Biomol Screen* 2013,\n", " [doi:10.1177/1087057112469257](https://doi.org/10.1177/1087057112469257)).\n", " `covariance=\"robust\"` gives BioProfiling.jl's robust variant; it refits the fits on every permutation and takes far\n", " longer, so it is left out here.\n", "\n", " The statistic is one minus the overlap of the two Gaussians. That overlap is a term for how much the two\n", " spreads differ, times exp(-Δμᵀ S⁻¹ Δμ / 8), where Δμ is the difference between the means and S is the average of\n", " the two covariances. The second factor is the squared Mahalanobis distance in the pooled covariance, the form the\n", " multidimensional perturbation value uses, so a pure shift in the mean still counts, and a change in spread adds to it.\n", "- `\"mp_value\"`: that squared Mahalanobis distance alone, between the gene's centroid and the controls', in the\n", " average of the two covariances. This is the multidimensional perturbation value itself. Only the average has to be\n", " positive definite, so a gene needs two rows rather than twice as many rows as features." ] }, { "cell_type": "code", "execution_count": 3, "id": "9050d9bf", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T03:56:30.978814Z", "iopub.status.busy": "2026-11-11T03:56:30.978656Z", "iopub.status.idle": "2026-11-11T03:58:12.552802Z", "shell.execute_reply": "2026-11-11T03:58:12.551037Z" } }, "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", "
positives calledpositives missednegatives callednegatives not calledother genes calledseconds
mahalanobis23111467.1
ks231115523.9
hellinger23012566.8
mp_value41931103.2
\n", "
" ], "text/plain": [ " positives called positives missed negatives called \\\n", "mahalanobis 2 3 1 \n", "ks 2 3 1 \n", "hellinger 2 3 0 \n", "mp_value 4 1 9 \n", "\n", " negatives not called other genes called seconds \n", "mahalanobis 11 46 7.1 \n", "ks 11 55 23.9 \n", "hellinger 12 5 66.8 \n", "mp_value 3 110 3.2 " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cell_calls, cell_seconds = call(\n", " cells,\n", " {\n", " \"mahalanobis\": {\"replicate\": \"Metadata_sgRNA\"},\n", " \"ks\": {\"method\": \"ks\", \"replicate\": \"Metadata_sgRNA\"},\n", " \"hellinger\": {\"method\": \"hellinger\", \"replicate\": \"Metadata_sgRNA\"},\n", " \"mp_value\": {\"method\": \"mp_value\", \"replicate\": \"Metadata_sgRNA\"},\n", " },\n", ")\n", "score(cell_calls).assign(seconds=cell_seconds)" ] }, { "cell_type": "code", "execution_count": 4, "id": "b043e460", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T03:58:12.556170Z", "iopub.status.busy": "2026-11-11T03:58:12.555718Z", "iopub.status.idle": "2026-11-11T03:58:12.572836Z", "shell.execute_reply": "2026-11-11T03:58:12.571524Z" } }, "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", "
mahalanobiskshellingermp_value
mahalanobis1.000.910.740.53
ks0.911.000.680.58
hellinger0.740.681.000.27
mp_value0.530.580.271.00
\n", "
" ], "text/plain": [ " mahalanobis ks hellinger mp_value\n", "mahalanobis 1.00 0.91 0.74 0.53\n", "ks 0.91 1.00 0.68 0.58\n", "hellinger 0.74 0.68 1.00 0.27\n", "mp_value 0.53 0.58 0.27 1.00" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "agreement(cell_calls)" ] }, { "cell_type": "markdown", "id": "32fcb518", "metadata": {}, "source": [ "`\"mahalanobis\"`, `\"ks\"` and `\"hellinger\"` each call two of the five positives and miss the other three.\n", "`\"mahalanobis\"` and `\"ks\"` each call one of the twelve negatives, and `\"hellinger\"` calls none. `\"mp_value\"` calls four\n", "of the five positives, but also nine of the twelve negatives. The negatives are guides that cut genes the cells do not\n", "express, so a method that calls them is reading a difference from the non-targeting guides that has nothing to do with\n", "the target. `\"mp_value\"` tests the centre of all of a gene's cells directly, which makes it the most likely of the\n", "four to pick up a small shift that every cutting guide shares (see the end of this page).\n", "\n", "The methods differ most among the genes without a known answer: `\"hellinger\"` calls 5 of them, `\"mahalanobis\"` 46,\n", "`\"ks\"` 55 and `\"mp_value\"` 110. The truth set cannot say which of those calls are right, but it shows that\n", "`\"hellinger\"` returns the shortest list at the same recall on the known positives. It is also the slowest of the four,\n", "and `\"mp_value\"` the fastest. `\"mahalanobis\"` and `\"ks\"` agree on 91% of the genes, and `\"hellinger\"` agrees with them\n", "on 74% and 68%. `\"mp_value\"` agrees with them on 53% and 58%, and with `\"hellinger\"` on 27%." ] }, { "cell_type": "markdown", "id": "d754bdde", "metadata": {}, "source": [ "## Guide profiles\n", "\n", "The same normalized cells, collapsed to one median profile per guide with {func}`~mantispy.tl.aggregate`, so that the\n", "two stages see the same measurements. A gene now has a handful of rows, its guides. `\"hellinger\"` needs a covariance\n", "within each group, and a covariance in nine features needs at least twice as many rows as features on each side, so it\n", "does not apply at this stage. `\"mp_value\"` needs only the average of the gene's and the controls' covariances to be\n", "positive definite, so it runs here alongside the other two, and each guide profile is its own unit." ] }, { "cell_type": "code", "execution_count": 5, "id": "c89eecba", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T03:58:12.575185Z", "iopub.status.busy": "2026-11-11T03:58:12.575014Z", "iopub.status.idle": "2026-11-11T03:58:16.541746Z", "shell.execute_reply": "2026-11-11T03:58:16.539961Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":26: UserWarning: 157 group(s) were left unscored because the group had fewer than 2 complete rows, the reference no more than features, or the averaged covariance was singular: AMBRA1, AMD1, ARID1A, BARD1, BRCA1, .... Score a reduced representation with use_rep, e.g. 'X_pca'.\n", " mt.tl.hit_calling(adata, groupby=\"Metadata_Gene\", n_permutations=1000, key_added=name, **keywords)\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", "
positives calledpositives missednegatives callednegatives not calledother genes calledseconds
mahalanobis1401230.1
ks23012110.1
mp_value1401211.6
\n", "
" ], "text/plain": [ " positives called positives missed negatives called \\\n", "mahalanobis 1 4 0 \n", "ks 2 3 0 \n", "mp_value 1 4 0 \n", "\n", " negatives not called other genes called seconds \n", "mahalanobis 12 3 0.1 \n", "ks 12 11 0.1 \n", "mp_value 12 1 1.6 " ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "guides = mt.tl.aggregate(cells, by=(\"Metadata_Gene\", \"Metadata_sgRNA\"), func=\"median\")\n", "guide_calls, guide_seconds = call(\n", " guides, {\"mahalanobis\": {}, \"ks\": {\"method\": \"ks\"}, \"mp_value\": {\"method\": \"mp_value\"}}\n", ")\n", "score(guide_calls).assign(seconds=guide_seconds)" ] }, { "cell_type": "code", "execution_count": 6, "id": "aa473707", "metadata": { "execution": { "iopub.execute_input": "2026-11-11T03:58:16.544897Z", "iopub.status.busy": "2026-11-11T03:58:16.544488Z", "iopub.status.idle": "2026-11-11T03:58:16.559464Z", "shell.execute_reply": "2026-11-11T03:58:16.558188Z" } }, "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", "
mahalanobisksmp_value
mahalanobis1.000.930.96
ks0.931.000.92
mp_value0.960.921.00
\n", "
" ], "text/plain": [ " mahalanobis ks mp_value\n", "mahalanobis 1.00 0.93 0.96\n", "ks 0.93 1.00 0.92\n", "mp_value 0.96 0.92 1.00" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "agreement(guide_calls)" ] }, { "cell_type": "markdown", "id": "102c65bc", "metadata": {}, "source": [ "On guide profiles no method calls a negative. `\"mahalanobis\"` finds one positive and `\"ks\"` two, along with 3 and 11\n", "other genes, and each runs in a fraction of a second because a gene now has four rows instead of eight hundred.\n", "`\"mp_value\"` leaves 156 of the 158 genes unscored, with a warning. A spot count's median is a whole number, and one of\n", "them, `Nuclei_Spots_Count_ESR1`, takes the same value for every non-targeting guide, so the controls do not vary in it\n", "and the averaged covariance is singular. It calls both genes it could score, one of them a positive.\n", "\n", "## Cells or profiles\n", "\n", "On this screen `\"mahalanobis\"` recovers more positives from the cells than from the guide medians, and `\"ks\"` recovers\n", "two at both stages. Guide profiles are quick to test and need nothing beyond {func}`~mantispy.tl.aggregate`.\n", "Cells take longer, but they are the only stage where `\"hellinger\"` applies and where a change in the spread of a gene's\n", "cells can show at all. `\"mp_value\"` finds the most positives on cells and calls most of the negatives with them, and on\n", "guide medians of these counts it cannot score most genes. `\"mahalanobis\"` stays the default because here it finds as\n", "many positives as `\"ks\"` or `\"hellinger\"` at a fraction of their runtime, without the negatives `\"mp_value\"` calls.\n", "`\"hellinger\"` is the choice when a shorter list, built on the shape of the cell cloud as well as its centre, matters\n", "more than speed.\n", "\n", "Every targeting guide cuts DNA and a non-targeting guide does not, so against non-targeting controls the cutting\n", "response alone can make most knockouts look active, here and in `tl.map(mode=\"activity\")`. Guides against genes the\n", "cells do not express, such as the olfactory receptors above, cut without changing a target, so they separate a gene's\n", "effect from cutting. {func}`~mantispy.io.unexpressed_genes` lists those genes for a cell line, and\n", "{func}`~mantispy.tl.empirical_fdr` calibrates the calls against them." ] } ], "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 }