---
title: "QC on demultiplexed scRNA-seq preadipocyte data (day-0)"
author: Anushka Gupta
date: February 16th 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)

# We demultiplexed scRNA-seq preadipocyte dataset to annotate barcodes as singlets, doublets or empty in demultiplexing.htos.Rmd notebook. Singlets were further annotated as white or brown in demultiplexing.htos.Rmd notebook. These barcodes can be found in rdata folder as .rds files
# The gene count matrix for singlet barcodes only 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/preadipocytes_post_demultiplexing.h5ad")

# Add latent space inferred using scVI to Seurat object.
latent <- read.csv("additional_files/latent_space/preadipocytes_post_demultiplexing.csv", row.names=1)
latent=as.matrix(latent)
rownames(latent)=rownames(dataset@meta.data)
dataset[["scVI"]] <- CreateDimReducObject(embeddings =latent,key = "scVI_")
dataset=FindNeighbors(dataset,reduction="scVI",graph.name="RNA_snn_scVI",dims=1:20)
dataset=RunUMAP(dataset,reduction="scVI",dims = 1:20)
UMAPPlot(dataset,group.by='type',label=TRUE,label.size=6)+theme_bw(base_size = 20)+theme(legend.position = "none")+scale_color_manual(values=c("brown", "grey"))

# Find Clusters 
dataset=FindClusters(dataset,graph.name = "RNA_snn_scVI",resolution = .4,algorithm = 1)
UMAPPlot(dataset,label=TRUE,cols=brewer.pal(n = 4,name = "Dark2"))+theme_bw(base_size = 20)
#### Perform QC on clusters to remove low quality clusters 
dataset=NormalizeData(dataset,verbose = FALSE)
dataset <- PercentageFeatureSet(dataset, pattern = "^MT-", col.name = "percent.mt")
dataset@meta.data %>% group_by(seurat_clusters) %>% summarise(mean(percent.mt))

# Cluster 3 is the only cluster which has > 5% MT-content.
# Check if DE also identifies MT genes
markers.3=FindMarkers(object = dataset,ident.1 = 3,assay = "RNA",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
markers.3=subset.data.frame(markers.3,subset = p_val_adj<.05)

# Only 21 marker genes with top DE genes for MT genes
# Therefore, remove cluster 3 from downstream analyses
dataset.subset=subset(dataset, cells =rownames(dataset@meta.data)[dataset@meta.data$seurat_clusters!=3])
UMAPPlot(dataset.subset,group.by='type',pt.size=1.5)+theme_void()+theme(legend.position = "none")+scale_color_manual(values=c("grey", "brown"))
UMAPPlot(dataset.subset,group.by='seurat_clusters',label=TRUE,label.size=6,cols=brewer.pal(n = 3,name = "Dark2"))+theme_bw(base_size = 20)+theme(legend.position = "none")
dataset.subset=SCTransform(dataset.subset,verbose = FALSE)
save(dataset.subset,file = "rdata/seurat.preadipocytes.post.demultiplexingAndQC.Rdata")
```
