```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.path = '.images/2/')
```

```{r, results='hide', message=FALSE, warning=FALSE}
library(tidyverse)
library(Seurat)

library(BUSpaRse)
library(SeuratWrappers)
library(DropletUtils)
```

We get two matrices here, for both spliced and unspliced counts. Note that these are (luckily) loaded in as sparse matrices.

```{r}
sample_folder <- function(x) {paste0("data/adipose_differentiation-merick/SRR8906269/", x)}
velocity_output <- read_velocity_output(
  spliced_dir = sample_folder("introns"),
  spliced_name = "s",
  unspliced_dir = sample_folder("cdna"),
  unspliced_name = "u"
)
spliced <- velocity_output$spliced
unspliced <- velocity_output$unspliced
```

```{r}
dim(spliced)
```

```{r}
dim(unspliced)
```

The spliced matrix contains a lot more cells than the unspliced matrix. A first indication that we have much more spliced UMIs than unspliced ones.

# Filtering empty droplets

The kallisto/BUStools pipeline does not do any filtering on the features or cells (except for zero counts), so it makes sense that most droplets contain only a very tiny number of UMIs.

```{r}
qplot(Matrix::colSums(spliced) + 1) +
  scale_x_log10("# UMIs", labels = function(x) as.integer(x) - 1) +
  theme_classic()
```

We will filter these empty droplets here by finding the inflection point within the total UMI count distribution.

```{r}
bc_rank_spliced <- DropletUtils::barcodeRanks(spliced)
bc_rank_unspliced <- DropletUtils::barcodeRanks(unspliced)
```

This function gives us the total UMI counts at which the curves makes an inflection or knee

```{r}
metadata(bc_rank_spliced)$inflection
metadata(bc_rank_spliced)$knee
```


```{r}
tibble(rank = bc_rank_spliced$rank, total = bc_rank_spliced$total, matrix = "spliced") %>%
  bind_rows(tibble(rank = bc_rank_unspliced$rank, total = bc_rank_unspliced$total, matrix = "unspliced")) %>%
  distinct() %>%
  ggplot(aes(total, rank, color = matrix)) +
  geom_line() +
  scale_color_manual(values = c(spliced = "#0074D9", unspliced = "#FF4136")) +
  geom_vline(xintercept = metadata(bc_rank_spliced)$inflection, color = "#FF4136", linetype = 3) +
  geom_vline(xintercept = metadata(bc_rank_unspliced)$inflection, color = "#0074D9", linetype = 3) +
  scale_x_log10() +
  scale_y_log10() +
  labs(y = "Rank", x = "Total UMI counts") +
  theme_classic()
```

(Note that this rank vs total count plot is transposed compared to the plot shown typically shown by e.g. CellRanger)

In this sample, it is clear that a sizeable group of cells has large spliced UMI counts (1000 - 10000). Note that although the total UMI counts are different between spliced and unspliced counts, the y location of the inflection point does not change that much. Moreover, about same empty droplets are retained in both analyses:

```{r}
bcs_filtered_spliced <- colnames(spliced)[colSums(spliced) > metadata(bc_rank_spliced)$inflection]
bcs_filtered_unspliced <- colnames(unspliced)[colSums(unspliced) > metadata(bc_rank_spliced)$inflection]

length(intersect(bcs_filtered_unspliced, bcs_filtered_spliced)) / length(union(bcs_filtered_unspliced, bcs_filtered_spliced))
```

We'll filter using the spliced inflection point
```{r}
spliced_filtered <- spliced[, bcs_filtered_spliced]
unspliced_filtered <- unspliced[, bcs_filtered_spliced]

features_filtered <- rownames(spliced_filtered)[Matrix::rowSums(spliced_filtered) > 0]

spliced_filtered <- spliced_filtered[features_filtered,]
unspliced_filtered <- unspliced_filtered[features_filtered,]
```

```{r}
dim(spliced_filtered)
```

# RNA velocity QC

We'll have a quick look at some statistics about the spliced vs unspliced reads in the sample, to confirm whether the separate quantification was succesful.

## Proportion of unspliced reads

```{r}
overall_ratio <- sum(unspliced_filtered@x) / (sum(unspliced_filtered@x) + sum(spliced_filtered@x))
overall_ratio
```

The proportion of unspliced reads is typically around 20%, irregardless of the technology. This was already described in [@mannoRNAVelocitySingle2018]. Although slight deviations between 10% and 50% are possible (especially outside of mammals), extreme values might indicate that something went wrong in the detection of unspliced reads.

## Gene-wise proportion

```{r}
calculate_ratio <- function(x, y) {
  ratios <- x / y
  ratios[x == 0] <- 0
  ratios[y == 0] <- Inf
  ratios[x == 0 & y == 0] <- NA
  ratios
}
ratios <- calculate_ratio(Matrix::rowSums(unspliced_filtered), Matrix::rowSums(spliced_filtered))
ratios <- tibble(ratio = ratios, gene = rownames(unspliced_filtered))
```

```{r}
ratios %>% 
  filter(is.finite(ratio) & ratio > 0) %>% 
  ggplot(aes(x = ratio)) +
    geom_density() +
    scale_x_log10("Ratio per gene", breaks = c(1e-4, 1e-3, 1e-2, 1e-1, 1, 1e1, 1e2, 1e3), labels = identity) +
    geom_vline(xintercept = 1, linetype = "dashed") +
    geom_vline(xintercept = overall_ratio, color = "red") +
    theme_classic() +
    scale_y_continuous("Gene density", expand = c(0, 0)) +
    annotate(geom = "text", label = "More unspliced", y = Inf, x = Inf, hjust = 1, vjust = 1) +
    annotate(geom = "text", label = "More spliced", y = Inf, x = 0, hjust = 0, vjust = 1)
```

The proportions spliced vs unspliced can be very variable between genes. This makes sense, given the (probable) variability of splicing efficiency and the presence and size of introns. Some of this variability might also be caused by incorrect annotations.

```{r}
ratios %>% 
  mutate(ratio = ifelse(ratio > 0 & ratio < Inf, 1, ratio)) %>% 
  pull(ratio) %>% 
  table(useNA = "ifany")
```

0 = no unspliced  
1 = both spliced and unspliced  
Inf = no spliced  
NA = no spliced or unspliced  

Ideally, a majority of genes should have both spliced and unspliced UMIs. Of course, a sizeable set of genes do not have any unspliced UMIs.

## Cell-wise proportion

```{r}
common_cells <- intersect(colnames(unspliced_filtered), colnames(spliced_filtered))
ratios <- calculate_ratio(Matrix::colSums(unspliced_filtered[, common_cells]), Matrix::colSums(spliced_filtered[, common_cells]))
ratios <- tibble(ratio = ratios, cell = common_cells)
```

```{r}
ratios %>% 
  filter(is.finite(ratio) & ratio > 0) %>% 
  ggplot(aes(x = ratio)) +
    geom_density() +
    scale_x_log10("Ratio per cell", breaks = c(1e-4, 1e-3, 1e-2, 1e-1, 1, 1e1, 1e2, 1e3), labels = identity) +
    geom_vline(xintercept = 1, linetype = "dashed") +
    geom_vline(xintercept = overall_ratio, color = "red") +
    theme_classic() +
    scale_y_continuous("Cell density", expand = c(0, 0)) +
    annotate(geom = "text", label = "More unspliced", y = Inf, x = Inf, hjust = 1, vjust = 1) +
    annotate(geom = "text", label = "More spliced", y = Inf, x = 0, hjust = 0, vjust = 1)
```

The proportions spliced vs unspliced per cell is typically more centered around the "average" proportion. If you see some weird distributions here, this might indicate that the filtering of cells was not stringent enough, or that there are batch effects.

```{r}
ratios %>% 
  mutate(ratio = ifelse(ratio > 0 & ratio < Inf, 1, ratio)) %>% 
  pull(ratio) %>% 
  table(useNA = "ifany")
```

0 = no unspliced  
1 = both spliced and unspliced  
Inf = no spliced  
NA = no spliced or unspliced 

Ideally, most of the cells should have both spliced and unspliced UMIs.

## Normalization

Let's do a basic log2 and size factors normalization. SCTransform [@hafemeisterNormalizationVarianceStabilization2019] can also be used here.

```{r}
seu <- CreateSeuratObject(spliced_filtered, assay = "spliced", min.cells = 3, min.features = 200) %>%
  NormalizeData() %>%
  FindVariableFeatures() %>% 
  ScaleData()

seu[["unspliced"]] <- CreateAssayObject(unspliced_filtered[rownames(seu@assays$spliced@counts), colnames(seu@assays$spliced@counts)]) %>%
  NormalizeData()
```

We only run ScaleData here on the spliced counts, as this is required to run PCA.

Let's do some basic clustering and dimensionality reduction.

```{r}
DefaultAssay(seu) <- "spliced"
seu <- RunPCA(seu, verbose = FALSE, npcs = 70)
```

```{r}
set.seed(1)

seu <- FindNeighbors(seu, verbose = FALSE) %>% 
  FindClusters(resolution = 0.1, verbose = FALSE)
```

```{r}
DimPlot(seu, reduction = "pca",pt.size = 0.5, label = TRUE, repel = TRUE)
```


```{r}
set.seed(1)

seu <- RunUMAP(seu, dims = 1:50, umap.method = "uwot")
```

```{r}
DimPlot(seu, reduction = "umap",pt.size = 0.5, label = TRUE, repel = TRUE)
```

Do you see the trajectory appearing? Most TI methods assume that only one trajectory is present in the dataset, with all cells being connected to this trajectory. We therefore filter those cells out, although you may want to try a different method, such as PAGA, if you don't want to bias the analysis.

```{r}
cells_trajectory <- seu@meta.data %>%
  mutate(cell_id = Cells(seu)) %>%
  filter(seurat_clusters %in% c(0, 1, 2)) %>%
  pull(cell_id)

seu_trajectory <- seu[, cells_trajectory]
```

We'll have a look at RNA velocity and trajectories in the next notebooks.

```{r}
write_rds(seu_trajectory, "data/adipose_differentiation-merick/seu_trajectory.rds")
```

## References
