prep_filtration.Rmd

---
title: "Pre-processing of new datasets"
output: html_notebook
date: "`r format(Sys.time(), '%d %B, %Y')`"
---

```{r, message=FALSE, warning=FALSE}
# library(Epilepsy19)
library(ggplot2)
library(magrittr)
library(Matrix)
library(pbapply)
library(conos)

devtools::load_all()

theme_set(theme_bw())
```

Load data:

```{r}
cms <- list.files(DataPath(""), pattern="*.h5") %>% setNames(., gsub(".h5", "", .)) %>% 
  pblapply(function(p) DataPath(p) %>% Seurat::Read10X_h5() %>% .[, colSums(.) > 10])
```

Filter small cells:

```{r, fig.width=12, fig.height=13}
thresholds <- c(Biopsy=3.05, GTS213=2.5, GTS217=1.0, GTS217_2=3.2, GTS219=2.5, GTS233=3.0, 
                HB26=3.4, HB51=3.2, HB52=2.85, HB53=3.1, HB56=2.75, HB65=2.85, NeuN=3.0, 
                CTR215=2.7, CTR240=3.0)

lapply(names(cms), function(n) 
  dropestr::PlotCellsNumberHist(colSums(cms[[n]]), estimate.cells.number=T, show.legend=F) + 
    geom_vline(aes(xintercept=thresholds[[n]]))) %>% 
  cowplot::plot_grid(plotlist=., ncol=3, labels=names(cms))
```

```{r}
cms <- names(cms) %>% setNames(., .) %>% 
  pblapply(function(n) cms[[n]] %>% .[, colSums(.) >= 10**thresholds[[n]]])
```

Append matrices from the previous runs:

```{r}
cm_paths_old <- list.files(DataPath(""), pattern=".*_p.*") %>% 
  setNames(., .) %>% sapply(DataPath)

cms_old <- pbapply::pblapply(cm_paths_old, pagoda2::read.10x.matrices, verbose=F, cl=10)

cms %<>% c(cms_old)
sapply(cms, ncol)
```

```{r}
sapply(cms, ncol) %>% .[!names(.) %in% c("GTS217", "NeuN")] %>% sum()
```

Mitochondrial fraction:

```{r}
mit_frac_per_dataset <- pblapply(cms, function(cm) 
  colSums(cm[grep("MT-", rownames(cm)), ]) / colSums(cm))
```

```{r, fig.width=12, fig.height=12}
lapply(mit_frac_per_dataset, function(fr) 
  qplot(fr[fr < 0.5], xlab="Mit. fraction", ylab="#Cells", xlim=c(-0.01, 0.5), bins=30)) %>% 
  cowplot::plot_grid(plotlist=., ncol=4, labels=names(cms))
```

```{r}
cms <- pbmapply(function(cm, frac) cm[, frac < 0.08], cms, mit_frac_per_dataset) %>% 
  setNames(names(cms))

cms$GTS213 <- NULL
cms$GTS217 <- NULL

sapply(cms, ncol)
```

Doublet detection:

```{r}
doublet_info <- pblapply(cms, GetScrubletScores, "~/mh/local/anaconda3/bin/python3.7", 
                         min.molecules.per.gene=50, cl=25)
```

```{r}
is_doublet <- lapply(doublet_info, `[[`, "score") %>% sapply(`>`, 0.25)
sapply(is_doublet, mean) %>% round(3)
```

## Filter matrices

```{r}
cms_filt <- cms
for (n in names(cms_filt)) {
  cms_filt[[n]] %<>% .[!grepl("^MT-", rownames(.)), !is_doublet[[n]]]
}

sapply(cms_filt, dim)
sapply(cms_filt, ncol) %>% .[names(.) != "NeuN"] %>% sum()
```

## Save data

```{r, eval=FALSE}
readr::write_rds(cms_filt, CachePath("count_matrices.rds"))
```

## Plot the doublet UMAP

**Note**: this part requires annotation, so it was run separately after performing annotation of the filtered data

```{r}
cms_renamed <- names(cms) %>% setNames(., .) %>% 
  lapply(function(n) set_colnames(cms[[n]], paste0(n, "_", colnames(cms[[n]]))))
```

```{r}
p2s <- lapply(cms_renamed, function(cm) cm %>% set_rownames(make.unique(rownames(.))) %>% 
                GetPagoda(clustering.type=NULL, embeding.type=NULL, build.graph=FALSE))
```

```{r}
con <- conos::Conos$new(p2s, n.cores=30)
con$buildGraph(k=50, k.self=5, k.self.weight=0.1)
con$embedGraph(method="UMAP", n.cores=45, min.prob.lower=1e-5, n.neighbors=30, 
               n.epochs=1000, spread=5, min.dist=1.0, verbose=T)
```

```{r, message=FALSE}
annotation <- readr::read_csv(MetadataPath("annotation.csv")) %$%
  setNames(l4, cell)

is_doublet_plt <- names(is_doublet) %>% lapply(function(n) 
  is_doublet[[n]] %>% setNames(paste0(n, "_", names(.)))) %>% 
  Reduce(c, .)

doublet_score_plt <- names(doublet_info) %>% lapply(function(n) 
  doublet_info[[n]]$score %>% setNames(paste0(n, "_", names(.)))) %>% 
  Reduce(c, .)
```

```{r, fig.width=7, fig.height=7}
gg_annot <- con$plotGraph(groups=annotation, shuffle.colors=T, raster=T)
color_scale <- scale_color_manual(values=c("#00913f", "#910000"))

gg_annot
ggsave(OutputPath("fig_filtration", "annotation.pdf"))

is_excluded <-setNames(!(names(con$getDatasetPerCell()) %in% names(annotation)), names(con$getDatasetPerCell()))
con$plotGraph(groups=is_excluded, mark.groups=F, show.legend=T, legend.pos=c(1, 1), size=0.1, alpha=0.2, raster=T) +
  gg_annot$layers[[2]] +
  scale_size_continuous(range=c(3, 5), trans='identity', guide='none') +
  guides(color=guide_legend(override.aes=list(size=3, alpha=1), title="Is excluded")) +
  color_scale

ggsave(OutputPath("fig_filtration", "excluded.pdf"))

con$plotGraph(groups=is_doublet_plt, mark.groups=F, show.legend=T, legend.pos=c(1, 1), size=0.1, alpha=0.2, raster=T) +
  gg_annot$layers[[2]] +
  scale_size_continuous(range=c(3, 5), trans='identity', guide='none') +
  guides(color=guide_legend(override.aes=list(size=3, alpha=1), title="Is doublet")) +
  color_scale

ggsave(OutputPath("fig_filtration", "doublets.pdf"))

con$plotGraph(colors=doublet_score_plt, show.legend=T, legend.title="Doublet score", legend.pos=c(1, 1), raster=T) +
  gg_annot$layers[[2]] +
  scale_size_continuous(range=c(3, 5), trans='identity', guide='none')

ggsave(OutputPath("fig_filtration", "doublet_scores.pdf"))
```

```{r, message=FALSE}
samp_aliases <- readr::read_csv(MetadataPath("sample_info.csv")) %$% setNames(Alias, Name) %>% c(NeuN="NeuNNeg")
doublet_frac <- sapply(is_doublet, mean)
excluded_frac <- split(is_excluded, con$getDatasetPerCell()[names(is_excluded)]) %>% sapply(mean)
excluded_frac["NeuN"] <- NA

doublet_report_df <- tibble(Sample=samp_aliases[names(doublet_frac)], `Doublet Fraction`=round(doublet_frac, 3), 
                            `Excluded Fraction`=round(excluded_frac[names(doublet_frac)], 3), `Doublet Threshold`=0.25)

readr::write_csv(doublet_report_df, OutputPath("fig_filtration", "doublet_report.csv"))

doublet_report_df
```