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

# This notebook has commands to recreate results discussed in Results-section-5: Integration of snRNA-seq and scRNA-seq datasets 
```{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

# Four datasets were integrated using scVI
#1. scRNA-seq white preadipocyte
#2. snRNA-seq white preadipocyte
#3. scRNA-seq white adipocyte
#4. snRNA-seq white adipocyte

# Read seurat objects for each of these 4 datasets
load(file = "rdata/seurat.preadipocytes.post.demultiplexingAndQC.RData")
white_cell_singlets=read.delim("additional_files/barcodes/singlets.high.quality.white.cells.txt",header = FALSE)
white_cells=subset(dataset.subset,cells = white_cell_singlets$V1)

white_nuclei_0=ReadH5AD('additional_files/adata/white_00.h5ad')
white_nuclei_singlets = read.delim("additional_files/barcodes/singlets.high.quality.white.nuclei.txt",header=FALSE)
white_nuclei_0=subset(white_nuclei_0,cells = white_nuclei_singlets$V1)

load(file = "rdata/seurat.cellenone.postQC.RData")
load(file = "rdata/seurat.white20.postQC.RData")

# Merge all 4 seurat objects
swhite_0_20=merge(white_cells, y = c(scellenone, white_nuclei_0,swhite_20), add.cell.ids = c("0", "1", "2","3"))
# Add batch level information
swhite_0_20$batch=substr(colnames(swhite_0_20),start = 1,stop = 1)
Idents(swhite_0_20)="batch"

# Add latent space inferred from scVI with NO batch correction and generate UMAP
lwhite_0_20 <- read.csv("additional_files/latent_space/no_integration.csv", row.names=1)
lwhite_0_20=as.matrix(lwhite_0_20)
rownames(lwhite_0_20)=rownames(swhite_0_20@meta.data)
swhite_0_20[["scVI"]] <- CreateDimReducObject(embeddings =lwhite_0_20,key = "scVI_")
swhite_0_20=RunUMAP(swhite_0_20,reduction="scVI",dims = 1:10)
UMAPPlot(swhite_0_20)

# Add cell-type information 
white_cell_singlets=white_cell_singlets$V1
white_nuclei_singlets=white_nuclei_singlets$V1
white_cell_singlets=paste("0",white_cell_singlets,sep = "_")
white_nuclei_singlets=paste("2",white_nuclei_singlets,sep = "_")
white.cells.20.0=paste("1",rownames(scellenone@meta.data)[scellenone@meta.data$seurat_clusters==0],sep = "_")
white.cells.20.1=paste("1",rownames(scellenone@meta.data)[scellenone@meta.data$seurat_clusters==1],sep = "_")
white.nuclei.20.0=paste("3",rownames(swhite_20@meta.data)[swhite_20@meta.data$seurat_clusters==0],sep = "_") #0
white.nuclei.20.1=paste("3",rownames(swhite_20@meta.data)[swhite_20@meta.data$seurat_clusters==1],sep = "_") #1
white.nuclei.20.6=paste("3",rownames(swhite_20@meta.data)[swhite_20@meta.data$seurat_clusters==6],sep = "_") #2
swhite_0_20@meta.data$type=factor("day.0.cells.pread",levels = c("day.0.cells.pread","day.0.nuclei.pread","day.20.cells.pread","day.20.cells.adipo","day.20.nuclei.pread","day.20.nuclei.adipo","day.20.nuclei.mitotic"))
swhite_0_20@meta.data[white_nuclei_singlets,"type"]="day.0.nuclei.pread"
swhite_0_20@meta.data[white.cells.20.0,"type"]="day.20.cells.pread"
swhite_0_20@meta.data[white.cells.20.1,"type"]="day.20.cells.adipo"
swhite_0_20@meta.data[white.nuclei.20.0,"type"]="day.20.nuclei.pread"
swhite_0_20@meta.data[white.nuclei.20.1,"type"]="day.20.nuclei.adipo"
swhite_0_20@meta.data[white.nuclei.20.6,"type"]="day.20.nuclei.mitotic"
Idents(swhite_0_20)="type"
UMAPPlot(swhite_0_20,label=TRUE,label.size=6,pt.size=1)+theme_classic(base_size = 20)+theme(legend.position = "top")

# Plot cluster tree without integration
swhite_0_20@reductions$pca=swhite_0_20@reductions$scVI
swhite_0_20=BuildClusterTree(object = swhite_0_20,reorder = TRUE,dims = 1:10)
PlotClusterTree(swhite_0_20)
```
```{r}
# Performing same analysis but WITH integration
# Add latent space inferred from scVI with batch correction and generate UMAP
lwhite_0_20 <- read.csv("additional_files/latent_space/with_integration.csv", row.names=1)
lwhite_0_20=as.matrix(lwhite_0_20)
rownames(lwhite_0_20)=rownames(swhite_0_20@meta.data)
swhite_0_20[["scVI"]] <- CreateDimReducObject(embeddings =lwhite_0_20,key = "scVI_")
swhite_0_20=RunUMAP(swhite_0_20,reduction="scVI",dims = 1:10)
UMAPPlot(swhite_0_20)

# Plot cluster tree
swhite_0_20@reductions$pca=swhite_0_20@reductions$scVI
swhite_0_20=BuildClusterTree(object = swhite_0_20,reorder = TRUE,dims = 1:10)
PlotClusterTree(swhite_0_20)
```

```{r}
# Investigating localization of preadipocytes from both day-0-snRNA-seq and day-0-scRNA-seq datasets into two distinct groups
top=intersect(rownames(swhite_0_20@reductions$umap@cell.embeddings)[swhite_0_20@reductions$umap@cell.embeddings[,2]>1.25],c(white_cell_singlets,white_nuclei_singlets))
UMAPPlot(swhite_0_20,cells.highlight=top)
bottom=setdiff(c(white_cell_singlets,white_nuclei_singlets),top)
UMAPPlot(swhite_0_20,cells.highlight=bottom)

# Manually assigning clusters
swhite_0_20@meta.data$manual_cluster=factor(0,levels = c(0,1,2,3,4))
swhite_0_20@meta.data[WhichCells(object = swhite_0_20,idents = c("day.20.nuclei.adipo","day.20.cells.adipo")),"manual_cluster"]=1
swhite_0_20@meta.data[top,"manual_cluster"]=2
swhite_0_20@meta.data[bottom,"manual_cluster"]=3
swhite_0_20@meta.data[WhichCells(object = swhite_0_20,idents = "day.20.nuclei.mitotic"),"manual_cluster"]=4
Idents(swhite_0_20)="manual_cluster"
UMAPPlot(swhite_0_20)

s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes
swhite_0_20 <- CellCycleScoring(swhite_0_20, s.features = s.genes, g2m.features = g2m.genes, set.ident = F)
UMAPPlot(swhite_0_20,label=TRUE,label.size=6,pt.size=1.25,cols=brewer.pal(n = 4,name = "Set3"),group.by="Phase")+theme_classic(base_size = 20)

UMAPPlot(swhite_0_20,label=TRUE,label.size=5,pt.size=.5,cols=brewer.pal(n = 5,name = "Dark2"))+theme_void()+theme(legend.position = "none")+theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank(),legend.text = element_text(size=size))
VlnPlot(swhite_0_20,"PLK1",pt.size=0,cols = brewer.pal(n = 5,name = "Dark2"))+theme_classic(base_size =size)
p3=VlnPlot(swhite_0_20,"MYBL2",pt.size=0,cols = brewer.pal(n = 5,name = "Dark2"))+theme_classic(base_size =size)+theme(legend.position = "none")
p4=VlnPlot(swhite_0_20,"BUB1",pt.size=0,cols = brewer.pal(n = 5,name = "Dark2"))+theme_classic(base_size =size)+theme(legend.position = "none")
p5=VlnPlot(swhite_0_20,"MKI67",pt.size=0,cols = brewer.pal(n = 5,name = "Dark2"))+theme_classic(base_size =size)+theme(legend.position = "none")
p6=VlnPlot(swhite_0_20,"CDK1",pt.size=0,cols = brewer.pal(n = 5,name = "Dark2"))+theme_classic(base_size =size)+theme(legend.position = "none")
p7=VlnPlot(swhite_0_20,"CCNB1",pt.size=0,cols = brewer.pal(n = 5,name = "Dark2"))+theme_classic(base_size =size)+theme(legend.position = "none")
```

# Seurat Integration
---
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)
library(Matrix)
library(ggplot2)
library(RColorBrewer)
library(viridis)
library(tibble)
library(clusterCrit)
library(dplyr)
library(ComplexHeatmap)
library("circlize")
library(DropletUtils)
library(VISION)

size=10

# Four datasets were integrated using scVI
#1. scRNA-seq white preadipocyte
#2. snRNA-seq white preadipocyte
#3. scRNA-seq white adipocyte
#4. snRNA-seq white adipocyte

# Read seurat objects for each of these 4 datasets
load(file = "rdata/seurat.preadipocytes.post.demultiplexingAndQC.RData")
white_cell_singlets=read.delim("additional_files/barcodes/singlets.high.quality.white.cells.txt",header = FALSE)
white_cells=subset(dataset.subset,cells = white_cell_singlets$V1)

white_nuclei_0=ReadH5AD('additional_files/adata/white_00.h5ad')
white_nuclei_singlets = read.delim("additional_files/barcodes/singlets.high.quality.white.nuclei.txt",header=FALSE)
white_nuclei_0=subset(white_nuclei_0,cells = white_nuclei_singlets$V1)

load(file = "rdata/seurat.cellenone.postQC.RData")
load(file = "rdata/seurat.white20.postQC.RData")
```

```{r}
ifnb.list=vector(mode = "list",length = 4)
ifnb.list[[1]]=white_cells
ifnb.list[[2]]=scellenone
ifnb.list[[3]]=white_nuclei_0
ifnb.list[[4]]=swhite_20
save(ifnb.list,file = "~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/revisions/rdata/ifnb.list.RData")
# The ifnb.list object was saved and integrated on AWS EC@ using the code chunk below.
```


```{r}
# This code chunk was runk on AWS EC2 with t3.2xlarge
library(Seurat)
load("~/ifnb.list.RData")
# normalize and identify variable features for each dataset independently
ifnb.list <- lapply(X = ifnb.list, FUN = function(x) {
    DefaultAssay(object = x) <- "RNA"
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = ifnb.list)
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features,k.filter = 147)
immune.combined <- IntegrateData(anchorset = immune.anchors)

DefaultAssay(immune.combined) <- "integrated"
# Run the standard workflow for visualization and clustering
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindClusters(immune.combined, resolution = 0.5)
UMAPPlot(immune.combined,label=T)
save(immune.combined,file = "~/immune.combined.RData")
```

```{r}
load(file = "~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/revisions/rdata/immune.combined.RData")
immune.combined$batch=substr(colnames(immune.combined),start = nchar(colnames(immune.combined)),stop = nchar(colnames(immune.combined)))
Idents(immune.combined)="batch"
library(ggplot2)
UMAPPlot(immune.combined,label=T,pt.size=.5)+theme_void() # Bring 2 on top
immune.combined$sample=factor("cells",levels = c("nuclei","cells"))
immune.combined@meta.data[WhichCells(immune.combined,idents = c(3,4)),"sample"]="nuclei"
UMAPPlot(immune.combined,label=FALSE,label.size=10,pt.size=.75,group.by="sample",cols=c("#CBC9E2","#756BB1"))+theme_void()+theme(legend.position = "none")
immune.combined$day=factor(0,levels = c(0,20))
immune.combined@meta.data[WhichCells(immune.combined,idents = c(2,4)),"day"]=20
UMAPPlot(immune.combined,label=FALSE,label.size=10,pt.size=.75,group.by="day",cols=c("#BDD7E7","#2171B5"))+theme_void()+theme(legend.position = "none")
```

```{r}
# Add cell-type information 
white_cell_singlets=white_cell_singlets$V1
white_nuclei_singlets=white_nuclei_singlets$V1
white_cell_singlets=paste(white_cell_singlets,"1",sep = "_")
white_nuclei_singlets=paste(white_nuclei_singlets,"3",sep = "_")
white.cells.20.0=paste(rownames(scellenone@meta.data)[scellenone@meta.data$seurat_clusters==0],"2",sep = "_")
white.cells.20.1=paste(rownames(scellenone@meta.data)[scellenone@meta.data$seurat_clusters==1],"2",sep = "_")
white.nuclei.20.0=paste(rownames(swhite_20@meta.data)[swhite_20@meta.data$seurat_clusters==0],"4",sep = "_") #0
white.nuclei.20.1=paste(rownames(swhite_20@meta.data)[swhite_20@meta.data$seurat_clusters==1],"4",sep = "_") #1
white.nuclei.20.6=paste(rownames(swhite_20@meta.data)[swhite_20@meta.data$seurat_clusters==6],"4",sep = "_") #2
immune.combined@meta.data$type=factor("day.0.cells.pread",levels = c("day.0.cells.pread","day.0.nuclei.pread","day.20.cells.pread","day.20.cells.adipo","day.20.nuclei.pread","day.20.nuclei.adipo","day.20.nuclei.mitotic"))
immune.combined@meta.data[white_nuclei_singlets,"type"]="day.0.nuclei.pread"
immune.combined@meta.data[white.cells.20.0,"type"]="day.20.cells.pread"
immune.combined@meta.data[white.cells.20.1,"type"]="day.20.cells.adipo"
immune.combined@meta.data[white.nuclei.20.0,"type"]="day.20.nuclei.pread"
immune.combined@meta.data[white.nuclei.20.1,"type"]="day.20.nuclei.adipo"
immune.combined@meta.data[white.nuclei.20.6,"type"]="day.20.nuclei.mitotic"
Idents(immune.combined)="type"
UMAPPlot(immune.combined,label=TRUE,label.size=6,pt.size=.5)+theme_classic(base_size = 20)+theme(legend.position = "top")

# Plot cluster tree without integration
#immune.combined@reductions$pca=immune.combined@reductions$scVI
immune.combined=BuildClusterTree(object = immune.combined,reorder = TRUE,dims = 1:30)
PlotClusterTree(immune.combined)
```

```{r}
library(RColorBrewer)
library(ggplot2)
Idents(immune.combined)="type"
x=as.data.frame(immune.combined@reductions$umap@cell.embeddings[WhichCells(object = immune.combined,idents="day.20.cells.adipo"),])
y=as.data.frame(immune.combined@reductions$umap@cell.embeddings[WhichCells(object = immune.combined,idents="day.0.cells.pread"),])
#z=as.data.frame(immune.combined@reductions$umap@cell.embeddings[WhichCells(object = immune.combined,idents="day.20.nuclei.adipo"),])
#z2=as.data.frame(immune.combined@reductions$umap@cell.embeddings[WhichCells(object = immune.combined,idents="day.20.nuclei.pread"),])
#z3=as.data.frame(immune.combined@reductions$umap@cell.embeddings[WhichCells(object = immune.combined,idents="day.20.nuclei.mitotic"),])
z4=as.data.frame(immune.combined@reductions$umap@cell.embeddings[WhichCells(object = immune.combined,idents="day.20.cells.pread"),])

day.0.cells.pread="#1B9E77"
day.0.nuclei.pread="#D95F02"
day.20.cells.pread="#7570B3"
day.20.cells.adipo="#E7298A"
day.20.nuclei.pread="#66A61E"
day.20.nuclei.adipo="#E6AB02"
day.20.nuclei.mitotic="#A6761D"
UMAPPlot(immune.combined,label=FALSE,label.size=10,pt.size=.75,cols=c(day.0.cells.pread,day.0.nuclei.pread,day.20.cells.pread,day.20.cells.adipo,day.20.nuclei.pread,day.20.nuclei.adipo,day.20.nuclei.mitotic))+theme_void()+theme(legend.position = "top")+geom_point(data=y,aes(x=UMAP_1,y=UMAP_2),color="#1B9E77",size=.75)+geom_point(data=x,aes(x=UMAP_1,y=UMAP_2),color="#E7298A",size=.75)+geom_point(data=z4,aes(x=UMAP_1,y=UMAP_2),color="#7570B3",size=.75)+theme(legend.position = "none")#+geom_point(data=z,aes(x=UMAP_1,y=UMAP_2),color="#E6AB02",size=1.5)+geom_point(data=z2,aes(x=UMAP_1,y=UMAP_2),color="#66A61E",size=1.5)+geom_point(data=z3,aes(x=UMAP_1,y=UMAP_2),color="#A6761D",size=1.5)
```

# 2 batches only
# One with cells and one with nuclei 
```{r}
# This code chunk was runk on AWS EC2 with t3.2xlarge
load("~/ifnb.list.RData")
x=merge(ifnb.list[[1]],ifnb.list[[2]],add.cell.ids = c("1","2")) # Cells as one batch
y=merge(ifnb.list[[3]],ifnb.list[[4]],add.cell.ids = c("3","4")) # Nuclei as one batch
ifnb.list=vector(mode = "list",length = 2) 
ifnb.list[[1]]=x
ifnb.list[[2]]=y
# normalize and identify variable features for each dataset independently
ifnb.list <- lapply(X = ifnb.list, FUN = function(x) {
    DefaultAssay(object = x) <- "RNA"
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = ifnb.list)
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features,k.filter = 147)
immune.combined <- IntegrateData(anchorset = immune.anchors)

DefaultAssay(immune.combined) <- "integrated"
# Run the standard workflow for visualization and clustering
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindClusters(immune.combined, resolution = 0.5)
UMAPPlot(immune.combined,label=T)

immune.combined$batch=substr(colnames(immune.combined),start = 1,stop = 1)
Idents(immune.combined)="batch"
UMAPPlot(immune.combined,label=T) # the immune.combined object can be loaded by load(file = "~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/revisions/rdata/immune.combined.2batches.cells.nuclei.RData")
```

# 2 batches only
# One with day0cells and day20 nuclei and one with day0 nuclei and day20 cells 
```{r}
library(Seurat)
load("~/ifnb.list.RData")
x=merge(ifnb.list[[1]],ifnb.list[[4]],add.cell.ids = c("1","4")) # Cells as one batch
y=merge(ifnb.list[[3]],ifnb.list[[2]],add.cell.ids = c("3","2")) # Nuclei as one batch
ifnb.list=vector(mode = "list",length = 2) 
ifnb.list[[1]]=x
ifnb.list[[2]]=y
# normalize and identify variable features for each dataset independently
ifnb.list <- lapply(X = ifnb.list, FUN = function(x) {
    DefaultAssay(object = x) <- "RNA"
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = ifnb.list)
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features,k.filter = 147)
immune.combined <- IntegrateData(anchorset = immune.anchors)

DefaultAssay(immune.combined) <- "integrated"
# Run the standard workflow for visualization and clustering
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)
immune.combined$batch=substr(colnames(immune.combined),start = 1,stop = 1)
Idents(immune.combined)="batch"
UMAPPlot(immune.combined,label=T) # the immune.combined object can be loaded by load(file = "~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/revisions/rdata/immune.combined.2batches.RData")
```