Calling gene hits in a pooled screen

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.

This case study scores the guides of the cp_posh screen for phenotypic activity with guide_activity(), then calls genes with 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.

import mantispy as mt

guides = mt.ds.cp_posh(aggregated=True)
print(guides)
guides.obs["Metadata_Gene"].value_counts().head(4)

AnnData object with n_obs × n_vars = 1551 × 1278
    obs: 'Metadata_Gene', 'Metadata_sgRNA', 'Metadata_CellCount', 'Metadata_Perturbation', 'Metadata_Plate', 'Metadata_Control', 'Metadata_Perturbation_Type'
    var: 'object', 'feature_group', 'feature', 'channel', 'scale', 'angle', 'gray_levels', 'radial_bin', 'params', 'is_feature'
    uns: 'mantispy'
    layers: None (.X)

Metadata_Gene
nontargeting    200
intergenic      199
ACIN1            10
ACO2             10
Name: count, dtype: int64

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.

# Per-guide activity: how far each guide sits from the intergenic control centre, as a one-sided
# p-value small for an active guide. The intergenic class sets the scale and is left unscored.
mt.tl.guide_activity(guides, reference="intergenic")
scored = guides.obs["guide_activity"].notna().sum()
print(f"scored {scored} guides against the intergenic reference")
guides.obs.loc[guides.obs["Metadata_Gene"] == "nontargeting", "guide_activity"].describe()[["min", "50%", "max"]]

scored 1352 guides against the intergenic reference

min    0.035
50%    0.560
max    0.995
Name: guide_activity, dtype: float64

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 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.

mt.tl.aggregate_guides(
    guides,
    score="guide_activity",
    guide="Metadata_sgRNA",
    gene="Metadata_Gene",
    control="nontargeting",
    method="stouffer",
    n_null=20000,
    alpha=0.05,
    seed=0,
)
table = guides.uns["mantispy"]["gene_aggregation"]
print(f"{int(table['is_hit'].sum())} of {len(table)} genes called at q < 0.05")
table.head(12)

118 of 124 genes called at q < 0.05

	gene	n_guides	statistic	pvalue	qvalue	is_hit
0	TPR	10	3.290871	0.00005	0.000067	True
1	TPT1	9	3.499075	0.00005	0.000067	True
2	TUBA1C	10	4.647507	0.00005	0.000067	True
3	TUBB	10	5.741675	0.00005	0.000067	True
4	TUBG1	10	4.970730	0.00005	0.000067	True
5	TUBGCP2	10	3.138307	0.00005	0.000067	True
6	TUBGCP3	10	6.311109	0.00005	0.000067	True
7	TUBGCP4	9	4.660966	0.00005	0.000067	True
8	TUBGCP5	10	3.822896	0.00005	0.000067	True
9	TUT1	9	5.006442	0.00005	0.000067	True
10	VCP	8	6.406676	0.00005	0.000067	True
11	ZBTB17	10	4.003342	0.00005	0.000067	True

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.

# The dataset documents its own known-mechanism genes; check how the call recovers them.
documented = ["KIF18A", "PSMB1", "PSMD4", "MRPL43", "MRPS5", "ARPC4", "ACTR6", "COPE", "ARCN1"]
recovery = table[table["gene"].isin(documented)].sort_values("qvalue")
print(f"{int(recovery['is_hit'].sum())} of {len(documented)} documented genes called")
recovery[["gene", "n_guides", "qvalue", "is_hit"]]

8 of 9 documented genes called

	gene	n_guides	qvalue	is_hit
12	ARPC4	10	0.000067	True
16	COPE	10	0.000067	True
43	MRPL43	10	0.000067	True
50	MRPS5	10	0.000067	True
66	PSMB1	8	0.000067	True
83	PSMD4	8	0.000067	True
84	ARCN1	8	0.000067	True
107	ACTR6	10	0.003961	True
119	KIF18A	10	0.114539	False

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.

# The genes that do not clear the cut, weakest last.
table.tail(6)

	gene	n_guides	statistic	pvalue	qvalue	is_hit
118	DDX11	10	0.424183	0.105495	0.109927	False
119	KIF18A	10	0.399453	0.110844	0.114539	False
120	CCDC115	10	0.389919	0.112794	0.115591	False
121	MBTPS1	10	-0.582145	0.486026	0.493993	False
122	MBTPS2	10	-0.871702	0.622469	0.627530	False
123	KIF4A	10	-1.047566	0.701665	0.701665	False

The genes that fall out are the ones a general morphology assay is blind to. MBTPS1 and MBTPS2 cut SREBP processing, a lipid-regulatory step with little shape consequence here, and their combined statistic is negative, below the non-targeting spread. The method does not call everything: a hit-enriched screen leaves a handful of genes uncalled, in the right place.

# Fisher instead of Stouffer: combine the guides' p-values, which fires on a single strong guide.
fisher = guides.copy()
mt.tl.aggregate_guides(
    fisher,
    score="guide_activity",
    guide="Metadata_sgRNA",
    gene="Metadata_Gene",
    control="nontargeting",
    method="fisher",
    n_null=20000,
    alpha=0.05,
    seed=0,
)
ft = fisher.uns["mantispy"]["gene_aggregation"]
both = set(table.loc[table["is_hit"], "gene"]) & set(ft.loc[ft["is_hit"], "gene"])
print(f"stouffer hits: {int(table['is_hit'].sum())}, fisher hits: {int(ft['is_hit'].sum())}, shared: {len(both)}")

stouffer hits: 118, fisher hits: 119, shared: 118

Stouffer and Fisher agree almost exactly on this screen, because its hits are strong and every guide of a hit gene moves, so there is no lone-guide gene for Fisher to rescue or Stouffer to miss. 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.

Summary

A pooled screen gives a noisy phenotype per guide, so a gene call pools the guides. guide_activity() turns each guide into a one-sided activity p-value against a reference control class, and aggregate_guides() combines a gene’s guides against a non-targeting null matched to its guide count. Here they recover the microtubule, nuclear-pore and the dataset’s documented genes, and leave the morphologically silent ones uncalled.