{ "cells": [ { "cell_type": "markdown", "id": "cell-0", "metadata": {}, "source": [ "# Calling gene hits in a pooled screen\n", "\n", "A pooled optical CRISPR screen targets each gene with many single guides, read out one cell at a time. A per-guide phenotype is noisy and a few guides of any gene are duds, so a gene call has to pool a gene's guides rather than trust any one of them. Each gene call asks whether its guides, taken together, move further than non-targeting guides do.\n", "\n", "This case study scores the guides of the `cp_posh` screen for phenotypic activity with {func}`~mantispy.tl.guide_activity`, then calls genes with {func}`~mantispy.tl.aggregate_guides`, which combines a gene's guides and calibrates the combined statistic against random same-size groups of non-targeting-control guides. The screen carries two control classes, non-targeting and intergenic, so one sets the per-guide activity scale while the other is the null the gene call is judged against. The two have to be different classes, or the null is scored against itself." ] }, { "cell_type": "code", "execution_count": 1, "id": "cell-1", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T00:50:18.549946Z", "iopub.status.busy": "2026-10-21T00:50:18.549558Z", "iopub.status.idle": "2026-10-21T00:50:22.437054Z", "shell.execute_reply": "2026-10-21T00:50:22.435833Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AnnData object with n_obs × n_vars = 1551 × 1278\n", " obs: 'Metadata_Gene', 'Metadata_sgRNA', 'Metadata_CellCount', 'Metadata_Perturbation', 'Metadata_Plate', 'Metadata_Control', 'Metadata_Control_Type', 'Metadata_ObjectType', 'Metadata_Perturbation_Type'\n", " var: 'object', 'feature_group', 'feature', 'channel', 'scale', 'angle', 'gray_levels', 'radial_bin', 'params', 'is_feature', 'feature_kind'\n", " uns: 'mantispy'\n", " layers: None (.X)\n" ] }, { "data": { "text/plain": [ "Metadata_Gene\n", "nontargeting 200\n", "intergenic 199\n", "ACIN1 10\n", "ACO2 10\n", "Name: count, dtype: int64" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import mantispy as mt\n", "\n", "guides = mt.ds.cp_posh(aggregated=True)\n", "print(guides)\n", "guides.obs[\"Metadata_Gene\"].value_counts().head(4)" ] }, { "cell_type": "markdown", "id": "cell-2", "metadata": {}, "source": [ "One profile per guide, about ten guides per gene, and two control classes: `nontargeting` guides that carry a guide hitting nothing, and `intergenic` guides that cut a gene-free locus. Both should be phenotypically silent, which is what lets one play the activity reference and the other the null." ] }, { "cell_type": "code", "execution_count": 2, "id": "cell-3", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T00:50:22.440572Z", "iopub.status.busy": "2026-10-21T00:50:22.440377Z", "iopub.status.idle": "2026-10-21T00:50:22.515903Z", "shell.execute_reply": "2026-10-21T00:50:22.514761Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "scored 1352 guides against the intergenic reference\n" ] }, { "data": { "text/plain": [ "min 0.035\n", "50% 0.560\n", "max 0.995\n", "Name: guide_activity, dtype: float64" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Per-guide activity: how far each guide sits from the intergenic control centre, as a one-sided\n", "# p-value small for an active guide. The intergenic class sets the scale and is left unscored.\n", "mt.tl.guide_activity(guides, reference=\"intergenic\")\n", "scored = guides.obs[\"guide_activity\"].notna().sum()\n", "print(f\"scored {scored} guides against the intergenic reference\")\n", "guides.obs.loc[guides.obs[\"Metadata_Gene\"] == \"nontargeting\", \"guide_activity\"].describe()[[\"min\", \"50%\", \"max\"]]" ] }, { "cell_type": "markdown", "id": "cell-4", "metadata": {}, "source": [ "The non-targeting guides, scored against the intergenic reference, sit near the middle of the range: a control class carries no phenotype, so its activity p-values spread across the unit interval rather than piling up near zero. That is the null {func}`~mantispy.tl.aggregate_guides` matches each gene against. It combines a gene's guides by Stouffer's z, which rewards a consistent shift across the guides: a gene whose guides all move a little is called, a gene carried by one lucky guide is not." ] }, { "cell_type": "code", "execution_count": 3, "id": "cell-5", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T00:50:22.519208Z", "iopub.status.busy": "2026-10-21T00:50:22.519055Z", "iopub.status.idle": "2026-10-21T00:50:22.987843Z", "shell.execute_reply": "2026-10-21T00:50:22.986251Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "118 of 124 genes called 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", "
groupn_guidesstatisticpvalueqvalueis_hit
0TPR103.2908710.000050.000067True
1TPT193.4990750.000050.000067True
2TUBA1C104.6475070.000050.000067True
3TUBB105.7416750.000050.000067True
4TUBG1104.9707300.000050.000067True
5TUBGCP2103.1383070.000050.000067True
6TUBGCP3106.3111090.000050.000067True
7TUBGCP494.6609660.000050.000067True
8TUBGCP5103.8228960.000050.000067True
9TUT195.0064420.000050.000067True
10VCP86.4066760.000050.000067True
11ZBTB17104.0033420.000050.000067True
\n", "
" ], "text/plain": [ " group n_guides statistic pvalue qvalue is_hit\n", "0 TPR 10 3.290871 0.00005 0.000067 True\n", "1 TPT1 9 3.499075 0.00005 0.000067 True\n", "2 TUBA1C 10 4.647507 0.00005 0.000067 True\n", "3 TUBB 10 5.741675 0.00005 0.000067 True\n", "4 TUBG1 10 4.970730 0.00005 0.000067 True\n", "5 TUBGCP2 10 3.138307 0.00005 0.000067 True\n", "6 TUBGCP3 10 6.311109 0.00005 0.000067 True\n", "7 TUBGCP4 9 4.660966 0.00005 0.000067 True\n", "8 TUBGCP5 10 3.822896 0.00005 0.000067 True\n", "9 TUT1 9 5.006442 0.00005 0.000067 True\n", "10 VCP 8 6.406676 0.00005 0.000067 True\n", "11 ZBTB17 10 4.003342 0.00005 0.000067 True" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mt.tl.aggregate_guides(\n", " guides,\n", " score=\"guide_activity\",\n", " guide=\"Metadata_sgRNA\",\n", " gene=\"Metadata_Gene\",\n", " control=\"nontargeting\",\n", " method=\"stouffer\",\n", " n_null=20000,\n", " alpha=0.05,\n", " seed=0,\n", ")\n", "table = guides.uns[\"mantispy\"][\"gene_aggregation\"]\n", "print(f\"{int(table['is_hit'].sum())} of {len(table)} genes called at q < 0.05\")\n", "table.head(12)" ] }, { "cell_type": "markdown", "id": "cell-6", "metadata": {}, "source": [ "The top of the list is the cell's structural machinery: the gamma-tubulin ring complex that nucleates microtubules (`TUBG1`, `TUBGCP2` through `TUBGCP5`), the tubulins themselves (`TUBB`, `TUBA1C`), the nuclear-pore anchor `TPR`, and the AAA-ATPase `VCP`. Knocking any of these out derails the cell in a way a morphology readout cannot miss, so every one of their guides moves, and Stouffer combines them into a tiny p-value." ] }, { "cell_type": "code", "execution_count": 4, "id": "cell-7", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T00:50:22.991877Z", "iopub.status.busy": "2026-10-21T00:50:22.991559Z", "iopub.status.idle": "2026-10-21T00:50:23.005330Z", "shell.execute_reply": "2026-10-21T00:50:23.004297Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "8 of 9 documented genes called\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", "
groupn_guidesqvalueis_hit
12ARPC4100.000067True
16COPE100.000067True
43MRPL43100.000067True
50MRPS5100.000067True
66PSMB180.000067True
83PSMD480.000067True
84ARCN180.000067True
107ACTR6100.003961True
119KIF18A100.114539False
\n", "
" ], "text/plain": [ " group n_guides qvalue is_hit\n", "12 ARPC4 10 0.000067 True\n", "16 COPE 10 0.000067 True\n", "43 MRPL43 10 0.000067 True\n", "50 MRPS5 10 0.000067 True\n", "66 PSMB1 8 0.000067 True\n", "83 PSMD4 8 0.000067 True\n", "84 ARCN1 8 0.000067 True\n", "107 ACTR6 10 0.003961 True\n", "119 KIF18A 10 0.114539 False" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The dataset documents its own known-mechanism genes; check how the call recovers them.\n", "documented = [\"KIF18A\", \"PSMB1\", \"PSMD4\", \"MRPL43\", \"MRPS5\", \"ARPC4\", \"ACTR6\", \"COPE\", \"ARCN1\"]\n", "recovery = table[table[\"group\"].isin(documented)].sort_values(\"qvalue\")\n", "print(f\"{int(recovery['is_hit'].sum())} of {len(documented)} documented genes called\")\n", "recovery[[\"group\", \"n_guides\", \"qvalue\", \"is_hit\"]]" ] }, { "cell_type": "markdown", "id": "cell-8", "metadata": {}, "source": [ "The proteasome, mitochondrial ribosome, ARP2/3 and COPI genes the dataset names all come through, so the call recovers the screen's stated ground truth across unrelated pathways, not only the structural set. `KIF18A` is the one that does not clear the cut: its knockout acts on chromosome congression in mitosis, a phenotype a whole-cell morphology summary captures only weakly, and the activity used here, the distance from the control centre, averages it away. A more specific per-guide score would be needed to call it." ] }, { "cell_type": "code", "execution_count": 5, "id": "cell-9", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T00:50:23.008539Z", "iopub.status.busy": "2026-10-21T00:50:23.008340Z", "iopub.status.idle": "2026-10-21T00:50:23.014443Z", "shell.execute_reply": "2026-10-21T00:50:23.013425Z" } }, "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", "
groupn_guidesstatisticpvalueqvalueis_hit
118DDX11100.4241830.1054950.109927False
119KIF18A100.3994530.1108440.114539False
120CCDC115100.3899190.1127940.115591False
121MBTPS110-0.5821450.4860260.493993False
122MBTPS210-0.8717020.6224690.627530False
123KIF4A10-1.0475660.7016650.701665False
\n", "
" ], "text/plain": [ " group n_guides statistic pvalue qvalue is_hit\n", "118 DDX11 10 0.424183 0.105495 0.109927 False\n", "119 KIF18A 10 0.399453 0.110844 0.114539 False\n", "120 CCDC115 10 0.389919 0.112794 0.115591 False\n", "121 MBTPS1 10 -0.582145 0.486026 0.493993 False\n", "122 MBTPS2 10 -0.871702 0.622469 0.627530 False\n", "123 KIF4A 10 -1.047566 0.701665 0.701665 False" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The genes that do not clear the cut, weakest last.\n", "table.tail(6)" ] }, { "cell_type": "markdown", "id": "cell-10", "metadata": {}, "source": [ "Six genes are not called. `DDX11`, `KIF18A` and `CCDC115` have small positive statistics, while `MBTPS1`, `MBTPS2` and `KIF4A` have negative ones, below the non-targeting spread. The method does not call everything: even a hit-enriched screen leaves a handful of genes uncalled." ] }, { "cell_type": "code", "execution_count": 6, "id": "cell-11", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T00:50:23.017236Z", "iopub.status.busy": "2026-10-21T00:50:23.017061Z", "iopub.status.idle": "2026-10-21T00:50:23.380830Z", "shell.execute_reply": "2026-10-21T00:50:23.379646Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "stouffer hits: 118, fisher hits: 119, shared: 118\n" ] } ], "source": [ "# Fisher instead of Stouffer: combine the guides' p-values, which fires on a single strong guide.\n", "fisher = guides.copy()\n", "mt.tl.aggregate_guides(\n", " fisher,\n", " score=\"guide_activity\",\n", " guide=\"Metadata_sgRNA\",\n", " gene=\"Metadata_Gene\",\n", " control=\"nontargeting\",\n", " method=\"fisher\",\n", " n_null=20000,\n", " alpha=0.05,\n", " seed=0,\n", ")\n", "ft = fisher.uns[\"mantispy\"][\"gene_aggregation\"]\n", "both = set(table.loc[table[\"is_hit\"], \"group\"]) & set(ft.loc[ft[\"is_hit\"], \"group\"])\n", "print(f\"stouffer hits: {int(table['is_hit'].sum())}, fisher hits: {int(ft['is_hit'].sum())}, shared: {len(both)}\")" ] }, { "cell_type": "markdown", "id": "cell-12", "metadata": {}, "source": [ "Stouffer and Fisher agree almost exactly on this screen, because its hits are strong and every guide of a hit gene moves, so the two share all 118 Stouffer hits and Fisher adds one more gene. The two part company on a gene with one potent guide among duds: Fisher calls it, Stouffer does not. Which is right depends on how much a single reagent is trusted, so the method leaves the choice to the caller.\n", "\n", "## Summary\n", "\n", "A pooled screen gives a noisy phenotype per guide, so a gene call pools the guides. {func}`~mantispy.tl.guide_activity` turns each guide into a one-sided activity p-value against a reference control class, and {func}`~mantispy.tl.aggregate_guides` combines a gene's guides against a non-targeting null matched to its guide count. Here they recover the microtubule and nuclear-pore genes and 8 of the 9 documented genes, and leave six genes uncalled." ] } ], "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 }