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

# This notebook has commands to recreate results discussed in Results-section-5: snRNA-seq detects relevant transcriptional regulation during adipogenesis in white preadipocytes 
```{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(viridis)
library(tibble)
library(clusterCrit)
library(dplyr)
library(ComplexHeatmap)
library("circlize")
library(DropletUtils)
library(VISION)

size=10

# Read scRNA-seq and snRNA-seq adipocyte data
# See pre_processing folder for QC notebooks
load(file = "rdata/seurat.cellenone.postQC.RData")
load(file = "rdata/seurat.white20.postQC.RData")

# Visualizing clusters in snRNA-seq data
# See QC notebook to see details about clustering
# panel (A) in Fig. 5
UMAPPlot(scellenone,group.by="seurat_clusters",pt.size=1,cols=brewer.pal(n = 2,name = "Set1"),label=TRUE,label.size=5)+theme_void()+theme(legend.position = "top")+theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())
FeaturePlot(scellenone,"THY1",pt.size=1)+scale_color_viridis(option = "B")+theme_void(base_size = size)
FeaturePlot(scellenone,"ADIPOQ",pt.size=1)+scale_color_viridis(option = "B")+theme_void(base_size = size)

# Visualizing clusters in snRNA-seq data
# See QC notebook to see details about clustering
# Panel (A) in Fig. 5
UMAPPlot(swhite_20,group.by="seurat_clusters",pt.size=.5,cols=brewer.pal(n = 3,name = "Set1"),label=TRUE,label.size=5)+theme_void()+theme(legend.position = "top")+theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())
FeaturePlot(swhite_20,"THY1",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)
FeaturePlot(swhite_20,"ADIPOQ",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)

# Panel (B) in Fig. 5
Idents(scellenone)="seurat_clusters"
markers=FindAllMarkers(scellenone, only.pos = TRUE,assay = "RNA",slot = "data",logfc.threshold = .5,verbose = FALSE)
markers=subset.data.frame(markers,subset = p_val_adj<.05)
top.cells <- markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_logFC)
DoHeatmap(scellenone, features = top.cells$gene,group.colors =brewer.pal(n = 3,name = "Set1"),raster=FALSE)+scale_fill_viridis(option = "B")+theme(text = element_text(size = 10))

# Panel (C) in Fig. 5
Idents(swhite_20)="seurat_clusters"
swhite_20_sub=subset(swhite_20,cells = rownames(swhite_20@meta.data)[swhite_20@meta.data$seurat_clusters!=6])
markers=FindAllMarkers(swhite_20_sub, only.pos = TRUE,assay = "RNA",slot = "data",logfc.threshold = .5,verbose = FALSE)
markers=subset.data.frame(markers,subset = p_val_adj<.05)
top.nuclei <- markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_logFC)
swhite_20_sub=ScaleData(swhite_20_sub,features = top.nuclei$gene)
DoHeatmap(swhite_20_sub, features = top.nuclei$gene,group.colors =brewer.pal(n = 3,name = "Set1"),cells=sample(x = colnames(swhite_20_sub),size = 150),raster=FALSE)+scale_fill_viridis(option = "B")+theme(text = element_text(size = 10))

# GO analysis for cluster 0 and 1 in both scRNA-seq and snRNA-seq datasets
# scRNA-seq dataset
# Cluster 1
markers.1=FindMarkers(object = scellenone,ident.1 = 1,ident.2=0,assay = "RNA",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
markers.1=subset.data.frame(markers.1,subset = p_val_adj<.05)
write.table(x = rownames(markers.1),file = "additional_files/de_genes/celleone_adipo.txt",sep = "\n",row.names = FALSE,col.names = FALSE,quote = FALSE,eol = "\n")
# Plotting GO results. Results can be found in additional_files folder. geneontology.org was used for GO
celleone_adipo_GO <- read.delim("additional_files/gene ontology/results/celleone_adipo.txt", header=TRUE, comment.char="#")
celleone_adipo_GO$FDR.Plot=-log10(celleone_adipo_GO$upload_1..FDR)
ggplot(celleone_adipo_GO[1:10,],mapping = aes(x = reorder(GO.biological.process.complete, FDR.Plot), FDR.Plot))+geom_bar(stat = "identity",fill="darkred",color="black")+ coord_flip()+theme_bw(base_size = 15)+xlab("GO Biological Process")+ylab("-log(FDR)")

# Cluster 0
markers.0=FindMarkers(object = scellenone,ident.1 = 0,ident.2=1,assay = "RNA",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
markers.0=subset.data.frame(markers.0,subset = p_val_adj<.05)
write.table(x = rownames(markers.0),file = "additional_files/de_genes/celleone_pre.txt",sep = "\n",row.names = FALSE,col.names = FALSE,quote = FALSE,eol = "\n")
# Plotting GO results. Results can be found in additional_files folder. geneontology.org was used for GO
celleone_pre_GO <- read.delim("additional_files/gene ontology/results/celleone_pre.txt", header=TRUE, comment.char="#")
celleone_pre_GO$FDR.Plot=-log10(celleone_pre_GO$upload_1..FDR)
ggplot(celleone_pre_GO[1:10,],mapping = aes(x = reorder(GO.biological.process.complete, FDR.Plot), FDR.Plot))+geom_bar(stat = "identity",fill="darkred",color="black")+ coord_flip()+theme_bw(base_size = 15)+xlab("GO Biological Process")+ylab("-log(FDR)")

# snRNA-seq
# Cluster 1
i_2=FindMarkers(object = swhite_20,ident.1 = 1,ident.2=0,assay = "RNA",slot = "data",logfc.threshold = .5)
i_2=subset.data.frame(i_2,subset = p_val_adj<.05)
write.table(x = rownames(i_2),file = "additional_files/de_genes/white_20_adipo.txt",sep = "\n",row.names = F,col.names = F,quote = F,eol = "\n")
# Plotting GO results. Results can be found in additional_files folder. geneontology.org was used for GO
white_20_adipo_GO <- read.delim("additional_files/gene ontology/results/white_20_adipo.txt", header=TRUE, comment.char="#")
white_20_adipo_GO$FDR.Plot=-log10(white_20_adipo_GO$upload_1..FDR.)
library(ggplot2)
p1=ggplot(white_20_adipo_GO[1:10,],mapping = aes(x = reorder(GO.biological.process.complete, FDR.Plot), FDR.Plot))+geom_bar(stat = "identity",fill="#377EB8",color="black")+ coord_flip()+theme_classic(base_size = 15)+xlab("GO Biological Process")+ylab("-log(FDR)")

# Cluster 0
i_1=FindMarkers(object = swhite_20,ident.1 = 0,ident.2=1,assay = "RNA",slot = "data",logfc.threshold = .5)
i_1=subset.data.frame(i_1,subset = p_val_adj<.05)
write.table(x = rownames(i_1),file = "additional_files/de_genes/white_20_pre.txt",sep = "\n",row.names = F,col.names = F,quote = F,eol = "\n")
# Plotting GO results. Results can be found in additional_files folder. geneontology.org was used for GO
white_20_pre_GO <- read.delim("additional_files/gene ontology/results/white_20_pre.txt", header=TRUE, comment.char="#")
white_20_pre_GO$FDR.Plot=-log10(white_20_pre_GO$upload_1..FDR.)
library(ggplot2)
p2=ggplot(white_20_pre_GO[1:10,],mapping = aes(x = reorder(GO.biological.process.complete, FDR.Plot), FDR.Plot))+geom_bar(stat = "identity",fill="#E41E1C",color="black")+ coord_flip()+theme_classic(base_size = 15)+xlab("GO Biological Process")+ylab("-log(FDR)")

# What genes that are DE in scRNA-seq are also DE in snRNA-seq
adipo.genes <- rownames(markers.1)
adipo.genes=intersect(rownames(swhite_20),adipo.genes)
i_3=FindMarkers(object = swhite_20,features = adipo.genes,ident.1 = 1,ident.2=0,assay = "RNA",slot = "data",logfc.threshold = 0,verbose = FALSE,only.pos = TRUE)
i_3=subset.data.frame(i_3,subset = p_val_adj<.05)
adipo.not.detected=(setdiff(adipo.genes,intersect(rownames(i_3),adipo.genes)))

# Vision analysis 
# Perform transcriptional signature analysis using Vision. This analysis takes time, hence you could directly load the vision object using the command in line # 134  
# Read in expression counts (Genes X Cells)
counts <- swhite_20@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 <- swhite_20@reductions$umap@cell.embeddings
# Read in meta data (Cells x Vars)
meta = swhite_20@meta.data
# Read in signatures of cluster-0 and cluster-1 generated using scRNA-seq dataset
# Signatures were derived by performing DE test between cluster-0 and cluster-1 in scRNA-seq dataset, with enriched genes for respective cluster forming a signature
dir="additional_files/signatures/CellenOne/"
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.adipocyte.nuclei.RData")
# Or directly load Vision object
#load(file = "rdata/vision.adipocyte.nuclei.RData")

# recreating panel (D) and (E) in Fig. 5
swhite_20@meta.data$adipo=vis@SigScores[rownames(swhite_20@meta.data),"cellenone_adipo"]
FeaturePlot(swhite_20,"adipo",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)
swhite_20@meta.data$pre=vis@SigScores[,"cellenone_pre"]
FeaturePlot(swhite_20,"pre",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)
```

# Script for analyzing brown adipocyte day 20 dataset
---
title: "R Notebook"
output: html_notebook
---

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

```{r}
# load libraries
library(Seurat)

# Read cellenone seurat data
load(file = "rdata/seurat.cellenone.RData")
scellenone=seurat_singlet_object
rm(seurat_singlet_object)

# We initially ran the experiment for both white and brown adipocytes
# This Seurat object has brown adipocytes as well. However, doublets and negatives were removed based on the images of cells acquired during spotting.
# Therefore, sub-sample to white cells and extract the count and metadata for white cells
# Then, initiate a new Seurat object 
scellenone=subset(scellenone,Type=="Deep")
count=scellenone@assays$RNA@counts
meta=scellenone@meta.data
scellenone=CreateSeuratObject(counts = count,meta.data = meta)
scellenone=NormalizeData(scellenone)
scellenone=FindVariableFeatures(scellenone)
scellenone=ScaleData(scellenone)
scellenone=RunPCA(scellenone)
ElbowPlot(scellenone)

# Clustering
scellenone=RunUMAP(scellenone,dims = 1:20)
scellenone=FindNeighbors(object = scellenone,dims = 1:20)
scellenone=FindClusters(scellenone,resolution = 1.1)
UMAPPlot(scellenone,label=T)
```

```{r}
# Load VISION
library(VISION)
library(Matrix)
# Read in expression counts (Genes X Cells)
counts <- scellenone@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 <- scellenone@reductions$umap@cell.embeddings

# Read in meta data (Cells x Vars)
meta = scellenone@meta.data

# Read in signatures
Signatures=c("/Users/anushka/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/additional_files/signatures/geneset.gmt","/Users/anushka/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/additional_files/signatures/Nuclei_White_Brown/1.cells.gmt","/Users/anushka/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/additional_files/signatures/Nuclei_White_Brown/2.cells.gmt")

#
vis <- Vision(scaled_counts,signatures = Signatures,meta = meta)
vis <- addProjection(vis, "UMAP", projection)
vis <- analyze(vis)
viewResults(vis)
scellenone$adipo=vis@SigScores[,1]
FeaturePlot(scellenone,"adipo",pt.size=3)+scale_color_viridis()+theme_void()
UMAPPlot(scellenone,cells.highlight=WhichCells(scellenone,idents = 3),pt.size=3,sizes.highlight = 3,cols = "black")+theme_void()+theme(legend.position = "top")
```
```{r}
library(dplyr)
scellenone@meta.data %>% group_by(seurat_clusters) %>% summarise(mean(adipo))
```

```{r}
#df=data.frame(score=c(vis@SigScores[WhichCells(scellenone,idents = 3),1],vis@SigScores[WhichCells(scellenone,idents = 3),2],vis@SigScores[WhichCells(scellenone,idents = 3),3]),sig=c(replicate(n = 26,expr = "adipo"),replicate(n = 26,expr = "p1"),replicate(n = 26,expr = "p2")))
df=data.frame(score=c(vis@SigScores[WhichCells(scellenone,idents = 3),2],vis@SigScores[WhichCells(scellenone,idents = 3),3]),sig=c(replicate(n = 26,expr = "P1"),replicate(n = 26,expr = "P2")))
ggplot(df,aes(x=sig,y=score))+geom_boxplot()+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.75)+theme_classic(base_size = 20)+xlab("Cluster")+ylab("Signature Score")
wilcox.test(score~sig,data=df)
```
```{r}
FeaturePlot(scellenone,"S100A4",pt.size=2)+scale_color_viridis()+theme_void()
ggplot(as.data.frame(t(scellenone@assays$RNA@data)),aes(x=ADIPOQ,y=SEMA5A))+geom_point(size=3,alpha=0.3)+theme_classic(base_size = 20)
ggplot(as.data.frame(t(scellenone@assays$RNA@data)),aes(x=ADIPOQ,y=S100A4))+geom_point(size=3,alpha=0.3)+theme_classic(base_size = 20)
```



