integrated_create_seurat.Rmd

```{r eval=FALSE}

library(Seurat)
library(tidyverse)
library(cowplot)
library(Matrix)
library(viridis)
library(harmony)
library(ggpubr)
library(patchwork)
library(RColorBrewer)
library(ggrepel)
library(CellTrek)

setwd("/dfs7/swaruplab/smorabit/analysis/ADDS_2021/splitseq/integration/")

fig_dir <- "figures/"
data_dir <- "data/"

umap_theme <- theme(
  axis.line=element_blank(),
  axis.text.x=element_blank(),
  axis.text.y=element_blank(),
  axis.ticks=element_blank(),
  axis.title.x=element_blank(),
  axis.title.y=element_blank(),
  panel.background=element_blank(),
  panel.border=element_blank(),
  panel.grid.major=element_blank(),
  panel.grid.minor=element_blank(),
  plot.background=element_blank()
)



# load old split-seq object
seurat_adds <- readRDS("/dfs7/swaruplab/smorabit/analysis/ADDS_2021/splitseq/data/ADDS_integrated_scvi_old.rds")


# load other datasets
seurat_zhou <- readRDS("/dfs7/swaruplab/smorabit/analysis/scWGCNA/data/Zhou_2020.rds")
seurat_morabito <- readRDS("/dfs7/swaruplab/smorabit/analysis/scWGCNA/data/Swarup_2021.rds")
seurat_mathys <- readRDS("/dfs7/swaruplab/smorabit/analysis/scWGCNA/data/Mathys_2019.rds")

# load cell metadata
cell_meta <- read.delim("/dfs7/swaruplab/smorabit/analysis/ADDS_2021/splitseq/integration/data/ADDS_subclusters_final.tsv", sep='\t', header=TRUE)
cell_meta$bc <- substr(cell_meta$X, 1, nchar(cell_meta$X) - 2)
rownames(cell_meta) <- cell_meta$X
cell_meta <- dplyr::select(cell_meta, -c(X))

# subset mathys
tmp <- subset(cell_meta, Study == 'Mathys_2019')
seurat_mathys <- seurat_mathys[,colnames(seurat_mathys) %in% tmp$bc]
seurat_mathys <- seurat_mathys[,tmp$bc]
seurat_mathys@meta.data <- tmp

# subset zhou
tmp <- subset(cell_meta, Study == 'Zhou_2020')
seurat_zhou <- seurat_zhou[,colnames(seurat_zhou) %in% tmp$bc]
seurat_zhou <- seurat_zhou[,tmp$bc]
meta <- seurat_zhou@meta.data
seurat_zhou@meta.data <- tmp
all.equal(meta$cell_type, tmp$cell_type)
all.equal(meta$group, tmp$Diagnosis)


# subset morabito
tmp <- subset(cell_meta, Study == 'Morabito')
seurat_morabito <- seurat_morabito[,colnames(seurat_morabito) %in% tmp$bc]
seurat_morabito <- seurat_morabito[,tmp$bc]
seurat_morabito@meta.data <- tmp

# subset ADDS
tmp <- subset(cell_meta, Study == 'ADDS')
seurat_adds <- seurat_adds[,colnames(seurat_adds) %in% tmp$bc]
seurat_adds <- seurat_adds[,tmp$bc]
seurat_adds@meta.data <- tmp


####################################################################
# set up separate ADDS and AD objects
#
# * need to have two separate objects because some of the genes are
#   not in both datasets, but we want all of the genes for the
#   control vs condition DEG tests etc...
####################################################################


# add UMAP:
umap <- dplyr::select(seurat_adds@meta.data, c(UMAP_1, UMAP_2)) %>% as.matrix
seurat_adds@reductions$umap <- Seurat::CreateDimReducObject(
  embeddings = umap,
  key="UMAP",
  assay="RNA"
)

# add broad cell type info
cell_ids <- read.csv(paste0(data_dir, 'integrated_subclusters.csv'))

ix <- match(seurat_morabito$cell_identity, cell_ids$cell_identity)
seurat_morabito$cell_type <- cell_ids$cell_type[ix]

ix <- match(seurat_zhou$cell_identity, cell_ids$cell_identity)
seurat_zhou$cell_type <- cell_ids$cell_type[ix]

ix <- match(seurat_mathys$cell_identity, cell_ids$cell_identity)
seurat_mathys$cell_type <- cell_ids$cell_type[ix]

# merge AD datasets
seurat_AD <- merge(x = seurat_morabito, y = c(seurat_mathys, seurat_zhou))

# fix study:
seurat_AD$Study <- ifelse(
  seurat_AD$Study == 'Morabito','Morabito & Miyoshi', seurat_AD$Study
)
seurat_AD$Study <- ifelse(
  seurat_AD$Study == 'Morabito & Miyoshi','Morabito_Miyoshi', seurat_AD$Study
)
seurat_obj$Study <- ifelse(
  seurat_obj$Study == 'Morabito & Miyoshi','Morabito_Miyoshi', seurat_obj$Study
)

# add UMAP:
umap <- dplyr::select(seurat_AD@meta.data, c(UMAP_1, UMAP_2)) %>% as.matrix
seurat_AD@reductions$umap <- Seurat::CreateDimReducObject(
  embeddings = umap,
  key="UMAP",
  assay="RNA"
)

# fix diagnosis
tmp <- seurat_AD$Diagnosis
seurat_AD$Diagnosis <- ifelse(seurat_AD$Diagnosis == 'Alzheimer Disease','AD', seurat_AD$Diagnosis)
seurat_AD$Diagnosis <- ifelse(seurat_AD$Diagnosis == 'P','AD', seurat_AD$Diagnosis)
seurat_AD$Diagnosis <- ifelse(seurat_AD$Diagnosis == 'C','Control', seurat_AD$Diagnosis)
seurat_AD$Diagnosis <- ifelse(seurat_AD$Diagnosis == 'Other','Early-AD', seurat_AD$Diagnosis)

# get rid of NA cells:
# seurat_AD <- seurat_AD[,!is.na(seurat_AD$UMAP_1)]
# any(is.na(seurat_AD$UMAP_1))

any(is.na(seurat_adds$UMAP_1))

plot_list <- PlotEmbedding(
  seurat_AD,
  group.by = 'cell_identity',
  split.by = 'Diagnosis',
  raster_dpi = 100,
  plot_theme = umap_theme + NoLegend(),
  plot_under=TRUE
)

pdf(paste0(fig_dir, 'umap_seurat_AD_integrated.pdf'), width=18, height=6)
wrap_plots(plot_list, ncol=3) + plot_layout(guides='collect')
dev.off()

# write seurat obj
saveRDS(seurat_AD, file="/dfs7/swaruplab/smorabit/analysis/ADDS_2021/splitseq/integration/data/AD_integrated.rds" )
saveRDS(seurat_adds, file="/dfs7/swaruplab/smorabit/analysis/ADDS_2021/splitseq/integration/data/ADDS_integrated.rds" )


####################################################################
# set up joint ADDS + AD object
####################################################################

genes.keep <- intersect(rownames(seurat_adds), rownames(seurat_AD))

seurat_adds <- seurat_adds[genes.keep,]
seurat_AD <- seurat_AD[genes.keep,]


# merge into one big Seurat object
# TODO: check if normalization is preserved?
seurat_obj <- merge(seurat_AD,y=seurat_adds)


ix <- match(seurat_obj$cell_identity, cell_ids$cell_identity)
seurat_obj$cell_type <- cell_ids$cell_type[ix]


# set factor levels:
seurat_obj$Diagnosis <- factor(
  as.character(seurat_obj$Diagnosis),
  levels = c('Control', 'Early-AD', 'AD', 'DSAD')
)
table(seurat_obj$Diagnosis)

# add UMAP:
umap <- dplyr::select(seurat_obj@meta.data, c(UMAP_1, UMAP_2)) %>% as.matrix

seurat_obj@reductions$umap <- Seurat::CreateDimReducObject(
  embeddings = umap,
  key="UMAP",
  assay="RNA"
)

p <- PlotEmbedding(
  seurat_obj,
  group.by = 'cell_identity',
  point_size=0.0001,
  raster_dpi=200,
  plot_theme = umap_theme + NoLegend()
)

pdf(paste0(fig_dir, 'umap_seurat_integrated2.pdf'), width=10, height=10)
p
dev.off()

p <- PlotEmbedding(
  seurat_obj,
  group.by = 'cell_identity',
  point_size=0.0001,
  raster=FALSE,
  plot_theme = umap_theme + NoLegend()
)

pdf(paste0(fig_dir, 'umap_seurat_integrated_noraster.pdf'), width=10, height=10)
p
dev.off()

seurat_obj@meta.data %>% subset(cell_type == 'FBR') %>% .$Diagnosis %>% table
seurat_obj@meta.data %>% subset(cell_type == 'PER') %>% .$Diagnosis %>% table
seurat_obj@meta.data %>% subset(cell_type == 'END') %>% .$Diagnosis %>% table



cur_groups <- c('ASC1', 'MG2')
plot_list <- PlotEmbedding(
  seurat_obj,
  group.by = 'cell_type',
  #selected = cur_groups,
  split.by = 'Study',
  raster_dpi = 100,
  plot_theme = umap_theme + NoLegend(),
  plot_under=TRUE
)

pdf(paste0(fig_dir, 'umap_seurat_integrated.pdf'), width=16, height=16)
wrap_plots(plot_list, ncol=2) + plot_layout(guides='collect')
dev.off()



##############################################################
# load scANVI
############################################################e#

# load cell metadata
cell_meta <- read.delim("/dfs7/swaruplab/smorabit/analysis/ADDS_2021/splitseq/integration/data/ADDS_subclusters_final.tsv", sep='\t', header=TRUE)
cell_meta$bc <- substr(cell_meta$X, 1, nchar(cell_meta$X) - 2)
rownames(cell_meta) <- cell_meta$X
cell_meta <- dplyr::select(cell_meta, -c(X))
scanvi <- read.csv(paste0(data_dir, 'ADDS_AD_scANVI.csv'), row.names=1) %>% as.matrix

ix <- match(seurat_obj$bc, cell_meta$bc)
length(ix)
dim(scanvi)
scanvi <- scanvi[ix,]
all.equal(as.character(cell_meta$bc[ix]), as.character(seurat_obj$bc))
dim(scanvi)
rownames(scanvi) <- colnames(seurat_obj)

seurat_obj@reductions$scanvi <- Seurat::CreateDimReducObject(
  embeddings = scanvi,
  key="SCANVI",
  assay="RNA"
)




# write seurat obj
saveRDS(seurat_obj, file="/dfs7/swaruplab/smorabit/analysis/ADDS_2021/splitseq/integration/data/ADDS_AD_integrated.rds")

p <- PlotEmbedding(
  seurat_obj,
  split.by = 'Study',
  group.by = 'Sample',
  point_size=0.05,
  raster_scale=0.5,
  label=FALSE,
  plot_theme = umap_theme + NoLegend()
)

pdf(paste0(fig_dir, 'umap_seurat_integrated_samples.pdf'), width=10, height=10)
wrap_plots(p, ncol=2)
dev.off()



####################################################################
# set up the ADDS
#
# was having some weird issues where subseting the ADDS object
# was insanely slow.
####################################################################


# re-load ADDS seurat object:
seurat_adds <- readRDS("/dfs7/swaruplab/smorabit/analysis/ADDS_2021/splitseq/data/ADDS_integrated_scvi_old.rds")
tic(); dim(subset(seurat_adds, cell_type == 'ASC')); toc() # 8.52 seconds

seurat_adds_new <- readRDS(file="/dfs7/swaruplab/smorabit/analysis/ADDS_2021/splitseq/integration/data/ADDS_integrated.rds" )

seurat_adds <- seurat_adds[, colnames(seurat_adds_new)]
tic(); dim(subset(seurat_adds, cell_type == 'ASC')); toc()

seurat_adds$cell_type <- seurat_adds_new$cell_type
seurat_adds$cell_identity <- seurat_adds_new$cell_identity
seurat_adds$UMAP_1 <- seurat_adds_new$UMAP_1
seurat_adds$UMAP_2 <- seurat_adds_new$UMAP_2
seurat_adds$Tissue <- seurat_adds_new$Tissue
seurat_adds$Sample <- seurat_adds_new$Sample
seurat_adds$Diagnosis <- seurat_adds_new$Diagnosis

tic(); dim(subset(seurat_adds, cell_type == 'ASC')); toc()


umap <- dplyr::select(seurat_adds@meta.data, c(UMAP_1, UMAP_2)) %>% as.matrix
seurat_adds@reductions$umap <- Seurat::CreateDimReducObject(
  embeddings = umap,
  key="UMAP",
  assay="RNA"
)
tic(); dim(subset(seurat_adds, cell_type == 'ASC')); toc()


saveRDS(seurat_adds, file="/dfs7/swaruplab/smorabit/analysis/ADDS_2021/splitseq/integration/data/ADDS_integrated.rds" )




```