Using fishash to assign gRNAs in Perturbseq data
Revised: 2026-07-25
Package
fishash 0.99.3
Contents
1 Introduction
Fishash is a package for assigning guides to cells in Perturbseq experiments, which combine single-cell RNAseq with CRISPR perturbations. During preprocessing, we must determine the presence/absence of guides in each cell based on the CRISPR gRNA UMI counts. Fishash does this using a modified version of Fisher’s test.
Other packages for performing perturbseq guide assignment include:
We compare fishash with the methods above in our preprint (Kamm, Yeung, and Forrest 2026). We find fishash to be particularly well suited for large perturbseq datasets (e.g. genome-wide perturbation screens) in terms of accuracy, runtime, and memory usage.
2 Installation
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("fishash")3 Example TAPseq dataset
library(Matrix)
library(SummarizedExperiment)
library(SingleCellExperiment)
library(ggplot2)
library(dplyr)
library(fishash)We illustrate fishash on a Targeted Perturbseq (“TAPseq”) dataset (Schraivogel et al. 2020) dataset which is also bundled with the crispat package for gRNA assignment.
The data consists of ~22,000 K562 cells from 2 samples with gene expression for 72 genes (mostly on chromosome 8) and CRISPRi perturbations consisting of 86 gRNAs. Note that TAPseq amplifies transcripts from a panel of genes instead of sequencing the whole transcriptome.
We represent the data as a SingleCellExperiment, with the
main experiment containing the gene expression sequencing counts.
data(tapseq_diffex)
tapseq_diffex
#> class: SingleCellExperiment
#> dim: 72 21977
#> metadata(0):
#> assays(1): counts
#> rownames(72): ANGPT1 ANKRD46 ... ZHX1 ZNF706
#> rowData names(0):
#> colnames(21977): TAP1-AAACCTGAGCCTATGT TAP1-AAACCTGAGGCGTACA ...
#> TAP2-TTTGTCATCGCAGGCT TAP2-TTTGTCATCTACCTGC
#> colData names(2): cell_id sample
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(1): grnaThe SingleCellExperiment object also contains an alternative experiment with the sequencing counts of the gRNA barcodes, which we will use to call presence/absence of the gRNAs:
altExp(tapseq_diffex)
#> class: SingleCellExperiment
#> dim: 86 21977
#> metadata(0):
#> assays(1): counts
#> rownames(86): CCNE2-95907328 CCNE2-95907382 ... NTC-28 NTC-29
#> rowData names(4): grna_id grna_name target_gene target_element_type
#> colnames(21977): TAP1-AAACCTGAGCCTATGT TAP1-AAACCTGAGGCGTACA ...
#> TAP2-TTTGTCATCGCAGGCT TAP2-TTTGTCATCTACCTGC
#> colData names(0):
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):The CRISPRi perturbations target promoters of 10 genes and enhancers of 4 genes. Each gene target has 4 guides targeting its promoter or enhancers. Additionally, there are 30 non-targeting control (NTC) guides.
rowData(altExp(tapseq_diffex)) |>
as.data.frame() |>
group_by(target_gene, target_element_type) |>
summarize(n_guides = n(), .groups = "drop") |>
arrange(target_element_type)
#> # A tibble: 15 × 3
#> target_gene target_element_type n_guides
#> <chr> <chr> <int>
#> 1 NTC NTC 30
#> 2 GATA1 enhancer 4
#> 3 HS2 enhancer 4
#> 4 MYC enhancer 4
#> 5 ZFPM2 enhancer 4
#> 6 CCNE2 promoter 4
#> 7 CPQ promoter 4
#> 8 DSCC1 promoter 4
#> 9 FAM83A promoter 4
#> 10 LRRCC1 promoter 4
#> 11 OXR1 promoter 4
#> 12 PHF20L1 promoter 4
#> 13 RIPK2 promoter 4
#> 14 STK3 promoter 4
#> 15 UBR5 promoter 44 Visualizing the guide counts
Before assigning the guides, it is good to visualize the count matrix.
One useful plot, originally proposed by the pertpy project, is to represent the counts as a heatmap, with the cells sorted by their top guide, creating a band near the diagonal:
set.seed(12345)
diagonal_heatmap(counts(altExp(tapseq_diffex)), subsample_cells = 200)
Another useful plot is the histogram of the nonzero counts. We plot the histogram together with a weighted version of the histogram, which rescales the y-axis to represent UMIs instead of the number of matrix entries:
nonzero_histogram_with_weighted(counts(altExp(tapseq_diffex)))
The top histogram represents the distribution for the number of UMIs if we sampled a random nonzero matrix entry, whereas the bottom histogram represents the distribution if we sampled a random UMI and asked how many UMIs are in the same entry. If we think of the left mode as “noise” and the right mode as “signal”, the areas in the weighted histogram provide a visual estimate of the signal-to-noise ratio in the gRNA counts.
5 Guide assignment
Next, we call fishash() to call presence/absence of the guides in each cell, based on the sequencing counts of the gRNA barcodes. Internally, for each (cell, gRNA) pair the function constructs a 2x2 table for the UMI counts, and then uses a one-sided Fisher’s test to call presence/absence of the guide, depending on whether there is significant association of the cell and guide barcodes.
res_fishash <- fishash(counts(altExp(tapseq_diffex)), padj_cutoff = .05)
res_fishash
#> class: SummarizedExperiment
#> dim: 86 21977
#> metadata(2): log_pval_cutoff num_iter
#> assays(4): assigned log_pval odds_ratio odds_ratio_regularized
#> rownames(86): CCNE2-95907328 CCNE2-95907382 ... NTC-28 NTC-29
#> rowData names(0):
#> colnames(21977): TAP1-AAACCTGAGCCTATGT TAP1-AAACCTGAGGCGTACA ...
#> TAP2-TTTGTCATCGCAGGCT TAP2-TTTGTCATCTACCTGC
#> colData names(2): demux_type assignmentThe return value is a SummarizedExperiment. The assigned assay is a
sparse boolean matrix indicating whether we call the guide present
(| means TRUE while . means FALSE):
assay(res_fishash, "assigned")[1:10, 1:10]
#> 10 x 10 sparse Matrix of class "lgCMatrix"
#> [[ suppressing 10 column names 'TAP1-AAACCTGAGCCTATGT', 'TAP1-AAACCTGAGGCGTACA', 'TAP1-AAACCTGAGGGCTTGA' ... ]]
#>
#> CCNE2-95907328 . . . . . . . . . .
#> CCNE2-95907382 . . . . . . . . . .
#> CCNE2-95907406 . . . . . . . . . .
#> CCNE2-95907017 . . . . . . . . . .
#> CPQ-97657557 . | . . . . . . . .
#> CPQ-97657573 . . . . . . . . . .
#> CPQ-97657567 . . . . . . . . . .
#> CPQ-97657591 . . . . . . . . . .
#> DSCC1-120867694 . . . . . . . . . .
#> DSCC1-120868042 . . . . . . . . . .The colData shows the classification of each cell; the assignment
column is a comma-delimited string containing the guides called, while
the demux_type column indicates whether the cell contained 0, 1, or
2+ assigned guides (only the cells with a single guide will be kept
for downstream analysis):
colData(res_fishash)
#> DataFrame with 21977 rows and 2 columns
#> demux_type assignment
#> <character> <character>
#> TAP1-AAACCTGAGCCTATGT doublet NTC-25,NTC-27
#> TAP1-AAACCTGAGGCGTACA singlet CPQ-97657557
#> TAP1-AAACCTGAGGGCTTGA singlet MYC-D
#> TAP1-AAACCTGAGGGTGTTG singlet PHF20L1-133787636
#> TAP1-AAACCTGAGGTGCTAG singlet LRRCC1-86019483
#> ... ... ...
#> TAP2-TTTGTCAGTTGTACAC singlet NTC-10
#> TAP2-TTTGTCATCAAGATCC singlet NTC-28
#> TAP2-TTTGTCATCCGTCAAA singlet OXR1-107670159
#> TAP2-TTTGTCATCGCAGGCT singlet NTC-29
#> TAP2-TTTGTCATCTACCTGC doublet LRRCC1-86019818,UBR5..We inspect the number of cells with 0 guides (unassigned), 1 guide (singlet), or 2+ guides (doublet):
colData(res_fishash) |>
as.data.frame() |>
ggplot(aes(x = demux_type)) +
geom_bar() +
theme_bw(base_size = 16)
Alternatively, we can plot a histogram of the number of assigned guides per cell by using the column sums of the assignment matrix:
data.frame(
n_assigned = colSums(assay(res_fishash, "assigned"))
) |>
ggplot(aes(x = n_assigned)) +
geom_histogram() +
theme_bw(base_size = 16)
#> `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
We can also inspect the log p-values from the Fisher test. Since we use a one-sided test, this matrix is also sparse:
assay(res_fishash, "log_pval")[1:10, 1:10]
#> 10 x 10 sparse Matrix of class "dgCMatrix"
#> [[ suppressing 10 column names 'TAP1-AAACCTGAGCCTATGT', 'TAP1-AAACCTGAGGCGTACA', 'TAP1-AAACCTGAGGGCTTGA' ... ]]
#>
#> CCNE2-95907328 . . . . . . . . . .
#> CCNE2-95907382 . . . . . . . . . .
#> CCNE2-95907406 . . . . . . . . . .
#> CCNE2-95907017 . . . . . . . . . .
#> CPQ-97657557 . -118.8247 . . . . . . . .
#> CPQ-97657573 . . . . . . . . . .
#> CPQ-97657567 . . . . . . . . -0.7195586 .
#> CPQ-97657591 . . . . . . . . . .
#> DSCC1-120867694 . . . . . . . . . .
#> DSCC1-120868042 . . . . . . . . . .The metadata shows the p-value cutoff selected by the FDR procedure:
metadata(res_fishash)
#> $log_pval_cutoff
#> [1] -7.466158
#>
#> $num_iter
#> [1] 4We can inspect the classifier decision rule by plotting the test statistic (the negative log p-value). Ideally its histogram will look bimodal, as it does below. We also plot the test statistic against the UMI count to understand how they correspond; we see that guides with 2 counts are usually not assigned, and guides with 10+ counts are nearly always assigned.
nz <- counts(altExp(tapseq_diffex)) > 0
data.frame(
count = counts(altExp(tapseq_diffex))[nz],
log_pval = assay(res_fishash, "log_pval")[nz],
assigned = assay(res_fishash, "assigned")[nz]
) |>
ggplot(aes(x = -log_pval, y = count, color = assigned)) +
geom_point(
shape = 1, alpha = .5,
position = position_jitter(width = 0, height = .02)
) +
geom_vline(
xintercept = -metadata(res_fishash)$log_pval_cutoff,
lty = "dotted"
) +
scale_x_continuous(trans = "log1p", breaks = c(0, 10^(0:6))) +
scale_y_log10() +
theme_bw(base_size = 16) +
theme(legend.position = "left") ->
p
if (require(ggExtra)) {
p <- ggMarginal(
p,
type = "hist",
groupFill = TRUE,
xparams = list(bins = 50),
yparams = list(bins = 50)
)
}
#> Loading required package: ggExtra
p
6 Downstream analysis
Finally, we briefly illustrate how the gRNA assignments may be used in downstream analyses such as differential expression.
First, we subset to cells with a single guide detected.
# Add fishash results to the SingleCellExperiment
colData(tapseq_diffex) <- cbind(colData(tapseq_diffex), colData(res_fishash))
assay(altExp(tapseq_diffex), "assigned") <- assay(res_fishash, "assigned")
# Subset to cells with a single guide
singlets <- tapseq_diffex[, colData(tapseq_diffex)$demux_type == "singlet"]
# Annotate cells with the perturbed gene
colData(singlets)$target_gene <- rowData(altExp(singlets))[
colData(singlets)$assignment, "target_gene"
]Next, we use glmGamPoi to perform differential expression of the cells from each guide against the NTC cells, controlling for the sample. For the NTC guides, we compare them against the other NTC guides in leave-one-out fashion. For the cell size factors, we use the geometric mean (CLR normalization) due to the compositional nature of the data and low number of features (72 genes). We also shrink the log-fold-changes with a small ridge penalty equivalent to a Normal prior with variance 1.
if (require(glmGamPoi)) {
# Run differential expression in a loop over guides
lfc_prior_sd <- 1
de_list <- lapply(
rownames(altExp(singlets)),
function(g) {
# Subset to NTCs or cells with guide g
keep <- with(
colData(singlets),
target_gene == "NTC" | assignment == g
)
sce_sub <- singlets[, keep]
colData(sce_sub)$prtrb <- colData(sce_sub)$assignment == g
# Run differential expression
test_de(
glm_gp(
counts(sce_sub),
design = ~ sample + prtrb,
col_data = colData(sce_sub),
size_factors = exp(colMeans(log1p(counts(sce_sub)))),
ridge_penalty = log(2) / lfc_prior_sd / sqrt(ncol(sce_sub)),
overdispersion_shrinkage = FALSE
),
contrast = "prtrbTRUE"
)
}
)
names(de_list) <- rownames(altExp(singlets))
# convert differential expression results to SummarizedExperiment
de_summexp <- list()
for (a in c("lfc", "adj_pval")) {
de_summexp[[a]] <- sapply(
de_list,
function(x) setNames(x[, a], x[, "name"])
)
}
de_summexp <- SummarizedExperiment(
assays = de_summexp,
colData = rowData(altExp(singlets))
)
} else {
de_summexp <- NULL
}
#> Loading required package: glmGamPoi
#>
#> Attaching package: 'glmGamPoi'
#> The following object is masked from 'package:dplyr':
#>
#> vars
#> The following object is masked from 'package:ggplot2':
#>
#> varsWe then plot a heatmap of the log-fold changes; significant genes at FDR 5% are marked with a point.
if (!is.null(de_summexp) && require(ComplexHeatmap)) {
Heatmap(
t(assay(de_summexp, "lfc")),
right_annotation = rowAnnotation(
df = colData(de_summexp)[, "target_element_type", drop = FALSE],
col = list(target_element_type = setNames(
palette()[1:3],
c("NTC", "enhancer", "promoter")
))
),
layer_fun = function(j, i, x, y, w, h, col) {
signif <- t(assay(de_summexp, "adj_pval"))[cbind(i, j)] <= .05
grid.points(
x[signif], y[signif],
pch = 21, size = unit(3, "mm"),
gp = gpar(col = "white", fill = "black", lwd = 2)
)
},
cluster_rows = FALSE, cluster_columns = FALSE,
heatmap_legend_param = list(title = "logFC")
)
}
#> Loading required package: ComplexHeatmap
#> Loading required package: grid
#> ========================================
#> ComplexHeatmap version 2.29.0
#> Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
#> Github page: https://github.com/jokergoo/ComplexHeatmap
#> Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
#>
#> If you use it in published research, please cite either one:
#> - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
#> - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
#> genomic data. Bioinformatics 2016.
#>
#>
#> The new InteractiveComplexHeatmap package can directly export static
#> complex heatmaps into an interactive Shiny app with zero effort. Have a try!
#>
#> This message can be suppressed by:
#> suppressPackageStartupMessages(library(ComplexHeatmap))
#> ========================================
#> The automatically generated colors map from the minus and plus 99^th of
#> the absolute values in the matrix. There are outliers in the matrix
#> whose patterns might be hidden by this color mapping. You can manually
#> set the color to `col` argument.
#>
#> Use `suppressMessages()` to turn off this message.
Overall the results look reasonable, with similar differential expression results between the guides for the same target, and in particular with the target gene usually being knocked down (e.g. ZFPM2 guides have negative logFC for ZFPM2 gene expression). Some guides should be investigated as outliers – for example, the fourth UBR5 guide shares few DEGs with the other UBR5 guides and in particular fails to knockdown the UBR5 gene.
The number of DEGs among the NTCs can also give a rough idea of the empirical false positive rate or miscalibration of the results. Overall, most NTC guides have no DEGs, with some exceptions (such as NTC-22) which may be worth investigating as outliers.
Beyond the scope of this vignette, a more detailed analysis could improve the quality of the DE results, such as by QC filtering of cells, by controlling for confounders (such as cell state), by using more sophisticated shrinkage estimators of the log-fold-changes, or by improving estimates of the cell size factors or gene overdispersions.
References
Kamm, Jack, Jake Yeung, and Bill Forrest. 2026. “Fishash: A Contingency Table Approach to Perturb-Seq Guide Assignment.” bioRxiv, 2026–01.
Schraivogel, Daniel, Andreas R Gschwind, Jennifer H Milbank, Daniel R Leonce, Petra Jakob, Lukas Mathur, Jan O Korbel, Christoph A Merten, Lars Velten, and Lars M Steinmetz. 2020. “Targeted Perturb-Seq Enables Genome-Scale Genetic Screens in Single Cells.” Nature Methods 17 (6): 629–35.
Appendix
A Session info
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] grid stats4 stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] ComplexHeatmap_2.29.0 glmGamPoi_1.25.1
#> [3] ggExtra_0.11.0 fishash_0.99.3
#> [5] dplyr_1.2.1 ggplot2_4.0.3
#> [7] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
#> [9] Biobase_2.73.2 GenomicRanges_1.65.1
#> [11] Seqinfo_1.3.0 IRanges_2.47.2
#> [13] S4Vectors_0.51.6 BiocGenerics_0.59.11
#> [15] generics_0.1.4 MatrixGenerics_1.25.0
#> [17] matrixStats_1.5.0 Matrix_1.7-6
#> [19] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 farver_2.1.2
#> [3] S7_0.2.2 fastmap_1.2.0
#> [5] promises_1.5.0 digest_0.6.39
#> [7] mime_0.13 lifecycle_1.0.5
#> [9] cluster_2.1.8.3 Cairo_1.7-0
#> [11] magrittr_2.0.5 compiler_4.6.1
#> [13] rlang_1.3.0 sass_0.4.10
#> [15] tools_4.6.1 utf8_1.2.6
#> [17] yaml_2.3.12 knitr_1.51
#> [19] S4Arrays_1.13.0 labeling_0.4.3
#> [21] DelayedArray_0.39.4 RColorBrewer_1.1-3
#> [23] abind_1.4-8 miniUI_0.1.2
#> [25] withr_3.0.3 nnet_7.3-21
#> [27] beachmat_2.29.0 colorspace_2.1-3
#> [29] xtable_1.8-8 scales_1.4.0
#> [31] iterators_1.0.14 dichromat_2.0-1
#> [33] tinytex_0.60 cli_3.6.6
#> [35] crayon_1.5.3 rmarkdown_2.31
#> [37] otel_0.2.0 rjson_0.2.23
#> [39] DelayedMatrixStats_1.35.0 cachem_1.1.0
#> [41] parallel_4.6.1 BiocManager_1.30.27
#> [43] XVector_0.53.0 vctrs_0.7.3
#> [45] jsonlite_2.0.0 bookdown_0.47
#> [47] GetoptLong_1.1.1 patchwork_1.3.2
#> [49] clue_0.3-68 magick_2.9.1
#> [51] foreach_1.5.2 jquerylib_0.1.4
#> [53] glue_1.8.1 codetools_0.2-20
#> [55] shape_1.4.6.1 gtable_0.3.6
#> [57] later_1.4.8 tibble_3.3.1
#> [59] pillar_1.11.1 htmltools_0.5.9
#> [61] circlize_0.4.18 R6_2.6.1
#> [63] sparseMatrixStats_1.25.0 doParallel_1.0.17
#> [65] evaluate_1.0.5 shiny_1.14.0
#> [67] lattice_0.22-9 extraDistr_1.10.0.5
#> [69] png_0.1-9 httpuv_1.6.17
#> [71] bslib_0.12.0 Rcpp_1.1.2
#> [73] SparseArray_1.13.2 xfun_0.60
#> [75] GlobalOptions_0.1.4 pkgconfig_2.0.3