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

We first performed individual QC on snRNA-seq white and brown preadipocyte dataset which could be found in resepctive R Notebooks in pre-processing folder. We then integrated snRNA-seq white and brown preadipocyte datasets together using scVI and performed additional QC in the notebook below. 

```{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(org.Hs.eg.db)
library(GO.db)
library(limma)

# The gene count matrix for snRNA-seq white and brown preadipocyte dataset was used as input to scVI package to infer the latent space. 
# 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
s_nuclei <- ReadH5AD(file ="additional_files/adata/nuclei_post_QC.h5ad")

# Add latent space inferred using scVI to Seurat object
latent <- read.csv("additional_files/latent_space/nuclei_post_QC.csv", row.names=1)
latent=as.matrix(latent)
rownames(latent)=rownames(s_nuclei@meta.data)
s_nuclei[["scVI"]] <- CreateDimReducObject(embeddings =latent,key = "scVI_")
s_nuclei=RunUMAP(s_nuclei,reduction="scVI",dims = 1:20)
UMAPPlot(s_nuclei,group.by="batch.indices",cols=c("grey","brown"))+theme_bw(base_size = 20)

# Clustering
s_nuclei=FindNeighbors(s_nuclei,reduction="scVI",graph.name="RNA_snn_scVI",dims=1:20)
s_nuclei=FindClusters(s_nuclei,graph.name = "RNA_snn_scVI",resolution = .24)
UMAPPlot(s_nuclei,label=TRUE)+theme_bw(base_size = 20)

# Clustering identified a common cluster 3 between white and brown preadipocytes.
# Performing DE to find marker genes for this common cluster # 3
s_nuclei=NormalizeData(s_nuclei)
x=FindMarkers(s_nuclei,ident.1 = 3,logfc.threshold = .5,only.pos = T)
# Performing GO on BC1 genes
entrez.ids <- mapIds(org.Hs.eg.db, keys=rownames(x), column="ENTREZID", keytype="SYMBOL")
entrez.ids.universe=mapIds(org.Hs.eg.db, keys=rownames(s_nuclei), 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)

# The common cluster # 3 was enriched for translational terms (suggesting high cellular background contamination), and stress-response terms 
# Therefore, remove cluster 3 from downstream analysis
s_nuclei=subset(s_nuclei,cells=WhichCells(s_nuclei,idents = c(0,1,2)))
s_nuclei=SCTransform(s_nuclei)
save(s_nuclei,file = "rdata/seurat.nuclei.preadipocytes.post.QC.Rdata")
```

