---
title: "QC on snRNA-seq white adipocyte data (day-20)"
author: Anushka Gupta
date: March 29 2021
---

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

```{r}
# Load libraries
library(Seurat)
library(Matrix)
library(ggplot2)
library(RColorBrewer)
library(dplyr)
library(limma)
library(org.Hs.eg.db)
library(GO.db)
size=20

# The gene count matrix for snRNA-seq white adipocyte swhite_20 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
swhite_20 <- ReadH5AD(file ="additional_files/adata/white_20.h5ad")
swhite_20@meta.data$nCount_RNA=colSums(swhite_20@assays$RNA@counts)
swhite_20@meta.data$nFeatures_RNA=colSums(swhite_20@assays$RNA@counts>0)

# Add latent space inferred using scVI to Seurat object.
lwhite_20 <- read.csv("additional_files/latent_space/white_20.csv", row.names=1)
lwhite_20=as.matrix(lwhite_20)
rownames(lwhite_20)=rownames(swhite_20@meta.data)
swhite_20[["scVI"]] <- CreateDimReducObject(embeddings =lwhite_20,key = "scVI_")
swhite_20=RunUMAP(swhite_20,reduction="scVI",dims = 1:20)
UMAPPlot(swhite_20)

# Normalize Data
swhite_20=NormalizeData(swhite_20)
# Clustering 
swhite_20=FindNeighbors(swhite_20,reduction="scVI",graph.name="RNA_snn_scVI",dims=1:20)
swhite_20=FindClusters(swhite_20,resolution = .45,graph.name = "RNA_snn_scVI", verbose = FALSE)
UMAPPlot(swhite_20,label=TRUE)+theme_bw(base_size = 20)+theme(legend.position = "none")

# 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/white_20.csv", quote="\"", comment.char="")
cells.cellbender=cells.cellbender$V1
cells=colnames(swhite_20@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 3 predominantly consisted of empty barcodes
UMAPPlot(swhite_20,cells.highlight = cells.cellranger,label = TRUE,label.size = 6) + theme_bw(base_size = 20)+theme(legend.position = "none")

# Cluster 5 had high MT content
# Visualize percent mt content
swhite_20[["percent.mt"]] <- PercentageFeatureSet(swhite_20, pattern = "^MT-")
FeaturePlot(swhite_20,"percent.mt")

# Cluster 4 was enriched for genes suggesting cellular debris based on GO analysis. 
# Enrichemnt for terms associated with trnaslation, protein localization to ER was observed
# Cluster 4 also had ~ 11% mitochondrial content
i_2=FindMarkers(object = swhite_20,ident.1 = 4,ident.2=c(0,1),assay = "RNA",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
i_2=subset.data.frame(i_2,subset = p_val_adj<.05)
# GO analysis on cluster 4
# Showing results from go.db
entrez.ids <- mapIds(org.Hs.eg.db, keys=rownames(i_2), column="ENTREZID", keytype="SYMBOL")
entrez.ids.universe=mapIds(org.Hs.eg.db, keys=rownames(swhite_20), column="ENTREZID", keytype="SYMBOL")
go.out <- goana(unique(entrez.ids), species="Hs",universe=unique(entrez.ids.universe))
go.out <- go.out[order(go.out$P.DE),]
go.useful <- go.out[go.out$Ont=="BP" & go.out$N <= 600,]
head(go.useful, 20)

# Doublets were detected using doublet_detection
# Cluster 2 mostly consisted of doublets
# Visualize doublets from doublet detection. Output can be found in additional_files
dwhite_20 <- read.table("additional_files/doublet_white_20.txt", quote="\"", comment.char="")
swhite_20@meta.data$doublet=dwhite_20$V1  
dwhite_20=rownames(swhite_20@meta.data)[swhite_20@meta.data$doublet==1]
dwhite_20=dwhite_20[!is.na(dwhite_20)]
dwhite_20=union(dwhite_20,rownames(swhite_20@meta.data)[swhite_20@meta.data$RNA_snn_scVI_res.0.45==2])
UMAPPlot(swhite_20,cells.highlight = dwhite_20,label = TRUE,label.size = 6) + theme_bw(base_size = 20)+theme(legend.position = "none")
```
```{r}
# Filter for cells from cluster 2,3,4 and 5
bwhite_20=rownames(swhite_20@meta.data)[swhite_20@meta.data$RNA_snn_scVI_res.0.45!=2 & swhite_20@meta.data$RNA_snn_scVI_res.0.45!=3 & swhite_20@meta.data$RNA_snn_scVI_res.0.45!=4 & swhite_20@meta.data$RNA_snn_scVI_res.0.45!=5]

# Remove remaining doublets and low-quality barcodes
bwhite_20=intersect(bwhite_20,cells.common)
bwhite_20=setdiff(bwhite_20,dwhite_20)
UMAPPlot(swhite_20,cells.highlight=bwhite_20)
swhite_20=subset(swhite_20,cells=bwhite_20)
save(swhite_20,file = "rdata/seurat.white20.postQC.RData")
```
