5  Gene Ontology for each cluster

#pkg
library(Seurat)
library(SoupX)
library(scDblFinder)
library(Matrix)
library(rtracklayer)
library(GenomeInfoDb)
library(tidyverse)
library(forcats)
library(ggplot2)
library(here)
library(tidyverse)
library(readxl)
library(gt)
library(paletteer)
library(scales)
library(DT)
library(htmlwidgets)
library(rtracklayer)
library(GenomeInfoDb)
library(fs)
library(glmGamPoi)
library(furrr)
library(future)

if (requireNamespace("glmGamPoi", quietly = TRUE)) message("glmGamPoi available: SCT will be faster.")
if (requireNamespace("presto", quietly = TRUE)) message("presto available: marker finding will be faster.")

# src
source(here::here("src/function/stat_function/stat_analysis_main.R")) # for make plot 
source(here::here("src/function/fig_export.R")) # This function saves a given plot (plot_x) as both a PDF and a high-resolution PNG file at specified dimensions.

runs_base <- "/mnt/ibioldata/LBMC_data/cmaslard/data/ScBdLR_init_2025_private/cellranger/"
#runs_base <- "/Volumes/IBIOLdata_storage/LBMC_data/cmaslard/data/ScBdLR_init_2025_private/cellranger/"

# To del if needit #####
obj_dir <- here("data/output/snrnaseq_objects")

# Load merged object with SCT, PCA, UMAP, clusters, module scores
obj_all <- readRDS(file.path(obj_dir, "obj_all_SCT_PCA_UMAP_clusters.rds"))

DefaultAssay(obj_all) <- "SCT"
Idents(obj_all)       <- obj_all$seurat_clusters

markers <- readr::read_csv(file.path(obj_dir, "markers_res0.6.csv"))

# Recreate module-score based cell type guess if it is missing
if (is.null(obj_all$celltype_guess) || all(is.na(obj_all$celltype_guess))) {
  marker_sets <- list(
    Pericycle = c("BdiBd21-3.2G0092900.v1.2", "BdiBd21-3.1G0109300.v1.2"),
    Endodermis = c("BdiBd21-3.3G0792100.v1.2"),
    Cortex     = c("BdiBd21-3.4G0368600.v1.2")
  )

  DefaultAssay(obj_all) <- "RNA"
  ms_col_map <- character()
  for (nm in names(marker_sets)) {
    genes_ok <- intersect(marker_sets[[nm]], rownames(obj_all))
    if (length(genes_ok) < 5) next
    prev_cols <- colnames(obj_all@meta.data)
    obj_all <- AddModuleScore(
      obj_all,
      features = list(genes_ok),
      name     = paste0("MS_", nm),
      assay    = "RNA"
    )
    new_cols <- setdiff(colnames(obj_all@meta.data), prev_cols)
    if (length(new_cols)) ms_col_map[new_cols] <- nm
  }

  ms_cols <- names(ms_col_map)
  if (length(ms_cols)) {
    score_mat <- as.matrix(obj_all@meta.data[, ms_cols, drop = FALSE])
    score_mat[is.na(score_mat)] <- -Inf
    best_col <- ms_cols[max.col(score_mat, ties.method = "first")]
    obj_all$celltype_guess <- unname(ms_col_map[best_col])
  } else {
    obj_all$celltype_guess <- NA_character_
  }

  DefaultAssay(obj_all) <- "SCT"
}