---
title: "jingyeST_Exp_Analysis"
author: "Vinh Ton"
date: "8/14/2023"
output: html_document
---
Exploratory analysis of Jingye's ST data ST359399
Following steps from https://satijalab.org/seurat/articles/spatial_vignette#integration-with-single-cell-data

Other useful vignettes
https://satijalab.org/seurat/articles/future_vignette

```{r setup, include=FALSE}
# note that in rmd setwd is not persistent through chunks in .Rmd.
# https://bookdown.org/yihui/rmarkdown-cookbook/working-directory.html
# instead use knitr::opts_knit$set(root.dir = '/tmp') to set a persistent working directory
# NOTE that this doesn't modify the console wd!
knitr::opts_knit$set(root.dir = '/gpfs/data/yanailab/projects/vt2217/vinh_202308_jingye_ST/r_scripts/')

library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
# https://patchwork.data-imaginist.com/articles/guides/annotation.html
library(dplyr)
library(tidyr)
library(future)
# BiocManager::install("clusterProfiler")
# BiocManager::install("EnhancedVolcano")
# BiocManager::install("org.Mm.eg.db")
# BiocManager::install("MAST")
library(clusterProfiler)
library(org.Mm.eg.db)
library(msigdbr)
library(cowplot)
library(harmony)
library(MAST)
library(glmGamPoi)
library(raster)
library(RColorBrewer)
# modify the statement below to match the number of cores in my r session
plan("multisession", workers = 40)
options(future.globals.maxSize = 16000 * 1024^2)  # Set the limit to 16 GB, increases the maximum size for https://nbisweden.github.io/workshop-scRNAseq/labs/compiled/seurat/seurat_07_spatial.html#Integration
# globals exported for the future package.

#make sure you have at least 256 gb ram to run this whole code safely!

set.seed(100)

scripts_path = '/gpfs/data/yanailab/projects/vt2217/vinh_202308_jingye_ST/r_scripts/' # goes to the project directory
data_path = '/gpfs/data/yanailab/projects/vt2217/vinh_202308_jingye_ST/data/' # goes to raw data
save_path = '/gpfs/data/yanailab/projects/vt2217/vinh_202308_jingye_ST/output/' # go to new folder or something that you've already made

knitr::opts_chunk$set(echo = TRUE)
```

# Loading in data
```{r}
noMALDI.list <- list() # stores seurat objects for 399
# loading in 399- samples (no MALDI)
for (x in c("A","B","C","D")){
  noMALDI.list[[paste0(x,'1','_399')]]<-Load10X_Spatial(data.dir = paste0(data_path, 'ST359399/399-',x,"1",'/outs'))
}

postMALDI.list <- list() # stores seurat objects for 359
# loading in 359- samples (post-MALDI)
for (x in c("A","B","C","D")){
  postMALDI.list[[paste0(x,'1','_359')]]<-Load10X_Spatial(data.dir = paste0(data_path, 'ST359399/359-',x,"1",'/outs'))
}
```

# Data preprocessing, filtering out high percentage mitochondrial data and low UMI counts
```{r}
# noMALDI 399 data
for (x in c("A","B","C","D")){
  noMALDI.list[[paste0(x,'1','_399')]]<-PercentageFeatureSet(noMALDI.list[[paste0(x,'1','_399')]],
                                                         "^mt-",col.name="pct.mito",assay = "Spatial")
  noMALDI.list[[paste0(x,'1','_399')]] <-
    noMALDI.list[[paste0(x,'1','_399')]][,noMALDI.list[[paste0(x,'1','_399')]]$nFeature_Spatial>300 
                                                   & noMALDI.list[[paste0(x,'1','_399')]]$nCount_Spatial>200]
  noMALDI.list[[paste0(x,'1','_399')]]$orig.ident<-paste0(x,'1','_399')
  print(VlnPlot(noMALDI.list[[paste0(x,'1','_399')]], 
                features = c("nFeature_Spatial", "nCount_Spatial", "pct.mito"), ncol = 3))
  print(SpatialFeaturePlot(noMALDI.list[[paste0(x,'1','_399')]], 
                           features = c("nFeature_Spatial", "nCount_Spatial", "pct.mito"), ncol = 3))
}

# postMALDI 359 data
for (x in c("A","B","C","D")){
  postMALDI.list[[paste0(x,'1','_359')]]<-PercentageFeatureSet(postMALDI.list[[paste0(x,'1','_359')]],
                                                         "^mt-",col.name="pct.mito",assay = "Spatial")
  postMALDI.list[[paste0(x,'1','_359')]] <-
    postMALDI.list[[paste0(x,'1','_359')]][,postMALDI.list[[paste0(x,'1','_359')]]$nFeature_Spatial>300 
                                                   & postMALDI.list[[paste0(x,'1','_359')]]$nCount_Spatial>200]
  postMALDI.list[[paste0(x,'1','_359')]]$orig.ident<-paste0(x,'1','_359')
  print(VlnPlot(postMALDI.list[[paste0(x,'1','_359')]], 
                features = c("nFeature_Spatial", "nCount_Spatial", "pct.mito"), ncol = 3))
  print(SpatialFeaturePlot(postMALDI.list[[paste0(x,'1','_359')]], 
                           features = c("nFeature_Spatial", "nCount_Spatial", "pct.mito"), ncol = 3))
}

```

# Using SCTransform to normalize the data. Then run dimensionality reduction, finding neighbors, and clustering, and rename cells
Note that SCTransform() replaces NormalizeData(), ScaleData() and FindVariableFeatures()
https://satijalab.org/seurat/articles/sctransform_vignette

```{r}
# wrapper function for normalization, reducing dimensions, finding neighbors and clustering.
# recommended that values for umap_dim, neighbor_dim, and k.param match.
NormSTNeighClust<-function(srt_obj, npcs=30, umap_dim=30, neighbor_dim=30, k.param = 30, clust.res = 0.8){
  srt_obj <- SCTransform(srt_obj, variable.features.n = 3000, assay="Spatial", verbose = F, vst.flavor="v2") 
  srt_obj <- RunPCA(srt_obj, assay ="SCT", npcs=npcs, verbose = F)
  srt_obj <- FindNeighbors(srt_obj, reduction = "pca", dims = 1:neighbor_dim, k.param=k.param, verbose = F) 
  srt_obj <- FindClusters(srt_obj, verbose = F, resolution = clust.res)
  srt_obj<-RunUMAP(srt_obj, dims=1:umap_dim, verbose = F)
  return(srt_obj)
}
# noMALDI 399 data
for (x in c('A', 'B', 'C', 'D')) {
  noMALDI.list[[paste0(x,'1','_399')]] <- NormSTNeighClust(noMALDI.list[[paste0(x,'1','_399')]])
}
# postMALDI 359 data
for (x in c('A', 'B', 'C', 'D')) {
  postMALDI.list[[paste0(x,'1','_359')]] <- NormSTNeighClust(postMALDI.list[[paste0(x,'1','_359')]])
}

# rename cells to ensure that we can identify origin when merging noMALDI and postMALDI data.
for (x in c('A', 'B', 'C', 'D')) {
  noMALDI.list[[paste0(x,'1','_399')]] <- RenameCells(noMALDI.list[[paste0(x,'1','_399')]], 
                                                      add.cell.id = paste0(x,'1','_399'))
  postMALDI.list[[paste0(x,'1','_359')]] <- RenameCells(postMALDI.list[[paste0(x,'1','_359')]], 
                                                      add.cell.id = paste0(x,'1','_359'))
}

```

# Viewing initial clustering UMAP and Spatial Clusters. 

```{r}
# Code for changing cluster resolutions (this was copied up here since more cluster resolutions are calculated later) 
# noMALDI 399
for (x in c('A', 'B', 'C', 'D')) {
  Idents(noMALDI.list[[paste0(x,'1','_399')]]) <- "SCT_snn_res.0.8"
  noMALDI.list[[paste0(x,'1','_399')]]$seurat_clusters <- noMALDI.list[[paste0(x,'1','_399')]]$SCT_snn_res.0.8
  #since some plots are grouping by seurat_clusters set this too
}
# postMALDI 359 
for (x in c('A', 'B', 'C', 'D')) {
  Idents(postMALDI.list[[paste0(x,'1','_359')]]) <- "SCT_snn_res.0.8"
  postMALDI.list[[paste0(x,'1','_359')]]$seurat_clusters <- postMALDI.list[[paste0(x,'1','_359')]]$SCT_snn_res.0.8
  #since some plots are grouping by seurat_clusters set this too
}

# noMALDI 399
for (x in c("A","B","C","D")){
  p1 <-DimPlot(noMALDI.list[[paste0(x,'1','_399')]], reduction = "umap", label = T)
  p2 <- SpatialDimPlot(noMALDI.list[[paste0(x,'1','_399')]], label = T, label.size = 3)
  print(p1+p2+ plot_annotation(title = paste0('399-',x,'1')))
}

# postMALDI 359
for (x in c("A","B","C","D")){
  p1<-DimPlot(postMALDI.list[[paste0(x,'1','_359')]], reduction = "umap", label = T)
  p2<-SpatialDimPlot(postMALDI.list[[paste0(x,'1','_359')]], label = T, label.size = 3)
  print(p1+p2+plot_annotation(title = paste0('359-',x,'1')))
}

```

# Exploring clusters with different resolutions. 
It is safe to run FindClusters multiple times on the same seurat object.
This chunk evaluates multiple cluster resolution possibilities.
```{r find.cluster.resolutions}

# note that resolution = 0.8 has already been run before.
# noMALDI 399
for (x in c('A', 'B', 'C', 'D')) {
  noMALDI.list[[paste0(x,'1','_399')]] <- FindClusters(noMALDI.list[[paste0(x,'1','_399')]],
                                                             verbose = F, resolution = c(0.4, 0.6, 1.0, 1.4, 1.6, 1.8, 2.0))
}
# postMALDI 359 
for (x in c('A', 'B', 'C', 'D')) {
  postMALDI.list[[paste0(x,'1','_359')]] <- FindClusters(postMALDI.list[[paste0(x,'1','_359')]],
                                                             verbose = F, resolution = c(0.4, 0.6, 1.0, 1.4, 1.6, 1.8, 2.0))
}

```

## Running FindAllMarkers on each slide
Only run it once, don't rerun since it takes a while to rerun.
Generated and stored for future use, if desired. 
``` {r}
# MAST needs to be set to spatial assay otherwise it fails
# this was initially calculated based on 0.8 resolution.
for (x in c('A', 'B', 'C', 'D')) {
  Idents(noMALDI.list[[paste0(x,'1','_399')]]) <- "SCT_snn_res.0.8"
}
noMALDI.allmarkers.cres0.8.list <- list()
for (x in c("A","B","C","D")){
  noMALDI.allmarkers.cres0.8.list[[paste0(x,'1','_399')]] <- FindAllMarkers(noMALDI.list[[paste0(x,'1','_399')]], 
                                                                    test.use = "MAST", assay = "Spatial")
}

# resolution 1.4
for (x in c('A', 'B', 'C', 'D')) {
  Idents(noMALDI.list[[paste0(x,'1','_399')]]) <- "SCT_snn_res.1.4"
}
noMALDI.allmarkers.cres1.4.list <- list()
for (x in c("A","B","C","D")){
  noMALDI.allmarkers.cres1.4.list[[paste0(x,'1','_399')]] <- FindAllMarkers(noMALDI.list[[paste0(x,'1','_399')]], 
                                                                    test.use = "MAST", assay = "Spatial")
}

# resolution 1
for (x in c('A', 'B', 'C', 'D')) {
  Idents(noMALDI.list[[paste0(x,'1','_399')]]) <- "SCT_snn_res.1"
}
noMALDI.allmarkers.cres1.list <- list()
for (x in c("A","B","C","D")){
  noMALDI.allmarkers.cres1.list[[paste0(x,'1','_399')]] <- FindAllMarkers(noMALDI.list[[paste0(x,'1','_399')]], 
                                                                    test.use = "MAST", assay = "Spatial")
}

# resolution 0.6
for (x in c('A', 'B', 'C', 'D')) {
  Idents(noMALDI.list[[paste0(x,'1','_399')]]) <- "SCT_snn_res.0.6"
}
noMALDI.allmarkers.cres0.6.list <- list()
for (x in c("A","B","C","D")){
  noMALDI.allmarkers.cres0.6.list[[paste0(x,'1','_399')]] <- FindAllMarkers(noMALDI.list[[paste0(x,'1','_399')]], 
                                                                    test.use = "MAST", assay = "Spatial")
}

# resolution 0.4
for (x in c('A', 'B', 'C', 'D')) {
  Idents(noMALDI.list[[paste0(x,'1','_399')]]) <- "SCT_snn_res.0.4"
}
noMALDI.allmarkers.cres0.4.list <- list()
for (x in c("A","B","C","D")){
  noMALDI.allmarkers.cres0.4.list[[paste0(x,'1','_399')]] <- FindAllMarkers(noMALDI.list[[paste0(x,'1','_399')]], 
                                                                    test.use = "MAST", assay = "Spatial")
}

# resolution 1.6
for (x in c('A', 'B', 'C', 'D')) {
  Idents(noMALDI.list[[paste0(x,'1','_399')]]) <- "SCT_snn_res.1.6"
}
noMALDI.allmarkers.cres1.6.list <- list()
for (x in c("A","B","C","D")){
  noMALDI.allmarkers.cres1.6.list[[paste0(x,'1','_399')]] <- FindAllMarkers(noMALDI.list[[paste0(x,'1','_399')]], 
                                                                    test.use = "MAST", assay = "Spatial")
}

# resolution 1.8
for (x in c('A', 'B', 'C', 'D')) {
  Idents(noMALDI.list[[paste0(x,'1','_399')]]) <- "SCT_snn_res.1.8"
}
noMALDI.allmarkers.cres1.8.list <- list()
for (x in c("A","B","C","D")){
  noMALDI.allmarkers.cres1.8.list[[paste0(x,'1','_399')]] <- FindAllMarkers(noMALDI.list[[paste0(x,'1','_399')]], 
                                                                    test.use = "MAST", assay = "Spatial")
}

# resolution 2.0
for (x in c('A', 'B', 'C', 'D')) {
  Idents(noMALDI.list[[paste0(x,'1','_399')]]) <- "SCT_snn_res.2"
}
noMALDI.allmarkers.cres2.list <- list()
for (x in c("A","B","C","D")){
  noMALDI.allmarkers.cres2.list[[paste0(x,'1','_399')]] <- FindAllMarkers(noMALDI.list[[paste0(x,'1','_399')]], 
                                                                    test.use = "MAST", assay = "Spatial")
}


# postMALDI 359 
for (x in c('A', 'B', 'C', 'D')) {
  Idents(postMALDI.list[[paste0(x,'1','_359')]]) <- "SCT_snn_res.0.8"
}
postMALDI.allmarkers.cres0.8.list <- list()
for (x in c("A","B","C","D")){
  postMALDI.allmarkers.cres0.8.list[[paste0(x,'1','_359')]] <- FindAllMarkers(postMALDI.list[[paste0(x,'1','_359')]], 
                                                                    test.use = "MAST", assay = "Spatial")
}

# postMALDI 359 resolution 1.4
for (x in c('A', 'B', 'C', 'D')) {
  Idents(postMALDI.list[[paste0(x,'1','_359')]]) <- "SCT_snn_res.1.4"
}
postMALDI.allmarkers.cres1.4.list <- list()
for (x in c("A","B","C","D")){
  postMALDI.allmarkers.cres1.4.list[[paste0(x,'1','_359')]] <- FindAllMarkers(postMALDI.list[[paste0(x,'1','_359')]], 
                                                                    test.use = "MAST", assay = "Spatial")
}

```


## Identifying spatially variable marker genes independent of clustering
Finding top spatially variable features for all slides, without accounting for clusters.  

``` {r}
# find top spatially variable features independent of clustering. Used top 3000 features as consistent with prior processing.
for (x in c("A","B","C","D")){
  noMALDI.list[[paste0(x,'1','_399')]] <- FindSpatiallyVariableFeatures(noMALDI.list[[paste0(x,'1','_399')]],
                              assay = "SCT", features = VariableFeatures(noMALDI.list[[paste0(x,'1','_399')]])[1:3000],
                              selection.method = "moransi")
}

# generate plots for the top 6 spatially variable genes
noMALDI.top6.svarfeats.plots <- list()
for (x in c("A","B","C","D")){
  top.features <- head(SpatiallyVariableFeatures(noMALDI.list[[paste0(x,'1','_399')]],
                                                        selection.method = "moransi"), 6)
  noMALDI.top6.svarfeats.plots[[x]] <- SpatialFeaturePlot(noMALDI.list[[paste0(x,'1','_399')]], 
                                         features = top.features, ncol = 3, alpha = c(0.1,1)) + 
                                          plot_annotation(title = paste0('399-',x,'1 Top 6 Spatially Variable Genes'))
}

for (x in c("A","B","C","D")){print(noMALDI.top6.svarfeats.plots[[x]])}

#this redoes above for 359 slices.
# find top spatially variable features independent of clustering. Used top 3000 features as consistent with prior processing.
for (x in c("A","B","C","D")){
  postMALDI.list[[paste0(x,'1','_359')]] <- FindSpatiallyVariableFeatures(postMALDI.list[[paste0(x,'1','_359')]],
                              assay = "SCT", features = VariableFeatures(postMALDI.list[[paste0(x,'1','_359')]])[1:3000],
                              selection.method = "moransi")
}

# generate plots for the top 6 spatially variable genes
postMALDI.top6.svarfeats.plots <- list()
for (x in c("A","B","C","D")){
  top.features <- head(SpatiallyVariableFeatures(postMALDI.list[[paste0(x,'1','_359')]],
                                                        selection.method = "moransi"), 6)
  postMALDI.top6.svarfeats.plots[[x]] <- SpatialFeaturePlot(postMALDI.list[[paste0(x,'1','_359')]], 
                                         features = top.features, ncol = 3, alpha = c(0.1,1)) + 
                                          plot_annotation(title = paste0('359-',x,'1 Top 6 Spatially Variable Genes'))
}

for (x in c("A","B","C","D")){print(postMALDI.top6.svarfeats.plots[[x]])}

```

# ST Deconvolve Annotations
https://www.nature.com/articles/s41467-022-30033-z
https://github.com/JEFworks-Lab/STdeconvolve
## Adding st deconvolve data to the Seurat Object Metadata
``` {r store st deconvolve}

# Importing STDeconvolve data from offline. 
# STDeconvolve cannot be run on the cluster due to installation issues.
# STDecon.399 <- readRDS(paste0(data_path,"/STDeconvolve_399_20230829.Rds")) # commented since it has already been imported
# STDecon.359 <- readRDS(paste0(data_path,"/STDeconvolve_359_20230829.Rds")) # commented since it has already been imported.

# Some notes about the STDecon data structure. each STDecon list contains an object storing the results for that corresponding slice
# for instance, STDecon.399[[3]] contains the results for slice C1-399
# pos stores tissue coordinates (GetTissueCoordinates(srt_obj))
# counts stores some processed counts data (from as.matrix(srt_obj@assays$Spatial@counts))
# deconProp is the derived pixel cell-type proportions
# deconGexp is the derived cell-type gene expression profiles
# results contains beta and theta, which are identical to deconGexp and deconProp respectively.
# ldas is stored because its how the modeling was run, ignore it. enables you to select other models if desired without waiting for STDeconvolve to rerun.
# marker.genes contains the list of of marker genes in order of decreasing log2fc per cell topic (celltype)... its the result of a different chunk of code.


# Function for assigning STDeconvolve cell topic data into the corresponding Seurat object as metadata
assign_STDecon_topics <- function (seurat_obj, STDecon_obj, prefix) {
  temp<-STDecon_obj$results$theta
  colnames(temp)<-paste0(prefix, "_cell_topic_", colnames(temp))
  seurat_annot<-AddMetaData(seurat_obj, temp)
  return(seurat_annot)
}

for (x in c("A","B","C","D")){
  noMALDI.list[[paste0(x,'1','_399')]] <- assign_STDecon_topics(noMALDI.list[[paste0(x,'1','_399')]],
                                                                     STDecon.399[[paste0(x,'1','_399')]], paste0(x,'1','_399'))
}
for (x in c("A","B","C","D")){
  postMALDI.list[[paste0(x,'1','_359')]] <- assign_STDecon_topics(postMALDI.list[[paste0(x,'1','_359')]],
                                                                     STDecon.359[[paste0(x,'1','_359')]], paste0(x,'1','_359'))
}

```


## Plotting bar plots of composition of STDeconvolve topics by cluster
``` {r}
get.avg.topic.scores <- function(seurat_obj, cluster.res) {
  # Extract the cell-topic score columns
  prediction_columns <- grep("cell_topic_", colnames(seurat_obj@meta.data), value = TRUE)
  # Calculate the average score of each cell topic column for each cell cluster identity
  average_scores <- aggregate(seurat_obj@meta.data[prediction_columns], 
                              by = list(seurat_obj@meta.data[paste0("SCT_snn_res.",cluster.res)][,1]), 
                              FUN = mean) 
  average_scores <- average_scores %>% rename("Group.1" = paste0("ClusterID.", cluster.res))
  # The above provides the average cell topic scoring for each seurat cluster.  
  return(average_scores)
}

topic.scores.barplotter <- function(topic_score_obj) {
  temp.df<-pivot_longer(topic_score_obj, cols = 2:ncol(topic_score_obj),
                      names_to = "cell_topic",
                      values_to = "avg_score")
  return(ggplot(temp.df, aes_string(fill="cell_topic", y="avg_score", x=(colnames(temp.df)[1]))) + 
      geom_bar(position="stack", stat="identity"))
}


###############
# printing out average cell topic scoring
###############
avg.topic.scores.399 <- list()
avg.topic.scores.359 <- list()
for (x in c('A', 'B', 'C', 'D')) {
  avg.topic.scores.399[[paste0(x,'1','_399')]] <- get.avg.topic.scores(noMALDI.list[[paste0(x,'1','_399')]],1.4)
  avg.topic.scores.359[[paste0(x,'1','_359')]] <- get.avg.topic.scores(postMALDI.list[[paste0(x,'1','_359')]],1.4)
}

topic.vs.clusters.399 <- list()
topic.vs.clusters.359 <- list()
for (x in c('A', 'B', 'C', 'D')) {
  topic.vs.clusters.399[[paste0(x,'1','_399')]] <- topic.scores.barplotter(avg.topic.scores.399[[paste0(x,'1','_399')]]) +
    plot_annotation(title = paste0(x,'1','_399'))
  topic.vs.clusters.359[[paste0(x,'1','_359')]] <- topic.scores.barplotter(avg.topic.scores.359[[paste0(x,'1','_359')]]) +
    plot_annotation(title = paste0(x,'1','_359'))
}

for (x in 1:4) {
  print(topic.vs.clusters.399[[x]])
  print(topic.vs.clusters.359[[x]])
}

```


# Analyses on merged data, combining spots from all slices in 399.
## Below is merge; combining all the spots from the no-MALDI without any batch correction.

```{r}
# noMALDI merge
merged.399 <- merge(x = noMALDI.list[[1]], y = c(noMALDI.list[[2]],noMALDI.list[[3]], noMALDI.list[[4]]))
merged.399 <- SCTransform(merged.399, assay="Spatial", variable.features.n = 3000, verbose= F, vst.flavor="v2")
merged.399 <- RunPCA(merged.399, assay ="SCT", npcs=30, verbose = F)
merged.399 <- FindNeighbors(merged.399, reduction = "pca", dims = 1:30, k.param=30, verbose = F)
merged.399 <- FindClusters(merged.399, verbose = F, resolution = 0.8)
merged.399 <- RunUMAP(merged.399, dims=1:30, verbose = F)

DimPlot(merged.399, label = T, repel = T, label.size = 3, label.box = T, group.by=c("seurat_clusters","orig.ident"),ncol=1) +
  plot_annotation(title = "Merged 4 Slices (No Integration)")
SpatialDimPlot(merged.399, label = T, repel = T, label.size = 3, label.box = T, group.by=c("seurat_clusters"), ncol = 2) +
  plot_annotation(title = "Merged 4 Slices (No Integration)")

#only run below if needed, since it is resource intensive.
# merged.399.all.markers <- FindAllMarkers(merged.399, test.use = "MAST", assay = "Spatial")

```
It looks like the clusters between slice A and B are pretty decent; since they're sequential. However, slices C and D are more distinct. 

# Harmonizing no-MALDI spots
https://portals.broadinstitute.org/harmony/SeuratV3.html
https://portals.broadinstitute.org/harmony/articles/quickstart.html
https://hbctraining.github.io/scRNA-seq_online/lessons/06a_integration_harmony.html

```{r}
# check to make sure the merged dataset has original dataset ids, which we will use harmony to integrate out
# note that unlike with the seurat method, we have to use a single merged dataset with all the prior things run on it
# hence we use merged.399 to continue here.
harmony.399 <- RunHarmony(merged.399, group.by.vars = "orig.ident", plot_convergence = T)

DimPlot(harmony.399, reduction = "harmony", label = T, repel = T, label.size = 3, label.box = T, group.by=c("seurat_clusters","orig.ident"),ncol=1) # note that you can change reduction to "umap" after running umap with harmony reduction aka further down in this workflow
# VlnPlot(object = harmony.399, features = "harmony_1", group.by = "orig.ident", pt.size = .1)

# run harmony
harmony.399 <- harmony.399 %>% 
  RunUMAP(reduction = "harmony", dims = 1:30) %>% 
  FindNeighbors(reduction = "harmony", dims = 1:30) %>% 
  FindClusters(resolution = 0.8)

#View output
print(DimPlot(harmony.399, reduction = "umap", label = T, repel = T, label.size = 3, label.box = T, 
              group.by=c("seurat_clusters","orig.ident"),ncol=1) + plot_annotation(title = "Harmony Integrated 4 Slices"))
print(SpatialDimPlot(harmony.399, group.by = "seurat_clusters", label = T, repel = T, 
                     label.size = 3, label.box = T, ncol = 2) +
        plot_annotation(title = "Harmony Integrated 4 Slices"))

#only run below if needed, since it is resource intensive.
# harmony.399.all.markers <- FindAllMarkers(harmony.399, test.use = "MAST", assay = "Spatial") # MAST needs Spatial otherwise it fails

```


# Spot cell identity predictions using Ayushi's scRNA data on the merged and harmonized datasets.
Integration predictions are completed independently for the merged dataset and the one merged using harmony.
See: https://satijalab.org/seurat/articles/spatial_vignette#integration-with-single-cell-data
https://hbctraining.github.io/scRNA-seq_online/lessons/06a_integration_harmony.html

``` {r}
# reference sc data from ayushi.
ayushi_sc_ref <- readRDS("../data/plt2_rna_atac_seurat_rnafiltered_cellfiltered.rds")
Idents(ayushi_sc_ref) <- "cell.ident.rna"
ayushi_sc_ref <- ayushi_sc_ref %>% RunPCA(verbose = F) %>% RunUMAP(dims = 1:30)

Idents(ayushi_sc_ref) <- "cell.ident.rna"
print(DimPlot(ayushi_sc_ref, reduction = "umap") + plot_annotation(title = "Ayushi SC Data", subtitle = "Cell Annotation"))
Idents(ayushi_sc_ref) <- "infercnv_annotation"
print(DimPlot(ayushi_sc_ref, reduction = "umap") + plot_annotation(title = "Ayushi SC Data", subtitle = "inferCNV Annotation"))

# use cell.ident.rna.
cell.idents<-unique(ayushi_sc_ref$cell.ident.rna)
# integration onto the merged dataset
DefaultAssay(merged.399) <- "SCT"
anchors.merged.399 <- FindTransferAnchors(reference = ayushi_sc_ref, query = merged.399, normalization.method = "SCT")
pred.assay.merged.399 <-TransferData(anchorset = anchors.merged.399, refdata = ayushi_sc_ref$cell.ident.rna, 
                                     prediction.assay = TRUE, weight.reduction = merged.399[["pca"]], dims = 1:30)
merged.399[["pred.cell.ident"]] <- pred.assay.merged.399

DefaultAssay(merged.399) <- "pred.cell.ident"
merged.399.cell.ID.plots <- list()
merged.399.cell.ID.plots[["Epithelial-cancer"]] <- SpatialFeaturePlot(merged.399, features = "Epithelial-cancer", ncol = 2)
merged.399.cell.ID.plots[["Epithelial-cells"]] <- SpatialFeaturePlot(merged.399, features = "Epithelial-cells", ncol = 2)
merged.399.cell.ID.plots[["Stromal-cells"]] <- SpatialFeaturePlot(merged.399, features = "Stromal-cells", ncol = 2)
merged.399.cell.ID.plots[["Immune-cells"]] <- SpatialFeaturePlot(merged.399, features = "Immune-cells", ncol = 2)
merged.399.cell.ID.plots[["Endothelial-cells"]] <- SpatialFeaturePlot(merged.399, features = "Endothelial-cells", ncol = 2)
merged.399.cell.ID.plots[["Fibroblasts"]] <- SpatialFeaturePlot(merged.399, features = "Fibroblasts", ncol = 2)

merged.399.cell.ID.plots

# integration onto the harmonized dataset
DefaultAssay(harmony.399) <- "SCT"
anchors.harmony.399 <- FindTransferAnchors(reference = ayushi_sc_ref, query = harmony.399, normalization.method = "SCT")
pred.assay.harmony.399 <-TransferData(anchorset = anchors.harmony.399, refdata = ayushi_sc_ref$cell.ident.rna, 
                                     prediction.assay = TRUE, weight.reduction = harmony.399[["harmony"]], dims = 1:30)
harmony.399[["pred.cell.ident"]] <- pred.assay.harmony.399

DefaultAssay(harmony.399) <- "pred.cell.ident"
harmony.399.cell.ID.plots <- list()
harmony.399.cell.ID.plots[["Epithelial-cancer"]] <- SpatialFeaturePlot(harmony.399, features = "Epithelial-cancer", ncol = 2)
harmony.399.cell.ID.plots[["Epithelial-cells"]] <- SpatialFeaturePlot(harmony.399, features = "Epithelial-cells", ncol = 2)
harmony.399.cell.ID.plots[["Stromal-cells"]] <- SpatialFeaturePlot(harmony.399, features = "Stromal-cells", ncol = 2)
harmony.399.cell.ID.plots[["Immune-cells"]] <- SpatialFeaturePlot(harmony.399, features = "Immune-cells", ncol = 2)
harmony.399.cell.ID.plots[["Endothelial-cells"]] <- SpatialFeaturePlot(harmony.399, features = "Endothelial-cells", ncol = 2)
harmony.399.cell.ID.plots[["Fibroblasts"]] <- SpatialFeaturePlot(harmony.399, features = "Fibroblasts", ncol = 2)

harmony.399.cell.ID.plots

# doesn't look very good, seurat integration method was not great it seems.

```


# Merging and harmonizing 359 and 399 slices, individually between slices and running harmony between slices to see if they are comparable.
This section is basically examining for comparability between the 359 and 399 runs; aka are the ST data post MALDI really comparable. I do this by treating the 359 vs 399 as a different "batch" and using harmony to see if they integrate well. 
``` {r}
merged.by.slice <- list()
for (x in c('A', 'B', 'C', 'D')) {
  merged.by.slice[[paste0(x,'1')]]<-merge(x = noMALDI.list[[paste0(x,'1','_399')]], y = postMALDI.list[[paste0(x,'1','_359')]])
  merged.by.slice[[paste0(x,'1')]]<-SCTransform(merged.by.slice[[paste0(x,'1')]], assay="Spatial", 
                                                variable.features.n = 3000, verbose= F, vst.flavor="v2")
  merged.by.slice[[paste0(x,'1')]] <- RunPCA(merged.by.slice[[paste0(x,'1')]], assay ="SCT", npcs=30, verbose = F)
  merged.by.slice[[paste0(x,'1')]] <- FindNeighbors(merged.by.slice[[paste0(x,'1')]], reduction = "pca", dims = 1:30, 
                                                    k.param=30, verbose = F)
  merged.by.slice[[paste0(x,'1')]] <- FindClusters(merged.by.slice[[paste0(x,'1')]], verbose = F, resolution = 0.8)
  merged.by.slice[[paste0(x,'1')]] <- RunUMAP(merged.by.slice[[paste0(x,'1')]], dims=1:30, verbose = F)
}
# before running next lines, check the orig.ident

# plotting just the merged, not harmonized by slice data
for (x in c('A', 'B', 'C', 'D')) {
  print(DimPlot(merged.by.slice[[paste0(x,'1')]], label = T, repel = T, label.size = 3, label.box = T,
                group.by=c("seurat_clusters","orig.ident"),ncol=1) + plot_annotation(title = 
                                                                                       paste0("Merged Slices (No Integration) ", x,'1')))
  print(SpatialDimPlot(merged.by.slice[[paste0(x,'1')]], group.by = "seurat_clusters", label = T, repel = T, 
                     label.size = 3, label.box = T, ncol = 2) + plot_annotation(title = 
                                                                                  paste0("Merged Slices (No Integration) ", x,'1')))
}
# plots demonstrate that harmonizing is required.

# before running this, check the orig.ident. don't rerun if not needed.
harmonized.by.slice <- list()
for (x in c('A', 'B', 'C', 'D')) {
  harmonized.by.slice[[paste0(x,'1')]] <- RunHarmony(merged.by.slice[[paste0(x,'1')]], group.by.vars = "orig.ident", plot_convergence = T)
  harmonized.by.slice[[paste0(x,'1')]] <- harmonized.by.slice[[paste0(x,'1')]] %>% 
                                            RunUMAP(reduction = "harmony", dims = 1:30) %>% 
                                            FindNeighbors(reduction = "harmony", dims = 1:30) %>% 
                                            FindClusters(resolution = c(0.2,0.3,0.4,0.5,0.6,0.7,0.8))
}

#plotting just the harmonized by slice data

for (x in c('A', 'B', 'C', 'D')) {
  Idents(harmonized.by.slice[[paste0(x,'1')]]) <- "SCT_snn_res.0.8"
  harmonized.by.slice[[paste0(x,'1')]]$seurat_clusters <- harmonized.by.slice[[paste0(x,'1')]]$SCT_snn_res.0.8
  
}
for (x in c('A', 'B', 'C', 'D')) {
  print(DimPlot(harmonized.by.slice[[paste0(x,'1')]], reduction = "umap", label = T, repel = T, label.size = 3, label.box = T,
                group.by=c("seurat_clusters","orig.ident"),ncol=1) + plot_annotation(title = paste0("Harmony Slices ", x,'1')))
  print(SpatialDimPlot(harmonized.by.slice[[paste0(x,'1')]], group.by = "seurat_clusters", label = T, repel = T, 
                     label.size = 3, label.box = T, ncol = 2) + plot_annotation(title = paste0("Harmony Slices ", x,'1')))
  
}

```


## Merging all spots

```{r}
merged.allslice <- merge(x = noMALDI.list[[1]], y = c(noMALDI.list[[2]],noMALDI.list[[3]], noMALDI.list[[4]], postMALDI.list[[1]],
                                                      postMALDI.list[[2]], postMALDI.list[[3]], postMALDI.list[[4]]))
merged.allslice <- SCTransform(merged.allslice, assay="Spatial", variable.features.n = 3000, verbose= F, vst.flavor="v2")
merged.allslice <- RunPCA(merged.allslice, assay ="SCT", npcs=30, verbose = F)
merged.allslice <- FindNeighbors(merged.allslice, reduction = "pca", dims = 1:30, k.param=30, verbose = F)
merged.allslice <- FindClusters(merged.allslice, verbose = F, resolution = 0.8)
merged.allslice <- RunUMAP(merged.allslice, dims=1:30, verbose = F)

DimPlot(merged.allslice, label = T, repel = T, label.size = 3, label.box = T, group.by=c("seurat_clusters","orig.ident"),ncol=1) +
  plot_annotation(title = "Merged All Slices (No Integration)")
SpatialDimPlot(merged.allslice, label = T, repel = T, label.size = 3, label.box = T, group.by=c("seurat_clusters"), ncol = 2) +
  plot_annotation(title = "Merged All Slices (No Integration)")

# Don't run the below again if we don't need to because its resource intensive.
# merged.allslice.all.markers <- FindAllMarkers(merged.allslice, test.use = "MAST", assay = "Spatial")

```
It looks like the clustering between slice A and B is pretty close together; since they're sequential. However, slices C and D are more distinct. 

# Harmonizing all spots from all slices

```{r}
harmony.allslice <- RunHarmony(merged.allslice, group.by.vars = "orig.ident", plot_convergence = T)
DimPlot(harmony.allslice, reduction = "harmony", label = T, repel = T, label.size = 3, label.box = T, group.by=c("seurat_clusters","orig.ident"),ncol=1) # note that you can change reduction to "umap" after running umap with harmony reduction aka further down in this workflow

harmony.allslice <- harmony.allslice %>% 
  RunUMAP(reduction = "harmony", dims = 1:30) %>% 
  FindNeighbors(reduction = "harmony", dims = 1:30) %>% 
  FindClusters(resolution = c(0.2, 0.4, 0.6, 0.8))


#change idents to change plotted clusters
Idents(harmony.allslice) <- "SCT_snn_res.0.8"
harmony.allslice$seurat_clusters <- harmony.allslice$SCT_snn_res.0.8 #since later plots are group by seurat clusters, need to set this too

print(DimPlot(harmony.allslice, reduction = "umap", label = T, repel = T, label.size = 1.5, label.box = T,
              group.by=c("seurat_clusters","orig.ident"),ncol=1) + plot_annotation(title = "Harmony Integrated All Slices 0.8"))
print(SpatialDimPlot(harmony.allslice, group.by = "seurat_clusters", label = T, repel = T,
                     label.size = 3, label.box = T, ncol = 2) +
        plot_annotation(title = "Harmony Integrated All Slices"))
print(DimPlot(merged.allslice, label = T, repel = T, label.size = 3, label.box = T,
              group.by=c("seurat_clusters","orig.ident"),ncol=1) +
      plot_annotation(title = "Merged All Slices (No Integration)"))
print(SpatialDimPlot(merged.allslice, label = T, repel = T, label.size = 3, label.box = T,
                     group.by=c("seurat_clusters"), ncol = 2) +
        plot_annotation(title = "Merged All Slices (No Integration)"))

```


# Get cell type specific markers from Ayushi's single cell data, and manually annotate seurat clusters
``` {r}
# FindAllMarkers on ayushi's scRNA data to get marker genes for her annotated clusters. Only run the FindAllMarkers code once--takes a long time!
# ayushi_sc_ref_allmarkers <- FindAllMarkers(ayushi_sc_ref, test.use = "MAST", assay = "RNA", only.pos = T)
ayushi_sc_ref_top10markers <- ayushi_sc_ref_allmarkers %>% group_by(cluster) %>% top_n(10, avg_log2FC)

# from here, I viewed the corresponding marker genes in the allmarkers set or top10 markers set, and manually cross-compared to viewing Seurat's cluster marker genes to manually assign annotations for each cluster type.
# I also used marker genes from the papers and cross referenced to seurat cluster marker genes list to annotate those clusters.
# papers: Travaglini, K.J., Nabhan, A.N., Penland, L. et al. Nature 587, 619–625 (2020); and Kim, N., Kim, H.K., Lee, K. et al. Nat Commun 11, 2285 (2020)

B1.399.top.genes.cres1.4 <- noMALDI.allmarkers.cres1.4.list[[2]] %>% group_by(cluster) %>% top_n(10, avg_log2FC)
View(B1.399.top.genes.cres1.4)
View(noMALDI.allmarkers.cres1.4.list[[2]])
View(ayushi_sc_ref_allmarkers)

```




# Examining markers and module expression of lung cells in all slices
```{r}
assignLungCellScores <- function(spatial_seurat) {
  # Markers below from: https://panglaodb.se/markers.html?cell_type=%27choose%27
  # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8996976/ important cell types
  # this applies module scores for each proposed cell type using each set of marker genes 
  spatial_seurat <- AddModuleScore(spatial_seurat, features = list(markers = c("Cd74","Cd79a","Cd79b","Ighm")), 
                                   name = "B-Cells")
  spatial_seurat <- AddModuleScore(spatial_seurat, features = list(markers = c("Cd3d","Cd3g","Cd3e","Ltb","Il7r")), 
                                   name = "T-Cells")
  spatial_seurat <- AddModuleScore(spatial_seurat, features = list(markers = c("Cd68","Fcgr1","Lyz2","Clec4a2",
                                                                               "Cxcl16","C3ar1")), 
                           name = "Macrophages")
  
  spatial_seurat <- AddModuleScore(spatial_seurat, features = list(markers = c("Siglecf","Marco","Mrc1","Itgax",
                                                                               "Gngt2", "Ear1")), 
                           name = "Alveolar-Macrophages")
  spatial_seurat <- AddModuleScore(spatial_seurat, features = list(markers = c("Cldn18","Ager","Col4a4","Cyp4b1",
                                                                               "Clic5","Col4a3","Pdpn",
                                                               "Krt7","Emp2","Aqp5","Hopx","Cyp2b10")), 
                           name = "PneumocyteT1-")
  spatial_seurat <- AddModuleScore(spatial_seurat, features = list(markers = 
                                              c("Nkx2-1","Etv5","Napsa","Lpcat1","Cd36","Lamp3","Cxcl15",
                                                "Lyz1","Cbr2","Ppp1r14c")), 
                           name = "PneumocyteT2-")
  
  spatial_seurat <- AddModuleScore(spatial_seurat, features = list(markers = c("Aqp1","Adh7","Cdh1","Sec14l3")), 
                           name = "AirEpithelial")
  spatial_seurat <- AddModuleScore(spatial_seurat, features = list(markers = c("Muc16","Dusp4","Muc4","Galnt6",
                                                               "Galnt5","Muc20","Scgb3a1")), 
                           name = "Goblet")
  spatial_seurat <- AddModuleScore(spatial_seurat, features = list(markers = c("Scgb1a1","Cyp2e1","Sftpd","Sftpa1",
                                                               "Aldh1a1","Cyp4b1","Muc1",
                                                               "Bpifa1","Ern1","Lepr", "Mmp28","Ahr",
                                                               "Muc4","Scgb3a1","Sftpc",
                                                               "Scgb3a2","Cyp2f2","Foxm1")), 
                           name = "Clara")
  spatial_seurat <- AddModuleScore(spatial_seurat, features = list(markers = c("Foxj1","Tubb4b","Sec14l3")), 
                           name = "Ciliated")
  
  spatial_seurat <- AddModuleScore(spatial_seurat, features = list(markers = c("Nog","Acta2","Foxf1","Gata5")), 
                           name = "airSmMuscle") # no canonical markers
  spatial_seurat <- AddModuleScore(spatial_seurat, features = list(markers = c("Col6a2","Vim","Pdgfrb","Vtn","Lum",
                                                               "Col1A2","Postn","Aspn","Prrx1",
                                                               "Col1a1","Pdgfra","Cela1","Lox","P4ha1",
                                                               "Ccr2","Prkcd","Cthrc1")), 
                           name = "Fibroblast")
  return(spatial_seurat)
}

plotLungCellDistr <- function(spatial_scored_seurat, plottitle = "title") {
  print(SpatialFeaturePlot(spatial_scored_seurat, features=c("B-Cells1","T-Cells1","Macrophages1"))
                       + plot_annotation(plottitle, subtitle = "Immune Cell Infiltration"))
  print(SpatialFeaturePlot(spatial_scored_seurat, features=c("Alveolar-Macrophages1","PneumocyteT1-1","PneumocyteT2-1"))
                       + plot_annotation(plottitle, subtitle = "Lung Cell Distribution"))
  print(SpatialFeaturePlot(spatial_scored_seurat, features=c("AirEpithelial1","Goblet1","Ciliated1"))
                       + plot_annotation(plottitle, subtitle = "Other Lung Cell Distribution")) 
  print(SpatialFeaturePlot(spatial_scored_seurat, features=c("Clara1", "PneumocyteT2-1"))
                       + plot_annotation(plottitle, subtitle = "Possible Tumor Cell Origin Distribution"))
  print(SpatialFeaturePlot(spatial_scored_seurat, features=c("airSmMuscle1","Fibroblast1"))
                       + plot_annotation(plottitle, subtitle = "Other Cell Distribution"))
}

for (x in c('A', 'B', 'C', 'D')) {
  noMALDI.list[[paste0(x,'1','_399')]] <- assignLungCellScores(noMALDI.list[[paste0(x,'1','_399')]])
  postMALDI.list[[paste0(x,'1','_359')]] <- assignLungCellScores(postMALDI.list[[paste0(x,'1','_359')]])
  plotLungCellDistr(noMALDI.list[[paste0(x,'1','_399')]], plottitle = paste0(x,'1','_399'))
  plotLungCellDistr(postMALDI.list[[paste0(x,'1','_359')]], plottitle = paste0(x,'1','_359'))
}

```
