---
title: "Reproducing results discussed in Results-Section-2"
author: Anushka Gupta
date: March 24 2021
---

# This notebook has commands to recreate results discussed in Results-section-2: 2.2	snRNA-seq identifies the same preadipocyte populations as scRNA-seq and detects biologically relevant differential expression 
```{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(viridis)
library(dplyr)
library(limma)
library(org.Hs.eg.db)
library(VISION)
library(Matrix)
library(viridis)
size=20

# load seurat object containing white and brown preadipocyte snRNA-seq data after QC. See pre-processing folder for pre-processing notebooks 
load("rdata/seurat.nuclei.preadipocytes.post.QC.RData")
# Recreating panel (A) from Fig. 2
UMAPPlot(s_nuclei,group.by='batch.indices',pt.size=.5)+theme_void()+theme(legend.position = "none")+scale_color_manual(values=c("grey", "brown"))
UMAPPlot(s_nuclei,group.by="seurat_clusters",pt.size=.5,cols=brewer.pal(n = 3,name = "Dark2"),label=TRUE,label.size=5)+theme_void()+theme(legend.position = "none")

# Perform transcriptional signature analysis using Vision. This analysis takes time, hence you could directly load the vision object using the command in line # 49  
# Read in expression counts (Genes X Cells)
counts <- s_nuclei@assays$RNA@counts
# Scale counts within a sample
n.umi <- colSums(counts)
scaled_counts <- t(t(counts) / n.umi) * median(n.umi)
# Adding UMAP
projection <- s_nuclei@reductions$umap@cell.embeddings
# Read in meta data (Cells x Vars)
meta = s_nuclei@meta.data
# Read in signatures of white & brown and cluster-1 & cluster-2 preadipocytes generated using scRNA-seq dataset
# Signatures were derived by performing DE test between white and brown preadipocytes in scRNA-seq dataset, with enriched genes for respective preadipocyte type forming a signature
# Similarly, signatures were derived by performing DE test between cluster-1 and cluster-2 preadipocytes in scRNA-seq dataset, with enriched genes for respective cluster forming a signature
dir="additional_files/signatures/Nuclei_White_Brown/"
Signatures=list.files(path=dir,full.names = TRUE)
vis <- Vision(scaled_counts,signatures = Signatures,meta = meta)
vis <- addProjection(vis, "UMAP", projection)
vis <- analyze(vis)
save(vis,file = "rdata/vision.nuclei.preadipocytes.post.QC.RData")
# Or directly load Vision object
load(file = "rdata/vision.nuclei.preadipocytes.post.QC.RData")

# Plotting transcriptional signature score. Recerating panel (B) from Fig. 2
s_nuclei@meta.data$white_score=vis@SigScores[,"White "]
s_nuclei@meta.data$brown_score=vis@SigScores[,"Brown Brown"]
FeaturePlot(s_nuclei,"white_score",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)
FeaturePlot(s_nuclei,"brown_score",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)
s_nuclei@meta.data$one_cell_score=vis@SigScores[,"one_cells"]
FeaturePlot(s_nuclei,"one_cell_score",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)
s_nuclei@meta.data$two_cell_score=vis@SigScores[,"two_cells"]
FeaturePlot(s_nuclei,"two_cell_score",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)

# Expression of genes enriched in white-preadipocytes or brown-preadipocytes in scRNA-seq dataset
Idents(s_nuclei)="batch.indices"
# Enrichment of Tcf21, Pax3, and PDGFRa
FindMarkers(s_nuclei,ident.1 = 0,assay = "SCT",slot = "data",features = c("TCF21","PAX3","PDGFRA"),logfc.threshold = 0,only.pos = FALSE,min.pct = 0)
#Enrichment of ID1
FindMarkers(s_nuclei,ident.1 = 0,assay = "SCT",slot = "data",features = "ID1",logfc.threshold = 0,only.pos = FALSE,min.pct = 0)

#Enrichment of ANKRD1 and CTGF
FindMarkers(s_nuclei,ident.1 = 0,assay = "SCT",slot = "data",features = c("ANKRD1","CTGF"),logfc.threshold = 0,only.pos = FALSE,min.pct = 0)
#Enrichment of ACTCA2, TAGLN, MYH11, MYL9, CNN1
FindMarkers(s_nuclei,ident.1 = 0,assay = "SCT",slot = "data",features = c("ACTA2","TAGLN","MYH11","MYL9","CNN1"),logfc.threshold = 0,only.pos = FALSE,min.pct = 0)

# Expression of genes classified as markers of white and brown preadipocytes in scRNA-seq dataset
m_w_cells=read.table("additional_files/markers/white.cells.txt", quote="\"", comment.char="")
m_w_cells=m_w_cells$V1
m_b_cells=read.table("additional_files/markers/brown.cells.txt", quote="\"", comment.char="")
m_b_cells=m_b_cells$V1
FindMarkers(object = s_nuclei,features = m_w_cells,ident.1 = 0,assay = "SCT",slot = "data",verbose = FALSE,only.pos = FALSE,logfc.threshold = 0,min.pct = 0)
FindMarkers(object = s_nuclei,features = m_b_cells,ident.1 = 0,assay = "SCT",slot = "data",verbose = FALSE,only.pos = FALSE,logfc.threshold = 0,min.pct = 0)

# Are top genes that are DE in scRNA-seq also DE in snRNA-seq
load("rdata/seurat.preadipocytes.post.demultiplexingAndQC.RData")
n=50
Idents(dataset.subset)="cell_type"
deg_w_cells=FindMarkers(object = dataset.subset,ident.1 = "WAT",ident.2 = "BAT",assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
deg_w_cells=subset.data.frame(deg_w_cells,subset = p_val_adj<.05)
top_deg_w_cells=rownames(deg_w_cells)[order(deg_w_cells$avg_logFC,decreasing = TRUE)[1:n]]
test=FindMarkers(object = s_nuclei,features = top_deg_w_cells,ident.1 = 0,assay = "SCT",slot = "data",verbose = FALSE,only.pos = TRUE,logfc.threshold = 0,min.pct = 0)
test=subset.data.frame(test,subset = p_val_adj<.05)
nrow(test)

deg_b_cells=FindMarkers(object = dataset.subset,ident.1 = "BAT",ident.2 = "WAT",assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
deg_b_cells=subset.data.frame(deg_b_cells,subset = p_val_adj<.05)
top_deg_b_cells=rownames(deg_b_cells)[order(deg_b_cells$avg_logFC,decreasing = TRUE)[1:n]]
test=FindMarkers(object = s_nuclei,features = top_deg_b_cells,ident.1 = 1,assay = "SCT",slot = "data",verbose = FALSE,only.pos = TRUE,logfc.threshold = 0,min.pct = 0)
test=subset.data.frame(test,subset = p_val_adj<.05)
nrow(test)
```

```{r}
# Examining transcriptional differences between cluster-1 and cluster-2 preadipocytes
Idents(s_nuclei)="seurat_clusters"
deg_b1_nuclei=FindMarkers(object = s_nuclei,ident.1 = 1,ident.2 = 2,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
deg_b1_nuclei=subset.data.frame(deg_b1_nuclei,subset = p_val_adj<.05)
deg_b2_nuclei=FindMarkers(object = s_nuclei,ident.1 = 2,ident.2 = 1,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
deg_b2_nuclei=subset.data.frame(deg_b2_nuclei,subset = p_val_adj<.05)

write.table(x = rownames(deg_b1_nuclei),file = "additional_files/de_genes/1.vs.2.nuclei.txt",sep = "\n",row.names = FALSE,col.names = FALSE,quote = FALSE)
write.table(x = rownames(deg_b2_nuclei),file = "additional_files/de_genes/2.vs.1.nuclei.txt",sep = "\n",row.names = FALSE,col.names = FALSE,quote = FALSE)

# Performing gene ontology analysis and identifying common GO terms between scRNA-seq and snRNA-seq datasets
GO.1.vs.2.nuclei <- read.delim("additional_files/gene ontology/results/1.vs.2.nuclei.txt", comment.char="#")
rownames(GO.1.vs.2.nuclei)=GO.1.vs.2.nuclei$GO.biological.process.complete
GO.2.vs.1.nuclei <- read.delim("additional_files/gene ontology/results/2.vs.1.nuclei.txt", comment.char="#")
rownames(GO.2.vs.1.nuclei)=GO.2.vs.1.nuclei$GO.biological.process.complete

GO.1.vs.2.cells <- read.delim("additional_files/gene ontology/results/1.vs.2.cells.txt", comment.char="#")
GO.1.vs.2.cells=GO.1.vs.2.cells[order(GO.1.vs.2.cells$upload_1..FDR.,decreasing = F),]
GO.2.vs.1.cells <- read.delim("additional_files/gene ontology/results/2.vs.1.cells.txt", comment.char="#")
GO.2.vs.1.cells=GO.2.vs.1.cells[order(GO.2.vs.1.cells$upload_1..FDR.,decreasing = F),]

terms.1.vs.2=intersect(GO.1.vs.2.cells$GO.biological.process.complete[1:10],GO.1.vs.2.nuclei$GO.biological.process.complete)
terms.2.vs.1=intersect(GO.2.vs.1.cells$GO.biological.process.complete[1:10],GO.2.vs.1.nuclei$GO.biological.process.complete)

# Expression of BMP4 and TGFB2
FindMarkers(object = s_nuclei,features = c("BMP4","TGFB2"),ident.1 = 1,ident.2 = 2,assay = "SCT",slot = "data",logfc.threshold = 0,min.pct = 0,verbose = FALSE,only.pos = TRUE)
#### Expression of FST, MEST, NPPB and SEMA5A
FindMarkers(object = s_nuclei,features = c("FST","NPPB","MEST","SEMA5A"),ident.1 = 1,ident.2 = 2,assay = "SCT",slot = "data",logfc.threshold = 0,min.pct = 0,verbose = FALSE,only.pos = TRUE)

# Expression of genes classified as markers of cluster 1 preadipocytes in scRNA-seq dataset
m_b1_cells=read.table("additional_files/markers/1.cells.txt", quote="\"", comment.char="")
m_b1_cells=m_b1_cells$V1
FindMarkers(object = s_nuclei,features = m_b1_cells,ident.1 = 1,ident.2 = 2,assay = "SCT",slot = "data",verbose = FALSE,only.pos = FALSE,logfc.threshold = 0,min.pct = 0)

# Expression of top 50 DE genes in cluster 1 in scRNA-seq dataset
n=50
Idents(dataset.subset)="seurat_clusters"
deg_b1_cells=FindMarkers(object = dataset.subset,ident.1 = 1,ident.2 = 2,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
deg_b1_cells=subset.data.frame(deg_b1_cells,subset = p_val_adj<.05)
top_deg_b1_cells=rownames(deg_b1_cells)[order(deg_b1_cells$avg_logFC,decreasing = TRUE)[1:n]]
test=FindMarkers(object = s_nuclei,features = top_deg_b1_cells,ident.1 = 1,ident.2 = 2,assay = "SCT",slot = "data",verbose = FALSE,only.pos = TRUE,logfc.threshold = 0,min.pct = 0)
test=subset.data.frame(test,subset = p_val_adj<.05)
nrow(test)

# Expression of S100A4, MFAP5, ECM1, COL6A2, and ACAN
FindMarkers(object = s_nuclei,features = c("S100A4","MFAP5","ECM1","COL6A2","ACAN"),ident.1 = 1,ident.2 = 2,assay = "SCT",slot = "data",logfc.threshold = 0,min.pct = 0,verbose = FALSE,only.pos = F)

# Expression of genes classified as markers of cluster 2 preadipocytes in scRNA-seq dataset
m_b2_cells=read.table("additional_files/markers/2.cells.txt", quote="\"", comment.char="")
m_b2_cells=m_b2_cells$V1
FindMarkers(object = s_nuclei,features = m_b2_cells,ident.1 = 1,ident.2 = 2,assay = "SCT",slot = "data",verbose = FALSE,only.pos = FALSE,logfc.threshold = 0,min.pct = 0)

# Expression of top 50 DE genes in cluster 2 in scRNA-seq dataset
deg_b2_cells=FindMarkers(object = dataset.subset,ident.1 = 2,ident.2 = 1,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
deg_b2_cells=subset.data.frame(deg_b2_cells,subset = p_val_adj<.05)
top_deg_b2_cells=rownames(deg_b2_cells)[order(deg_b2_cells$avg_logFC,decreasing = TRUE)[1:n]]
test=FindMarkers(object = s_nuclei,features = top_deg_b2_cells,ident.1 = 2,ident.2 = 1,assay = "SCT",slot = "data",verbose = FALSE,only.pos = TRUE,logfc.threshold = 0,min.pct = 0)
test=subset.data.frame(test,subset = p_val_adj<.05)
nrow(test)
```

