---
title: "QC on snRNA-seq brown preadipocyte data (day-0)"
author: Anushka Gupta
date: March 24 2021
---

```{r setup}
knitr::opts_knit$set(root.dir = '~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/')
```

```{r}
# Load libraries
library(Seurat)
library(ggplot2)
library(RColorBrewer)
library(dplyr)
library(readr)

# The gene count matrix for snRNA-seq brown preadipocyte dataset was used as input to scVI package to infer the latent space. While using scVI, singlet barcodes with < 200 genes were removed. 
# We used this inferred latent space as input in this notebook, along with the anndata object used for analysis in scVI
# Find latent space in additional_files --> latent_space 
# Find anndata object in additional_files --> adata
dataset <- ReadH5AD(file ="additional_files/adata/brown_00.h5ad")

# Add latent space inferred using scVI to Seurat object.
latent <- read.csv("additional_files/latent_space/brown_00.csv", row.names=1)
latent=as.matrix(latent)
rownames(latent)=rownames(dataset@meta.data)
dataset[["scVI"]] <- CreateDimReducObject(embeddings =latent,key = "scVI_")
dataset=RunUMAP(dataset,reduction="scVI",dims = 1:20)
UMAPPlot(dataset)+theme_bw(base_size = 20)+theme_bw(base_size = 20)+theme(legend.position = "none")

# Clustering 
dataset=FindNeighbors(dataset,reduction="scVI",graph.name="RNA_snn_scVI",dims=1:20)
dataset=FindClusters(dataset,graph.name = "RNA_snn_scVI",resolution = .06,algorithm = 1)
UMAPPlot(dataset,label=TRUE)+theme_bw(base_size = 20)+theme(legend.position = "none")

# Cluster 2 has high MT content
dataset[["percent.mt"]] <- PercentageFeatureSet(dataset, pattern = "^MT-")
ggplot(dataset@meta.data, aes(x = as.factor(RNA_snn_scVI_res.0.06), y=percent.mt)) +geom_boxplot() + theme_bw(base_size = 20)

# We removed empty droplets using CellBender package. To remove barcodes that could be empty, we only worked with barcodes that were called as "cell-containing" by both CellRanger and CellBender.
# Reading barcodes called as "cells" by CellBender below 
cells.cellbender=read.table("additional_files/cellbender/brown_00.csv", quote="\"", comment.char="")
cells.cellbender=cells.cellbender$V1
cells=colnames(dataset@assays$RNA@data)
cells.common=intersect(cells.cellbender,cells)
cells.cellranger=setdiff(cells,cells.common)
# Visualizing barcodes called as "cells" only by cellranger, which are potentially empty droplets (as they were called empty by CellBender)
# Cluster 1 predominatly consisted of empty barcodes
UMAPPlot(dataset,cells.highlight = cells.cellranger,label = TRUE,label.size = 6) + theme_bw(base_size = 20)+theme(legend.position = "none")

# Therefore, removing all cluster 2 cells and cluster 1 cells and any barcodes in cluster 0 that were not called "cells" by CellBender
cells.cluster.1=rownames(dataset@meta.data)[dataset@meta.data$seurat_clusters==1]
cells.cluster.2=rownames(dataset@meta.data)[dataset@meta.data$seurat_clusters==2]
barcodes.high.quality=setdiff(cells.common,cells.cluster.1)
barcodes.high.quality=setdiff(barcodes.high.quality,cells.cluster.2)
s_cells=subset(dataset,cells = barcodes.high.quality)

# We used Scrublet to remove doublets
scrublet_output_table <- read_csv("additional_files/scrublet_brown_00.csv")
s_cells@meta.data$isAdoublet=NA
s_cells@meta.data[scrublet_output_table$Barcodes,"isAdoublet"]=scrublet_output_table$predicted_doublet
doublets=rownames(s_cells@meta.data)[s_cells@meta.data$isAdoublet=="TRUE"]
UMAPPlot(s_cells,cells.highlight = doublets,label = TRUE,label.size = 6)  +theme_bw(base_size = 20)+theme(legend.position = "none")

#### Removing doublets from barcodes.high.quality
singlets.high.quality=setdiff(barcodes.high.quality,doublets)
UMAPPlot(s_cells,cells=singlets.high.quality)
```