# ========================================================-
# Title:         GSVA Analysis on scRNAseq and snRNAseq datasets
# Description:   Code to generate SuppFig13
# Author:        Monah Abou Alezz
# Date:          2025-03-06
# ========================================================-
suppressPackageStartupMessages(library(Seurat))
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(patchwork))
suppressPackageStartupMessages(library(pheatmap))
suppressPackageStartupMessages(library(GSVA))
suppressPackageStartupMessages(library(cogena))
suppressPackageStartupMessages(library(purrr))
suppressPackageStartupMessages(library(presto))
suppressPackageStartupMessages(library(DESeq2))
suppressPackageStartupMessages(library(org.Hs.eg.db))
suppressPackageStartupMessages(library(org.Mm.eg.db))
## load hallmark objects
gs_human <- gmt2list("data/h.all.v7.2.symbols.gmt")
gs_mouse <- gmt2list("data//h.all.v7.2.symbols_mouse.gmt")
## load datasets
mouse.combined <- readRDS("data/mouse_striatum_snRNA_AAV9.rds")
Idents(mouse.combined) <- "Cell_Category"
d4_final <- readRDS("data/mixed_d4_final.rds")
d4_aav9_ut <- subset(d4_final, orig.ident %in% c("AAV9", "UT"))
d4_aav9_ut@meta.data[['orig.ident']] <- recode_factor(d4_aav9_ut@meta.data[["orig.ident"]],
                                                      "AAV9" = "d4_AAV9",
                                                      "UT" = "d4_UT")
d4_aav9_ut@meta.data[["Cell_Category"]] <- recode_factor(
  d4_aav9_ut@meta.data[["Cell_Type"]], 
  "Undifferentiated" = "Undifferentiated", 
  "Oligodendrocytes" = "d4_Oligodendrocytes", 
  "OPC" = "d4_Oligodendrocytes", 
  "Astrocytes" = "d4_Astrocytes", 
  "Immature.Astrocytes" = "d4_Astrocytes", 
  "Neurons" = "d4_Neurons",
  "Immature.Neurons" = "d4_Neurons")
d4_aav9_ut$Cell_Category <- factor(d4_aav9_ut$Cell_Category,
                                   levels = c("d4_Astrocytes",
                                              "d4_Neurons",
                                              "d4_Oligodendrocytes",
                                              "Undifferentiated"))
d4_aav9_ut <- subset(d4_aav9_ut, Cell_Category != "Undifferentiated")
d4_aav9_ut$Cell_Category <- droplevels(d4_aav9_ut$Cell_Category)
Idents(d4_aav9_ut) <- "Cell_Category"
d5_final <- readRDS("data/organoids_d5_final.rds")
d5_aav9_ut <- subset(d5_final, orig.ident %in% c("AAV9", "UT"))
d5_aav9_ut@meta.data[['orig.ident']] <- recode_factor(d5_aav9_ut@meta.data[["orig.ident"]],
                                                      "AAV9" = "d5_AAV9",
                                                      "UT" = "d5_UT")
d5_aav9_ut@meta.data[["Cell_Category"]] <- recode_factor(
  d5_aav9_ut@meta.data[["Cell_Type"]], 
  "Undifferentiated" = "Undifferentiated", 
  "Oligodendrocytes" = "d5_Oligodendrocytes", 
  "OPC" = "d5_Oligodendrocytes", 
  "Astrocytes" = "d5_Astrocytes", 
  "Immature.Astrocytes" = "d5_Astrocytes", 
  "Neurons" = "d5_Neurons",
  "Immature.Neurons" = "d5_Neurons")
d5_aav9_ut$Cell_Category <- factor(d5_aav9_ut$Cell_Category,
                                   levels = c("d5_Astrocytes",
                                              "d5_Neurons",
                                              "d5_Oligodendrocytes",
                                              "Undifferentiated"))
d5_aav9_ut <- subset(d5_aav9_ut, Cell_Category != "Undifferentiated")
d5_aav9_ut$Cell_Category <- droplevels(d5_aav9_ut$Cell_Category)
Idents(d5_aav9_ut) <- "Cell_Category"
## create a pseudobulk and normalize
datasets_pseudo <- c("mouse.combined", "d4_aav9_ut", "d5_aav9_ut")
for (i in datasets_pseudo) {
  obj_i <- get(i)
  data_collapsed <- presto::collapse_counts(obj_i@assays$RNA@counts, 
                                            obj_i@meta.data, 
                                            c('Cell_Category', 'orig.ident'))
  meta_data<- data_collapsed$meta_data
  mat <- data_collapsed$counts_mat
  colnames(mat)<- paste(meta_data$Cell_Category, meta_data$orig.ident, sep="_")
  rownames(meta_data) <- paste(meta_data$Cell_Category, meta_data$orig.ident, sep="_")
  pseudo_bulk_obj <- list(mat = mat, meta_data = meta_data)
  assign(paste0(i, "_mat"), mat)
  assign(paste0(i,"_pseduo_bulk_obj"), pseudo_bulk_obj)
  DESeq.ds <- DESeq2::DESeqDataSetFromMatrix(countData = pseudo_bulk_obj$mat,
                                             colData = pseudo_bulk_obj$meta_data,
                                             design = ~orig.ident)
  DESeq.ds <- DESeq.ds[rowSums(counts(DESeq.ds)) > 0, ]
  DESeq.ds <- DESeq(DESeq.ds)
  DESeq.rlog <- rlogTransformation(DESeq.ds, blind = TRUE)
  counts <- assay(DESeq.rlog)
  assign(paste0(i,"_counts"), counts)
  if (i != "mouse.combined") {
    counts_es <- ExpressionSet(assayData = as.matrix(counts), annotation = org.Hs.eg.db)
  } else {
    counts_es <- ExpressionSet(assayData = as.matrix(counts), annotation = org.Mm.eg.db)
  }
  assign(paste0(i,"_es"), counts_es)
}
## using matrices
gsva_d4 <- gsva(as.matrix(d4_aav9_ut_counts), gs_human, min.sz=5, max.sz=500)
gsva_d5 <- gsva(as.matrix(d5_aav9_ut_counts), gs_human, min.sz=5, max.sz=500)
gsva_mouse <- gsva(as.matrix(mouse.combined_counts), gs_mouse, min.sz=5, max.sz=500)
gsva_counts_combined <- cbind(gsva_d4_1707_counts, gsva_d5_1743_counts, gsva_mouse.combined_counts)
p <- pheatmap(gsva_counts_combined)