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

## ----install, eval=FALSE------------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# 
# BiocManager::install("fishash")

## ----load-packages, message = FALSE-------------------------------------------
library(Matrix)
library(SummarizedExperiment)
library(SingleCellExperiment)

library(ggplot2)
library(dplyr)

library(fishash)

## ----load-data, cache=TRUE----------------------------------------------------
data(tapseq_diffex)
tapseq_diffex

## ----print-altexp-------------------------------------------------------------
altExp(tapseq_diffex)

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

## ----diag-heatmap, fig.width=9, fig.height=8, fig.wide=TRUE-------------------
set.seed(12345)
diagonal_heatmap(counts(altExp(tapseq_diffex)), subsample_cells = 200)

## ----weighted-hist, fig.width=7, fig.height=7---------------------------------
nonzero_histogram_with_weighted(counts(altExp(tapseq_diffex)))

## ----run-fishash--------------------------------------------------------------
res_fishash <- fishash(counts(altExp(tapseq_diffex)), padj_cutoff = .05)
res_fishash

## ----show-assigned------------------------------------------------------------
assay(res_fishash, "assigned")[1:10, 1:10]

## ----show-coldata-------------------------------------------------------------
colData(res_fishash)

## ----bar-demux-type, fig.small=TRUE-------------------------------------------
colData(res_fishash) |>
    as.data.frame() |>
    ggplot(aes(x = demux_type)) +
    geom_bar() +
    theme_bw(base_size = 16)

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

## ----show-logpval-------------------------------------------------------------
assay(res_fishash, "log_pval")[1:10, 1:10]

## ----show-metadata------------------------------------------------------------
metadata(res_fishash)

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

## ----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"
]

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

## ----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")
    )
}

## ----session-info-------------------------------------------------------------
sessionInfo()