---
title: "<span style='font-size: 16px; color: grey'> crisprQTL proof-of-concept experiment in T cells </span> <br> <span style='font-size: 36px'> gRNA calling </span>"
author: "Ximena Ibarra-Soria"
date: '`r format(Sys.time(), "%d %B, %Y")`'
output:
  html_document:
    theme: paper
    keep_md: true
    fig_width: 5
    fig_height: 5
    fig_caption: yes
    df_print: paged
    code_folding: hide
    toc: true
    toc_depth: 3
    toc_float: 
      collapsed: false
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(scater)
library(scran)
library(data.table)
library(Matrix)
library(BiocParallel)
library(BiocSingular)
library(ggplot2)
library(ggpubr)
library(RColorBrewer)
library(BuenColors)
library(viridis)
library(inlmisc)
library(iSEE)

dir <- paste0(getwd(), "/")

th <- theme_bw() + theme(
  axis.text.x = element_text(size=10), 
  axis.title.x = element_text(size=12), 
  axis.text.y = element_text(size=12), 
  axis.title.y = element_text(size=12),
  panel.grid.major = element_blank(), 
  panel.grid.minor = element_blank(), 
  axis.line = element_line(colour = "black"), 
  panel.border = element_blank(), 
  plot.title = element_text(face="bold", hjust = 0.5))

## configure parallelisation options
bpp <- MulticoreParam(38)
```

We start from the set of ~250K good-quality cells retained after quality control. 

```{r data}
sce <- readRDS(paste0(dir, "results/01_sce_POC_Tcells.goodQual.RAW.Rds"))

c(n_cells=ncol(sce))
```

As the next step, we need to determine which gRNAs are present in each cell.

### gRNA detection in cDNA library

Before performing gRNA calling, we check the concordance in gRNA transcript detection from the cDNA and gRNA libraries, to ensure the PCR enrichment process doesn't introduce significant biases. The cDNA data was mapped against a modified reference genome containing the gRNA sequences from the plasmid library, allowing us to recover reads from gRNA transcripts directly.

We compare the readouts from both the cDNA and gRNA fractions. As expected, the vast majority of cells have no detection in the cDNA fraction. However, the majority of the cells with detection in both libraries identify the same gRNA as the most abundant one. The discordant calls are an upper bound, since some of these will be concordant after accounting for multiple gRNAs per cell (currently only looking at whether the highest detected is the same in both libraries).

```{r compare_grna_counts, fig.width=8, fig.height=5}
## matched counts for gRNAs in the `gRNA` and `cDNA` fractions in the two altExp slots
## find which is the most abundant guide in each
## if there are zero counts, 1 is returned by which.max
df <- data.frame(top_gRNA = apply(counts(altExp(sce, 'CRISPR')), 2, which.max),
                 top_cDNA = apply(counts(altExp(sce, 'gRNA_in_cDNA')), 2, which.max))
df$total_gRNA <- colSums(counts(altExp(sce, 'CRISPR')))
df$total_cDNA <- colSums(counts(altExp(sce, 'gRNA_in_cDNA')))
df$class <- ifelse(df$top_cDNA==1 & df$total_cDNA==0, "not_detected_cDNA",
                   ifelse(df$top_gRNA==1 & df$top_gRNA==0, "not_detected_gRNA",
                          ifelse(df$top_cDNA == df$top_gRNA, "concordant", "discordant")))
df[df$class=="not_detected_cDNA",]$top_cDNA <- -10

plots <- list()
plots[[1]] <- ggplot(df, aes(top_gRNA, top_cDNA, colour=class)) +
  geom_point(size=1, alpha = 0.1) +
  scale_color_manual(values = as.character(GetColors(3, scheme = "light"))) +
  xlab("gRNAs 1 --> 355 - gRNA library") +
  ylab("gRNAs 1 --> 355 - cDNA library") +
  guides(colour = guide_legend(override.aes = list(alpha=1))) +
  th + theme(axis.text.x = element_blank(),
             axis.text.y = element_blank(),
             axis.ticks = element_blank())

plots[[2]] <- ggplot(df, aes(1, fill=class)) +
  geom_bar() +
  scale_fill_manual(values = as.character(GetColors(3, scheme = "light"))) +
  xlab("") + ylab("number of cells") +
  annotate("text", x=1, y=c(0.69e5, 1.44e5, 2e5), 
           label=c("not_detected_cDNA", 
                   "discordant", "concordant"), size=4) +
  th + theme(axis.text.x = element_blank(),
             axis.ticks.x = element_blank())

ggarrange(plotlist = plots, ncol=2, widths = c(0.65, 0.35), 
          legend="none", align = "h")
```

```{r concordance}
round(prop.table(table(df$class))*100, 2)
```

This indicates that the PCR amplification step doesn't introduce spurious signal, but it greatly enhances gRNA transcript detection.

### gRNA calling 

We use a binomial test to determine which gRNAs are likely present in each cell. 

```{r calling_function}
call_gRNAs_binom <- function(sce=NULL, pval_threshold=0.001, 
                             max_to_test=5, total_guides=NULL,
                             prob=NULL,
                             min_count=3, min_fraction=0.05, 
                             plot=TRUE, res=200,
                             ann=NULL, columns=NULL, add_stats=TRUE){
  if(is.null(total_guides)){
    stop("Please provide the total number of gRNAs to use in the binomial test. 
         This should be the total number of gRNAs present in the library.")
  }
  
  ## create data table with detected gRNA-cell combinations
  gRNA_calls <- reshape2::melt(as.matrix(counts(altExp(sce, 'CRISPR'))))
  ## remove guide-cell pairs with no counts
  gRNA_calls <- gRNA_calls[gRNA_calls$value>0,]
  colnames(gRNA_calls) <- c("gRNA", "cell", "UMI_count")
  
  ## remove guides with too few UMIs
  gRNA_calls <- gRNA_calls[gRNA_calls$UMI_count > min_count,]
  
  ## annotate the targets associated with each gRNA
  if(!is.null(ann)){
    if(length(setdiff(columns, colnames(ann)))>0){
      stop(paste("Column(s)", setdiff(columns, colnames(ann)), "are not present in annotation data.frame"))
    }
    for(i in seq_along(columns)){
      gRNA_calls[,columns[i]] <- ann[match(gRNA_calls$gRNA, ann$ID), columns[i]]
    }
  }
  
  ## add total UMIs per cell, to restrict analysis to the largest counts per cell
  totals <- colSums(counts(altExp(sce, 'CRISPR')))
  gRNA_calls$cell_total_UMIs <- totals[gRNA_calls$cell]
  # compute % of total UMIs taken by each gRNA
  gRNA_calls$fraction_total_UMIs <- gRNA_calls$UMI_count/gRNA_calls$cell_total_UMIs
  
  ## convert to data.table
  gRNA_calls <- as.data.table(gRNA_calls)
  
  ## order decreasingly by UMI_count
  setorder(gRNA_calls, -UMI_count)
  # add column with order per cell based on UMI abundance
  gRNA_calls[, order := 1:.N, by = c("cell")] 
  
  # test top gRNAs per cell as determined by max_to_test
  # if more than max_to_test integrations are true, the rest won't be called
  # higher number in max_to_test results in higher computation time
  
  ## use a binomial test to determine if observed counts are greater than background
  # if no probabilitis provided, use equal probabilities
  if(is.null(prob)){
    gRNA_calls$prob <- 1/total_guides
  }else{
    # add proportions of guides in guide library to gRNA_calls
    # provided as a named vector, with names being gRNA ids
    prob <- as.data.table(prob, keep.rownames = TRUE)
    gRNA_calls <- merge(gRNA_calls, prob, by.x="gRNA", by.y="rn")
  }
  gRNA_calls[order <= max_to_test,
             binom := binom.test(UMI_count, cell_total_UMIs, prob, "greater")$p.value, 
             by = seq_len(nrow(gRNA_calls[order <= max_to_test]))]
  # bonferroni correction, per cell
  gRNA_calls[, p_adj := binom * .N,  by = "cell"] 
  gRNA_calls[p_adj > 1, p_adj := 1] # reset pvals that are larger than 1 after correction
  
  # count significant results
  gRNA_calls[, n_detected_gRNAs := length(which(p_adj <= pval_threshold)), by = "cell"]
  
  ## check if there are cells with at least max_to_test called hits
  # if so, retest a larger number
  while(max(gRNA_calls$n_detected_gRNAs) >= max_to_test){
    # get cells that need higher threshold
    done <- gRNA_calls[gRNA_calls$n_detected_gRNAs < max_to_test,]
    missing <- gRNA_calls[gRNA_calls$n_detected_gRNAs >= max_to_test,]
    
    # compute pvals for those cells with twice as many cells tested
    max_to_test <- max_to_test*2
    missing[order <= max_to_test,
            binom := binom.test(UMI_count, cell_total_UMIs, (1 / total_guides), "greater")$p.value, 
            by = seq_len(nrow(missing[order <= max_to_test]))]
    missing[, p_adj := binom * .N,  by = "cell"] 
    missing[p_adj > 1, p_adj := 1]
    
    # join with previous calls
    gRNA_calls <- rbind(done, missing)
    gRNA_calls[, n_detected_gRNAs := length(which(p_adj <= pval_threshold)), by = "cell"]
  }
  
  # for unique calls, also remove those where the fraction from total UMIs is too low
  gRNA_calls[gRNA_calls$n_detected_gRNAs == 1 &
               gRNA_calls$fraction_total_UMIs <= min_fraction,]$p_adj <- 1
  # count significant results again after adjustments above
  gRNA_calls[, n_detected_gRNAs := length(which(p_adj <= pval_threshold)), by = "cell"]
  
  # assign call based on number of significant gRNAs
  gRNA_calls$call <- ifelse(gRNA_calls$n_detected_gRNAs == 0, "unassigned",
                            ifelse(gRNA_calls$n_detected_gRNAs == 1, "unique", "multiple"))
  
  # annotate significant results
  gRNA_calls[, detected_gRNAs := .(paste0(gRNA[which(p_adj <= pval_threshold)], collapse="|")), by = "cell"]
  gRNA_calls[, detected_UMIs := .(paste0(UMI_count[which(p_adj <= pval_threshold)], collapse="|")), by = "cell"]
  gRNA_calls[, detected_fraction := .(paste0(fraction_total_UMIs[which(p_adj <= pval_threshold)], collapse="|")), by = "cell"]
  gRNA_calls[, detected_padj := .(paste0(p_adj[which(p_adj <= pval_threshold)], collapse="|")), by = "cell"]
  
  ## add annotation if required
  if(!is.null(ann)){
    for(i in seq_along(columns)){
      gRNA_calls[, tmp := .(paste0(get(columns[i])[which(p_adj <= pval_threshold)], collapse="|")), by = "cell"]
      colnames(gRNA_calls)[ncol(gRNA_calls)] <- paste0('detected_', columns[i])
    }
  }
  
  ## QC plot
  if(plot){
    ## remove overlapping points
    idx <- subsetPointsByGrid(X = gRNA_calls$fraction_total_UMIs,
                              Y = log10(gRNA_calls$UMI_count),
                              grouping = gRNA_calls$col,
                              resolution = res)
    df <- gRNA_calls[idx,]
    df$call <- factor(df$call, levels=c("unassigned", "multiple", "unique"))
    df <- df[order(df$call),]
    
    ## assign colors to unique/multiple/unassigned
    cols <- c(multiple="indianred3", unassigned="grey", unique="grey10")
    df$col <- cols[match(df$call, names(cols))]
    
    ## plot 
    p <- ggplot(df, aes(fraction_total_UMIs, log10(UMI_count), colour=call)) +
      scale_color_manual(values = cols) +
      geom_point(size=0.25) +
      xlab("% of total UMIs in cell") +
      ylab(expression('log'[10]*' UMIs in gRNA')) +
      geom_hline(yintercept = log10(min_count), lty=2, col="grey10") +
      geom_vline(xintercept = min_fraction, lty=2, col="grey10") +
      th
  }
  
  ## produce sparse matrix of cell-gRNA pairs that are significant
  # subset to significant
  gRNA_calls_sig <- gRNA_calls[gRNA_calls$p_adj < pval_threshold,]
  # gRNA indices
  rows <- match(gRNA_calls_sig$gRNA, rownames(altExp(sce, 'CRISPR')))
  # cell indices
  cols <- match(gRNA_calls_sig$cell, colnames(altExp(sce, 'CRISPR')))
  # create sparse matrix
  gRNA_cell_matrix <- sparseMatrix(i = rows, j = cols, dims = dim(altExp(sce, 'CRISPR')))
  row.names(gRNA_cell_matrix) <- row.names(altExp(sce, 'CRISPR'))
  colnames(gRNA_cell_matrix) <- colnames(altExp(sce, 'CRISPR'))
  # add to sce obect, copying over col/rowData from CRISPR counts
  altExp(sce, 'gRNA_calls') <- altExp(sce, 'CRISPR')
  assay(altExp(sce, 'gRNA_calls')) <- gRNA_cell_matrix
  
  ## produce matrix of assignments
  gRNA_calls <- gRNA_calls[order == 1, ] # keep only one entry per cell
  gRNA_calls <- gRNA_calls[, c("gRNA", "UMI_count", "cell_total_UMIs", "fraction_total_UMIs", "order", "binom", "p_adj"):=NULL]
  if(!is.null(columns)){
    gRNA_calls <- gRNA_calls[, (columns):=NULL]
  }
  ## add to sce colData
  if(add_stats == FALSE){
    # only add number of gRNAs called, and call annotation
    sel <- c("n_detected_gRNAs", "call", "detected_gRNAs")
    colData(sce) <- cbind(colData(sce), gRNA_calls[match(colnames(sce), gRNA_calls$cell), ..sel])
  }else{
    # add everything
    colData(sce) <- cbind(colData(sce), gRNA_calls[match(colnames(sce), gRNA_calls$cell),-1])
  }
  # remove NAs produced by cells with _no_ gRNA counts at all (missing from analysis above)
  if(sum(is.na(sce$n_detected_gRNAs))>0){
    colData(sce[,is.na(sce$n_detected_gRNAs)])[,'call'] <- "unassigned"
    colData(sce[,is.na(sce$n_detected_gRNAs)])[,'n_detected_gRNAs'] <- 0
  }
  
  ## return results
  if(plot == TRUE){
    return(list(sce = sce,
                calls = gRNA_calls,
                plot = p))
  }else{
    return(list(sce = sce,
                calls = gRNA_calls))
  }
}
```

We perform the test using the probability of selecting each gRNA based on their frequencies in the initial plasmid library. The vast majority of cells have a single gRNA detected, consistent with this being a low MOI experiment.

```{r call_gRNAs_cells, fig.width=10, fig.height=15}
## gRNA anntation - same library as previous pilots
gRNA_ann <- read.table(paste0(dir, "data/gRNA_library.tsv"), sep="\t", header = TRUE)

## gRNA representation in initial plasmid library
gRNA_prob <- read.table(paste0(dir, "data/gRNA_count_plasmid_library.tsv"), header = TRUE)
probs <- gRNA_prob$Count/sum(gRNA_prob$Count)
names(probs) <- gRNA_prob$gRNA_ID

## call gRNAs for each sample
sce.list <- sapply(unique(sce$Sample), function(x) sce[,sce$Sample == x])
gRNA_calls <- bplapply(sce.list, 
                       function(x) call_gRNAs_binom(sce = x, 
                                                    pval_threshold = 0.001,
                                                    prob = probs,
                                                    max_to_test = 5,
                                                    total_guides = nrow(altExp(sce)),
                                                    plot = TRUE, 
                                                    ann = gRNA_ann, 
                                                    columns = c("target", "expected_DE_gene", "class"),
                                                    min_count = 3, min_fraction = 0.05,
                                                    add_stats = FALSE),
                       BPPARAM = bpp)

## merge all samples
sce <- do.call(cbind, lapply(gRNA_calls, '[[', 1))
plots <- lapply(gRNA_calls, '[[', 3)
gRNA_calls <- lapply(gRNA_calls, '[[', 2)
gRNA_calls <- do.call(rbind, gRNA_calls)

## plot
ggarrange(plotlist = plots, ncol=5, nrow=9,
          labels = unique(sce$Sample), font.label = list(size=8),
          common.legend = TRUE, legend = "bottom")
```

Overall, around 63% of cells have a gRNA detected.

```{r prop_unique, fig.width=10, fig.height=5}
# proportion of cells with each call type
df <- round(prop.table(table(sce$Sample, sce$call),1)*100, 2)

par(mar=c(6,4,2,2))
barplot(t(df[,c(3,1,2)]), col=c("grey10", "indianred3", "grey"), 
        las=2, cex.names = 0.7,
        width = 0.8, space = 0.2, xlim=c(0,40))
legend("topright", legend = c("unassigned", "multiple", "unique"),
       col = c("grey", "indianred3", "grey10"), pch = 15, cex=0.7)
```

Moving forward we only use cells with a gRNA assignment. This leaves us with this many cells:

```{r remove_unass}
sce <- sce[,sce$call != "unassigned"]
c(n_cells=ncol(sce))
```


```{r save}
saveRDS(sce, paste0(dir, "results/02_sce_POC_Tcells.goodQual.calls.Rds"))
```


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

