---
title: Using fishash to assign gRNAs in Perturbseq data
author: Jack Kamm
date: "Revised: 2026-07-25"
output:
  BiocStyle::html_document
package: fishash
bibliography: ref.bib
vignette: >
  %\VignetteIndexEntry{Using fishash to assign gRNAs in Perturbseq data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup-opts, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)
```

# 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:

- R packages
  - `r Githubpkg("katsevich-lab/sceptre")`
- Python packages
  - [geomux](https://github.com/noamteyssier/geomux)
  - [crispat](https://github.com/velten-group/crispat)
  - [demuxEM](https://github.com/lilab-bcb/demuxEM)
- Command line tools
  - [CLEANSER](https://github.com/Gersbachlab-Bioinformatics/CLEANSER)

We compare fishash with the methods above in our preprint
[@kamm2026fishash].  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.

# Installation

```{r install, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install("fishash")
```

# Example TAPseq dataset

```{r load-packages, message = FALSE}
library(Matrix)
library(SummarizedExperiment)
library(SingleCellExperiment)

library(ggplot2)
library(dplyr)

library(fishash)
```

We illustrate fishash on a Targeted Perturbseq ("TAPseq") dataset
[@schraivogel2020targeted] dataset which is also bundled with the
[crispat](https://github.com/velten-group/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.

```{r load-data, cache=TRUE}
data(tapseq_diffex)
tapseq_diffex
```

The 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:

```{r print-altexp}
altExp(tapseq_diffex)
```

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.

```{r table-grna}
rowData(altExp(tapseq_diffex)) |>
    as.data.frame() |>
    group_by(target_gene, target_element_type) |>
    summarize(n_guides = n(), .groups = "drop") |>
    arrange(target_element_type)
```

# Visualizing the guide counts

Before assigning the guides, it is good to visualize the count matrix.

One useful plot, originally proposed by the
[pertpy](https://pertpy.readthedocs.io/en/stable/tutorials/notebooks/guide_rna_assignment.html)
project,
is to represent the counts as a heatmap, with the cells sorted by
their top guide, creating a band near the diagonal:

```{r diag-heatmap, fig.width=9, fig.height=8, fig.wide=TRUE}
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:

```{r weighted-hist, fig.width=7, fig.height=7}
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.

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

```{r run-fishash}
res_fishash <- fishash(counts(altExp(tapseq_diffex)), padj_cutoff = .05)
res_fishash
```

The 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`):

```{r show-assigned}
assay(res_fishash, "assigned")[1:10, 1:10]
```

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):

```{r show-coldata}
colData(res_fishash)
```

We inspect the number of cells with 0 guides (unassigned), 1 guide (singlet),
or 2+ guides (doublet):

```{r bar-demux-type, fig.small=TRUE}
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:

```{r hist-n-assigned}
data.frame(
    n_assigned = colSums(assay(res_fishash, "assigned"))
) |>
    ggplot(aes(x = n_assigned)) +
    geom_histogram() +
    theme_bw(base_size = 16)
```

We can also inspect the log p-values from the Fisher test. Since we
use a one-sided test, this matrix is also sparse:

```{r show-logpval}
assay(res_fishash, "log_pval")[1:10, 1:10]
```

The metadata shows the p-value cutoff selected by the FDR procedure:

```{r show-metadata}
metadata(res_fishash)
```

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

```{r xy-pvalues-cnt}
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)
    )
}

p
```

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

```{r subset-singlet, cache=TRUE}
# 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 `r Biocpkg("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.

```{r run-glmgp, cache=TRUE}
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
}
```

We then plot a heatmap of the log-fold changes; significant genes at FDR 5% are marked with a point.

```{r de-heatmap, fig.height=15, fig.width=15, fig.wide=TRUE}
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")
    )
}
```

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 {.unnumbered}

# Session info

```{r session-info}
sessionInfo()
```