suppressPackageStartupMessages({
  library(Seurat)
  library(stringr)
  library(sctransform)
  library(future)
  require(scales)
  library(RColorBrewer)
  library("readxl")
  library(dplyr)
  library(dendextend)
})


# After alignment with cellranger, (see cellrangercount.sh)
# the "filtered_feature_bc_matrix" per sample are read.

# List of the 6 intact (6 mo) samples and 3 untreated young mice (3 mo) and 3 untreated old mice (18 mo). 
# for mice sample IDs see supplementary file#1, sampleinfo_intact_wt. 

# Read sample info from CSV file
sample_info <- read.csv("/path/to/dir/sampleinfo_intact_wt.csv")

# List of sample IDs
sample_ids <- sample_info$`10X_Serial_ID`

base_dir <- "/common/path/string/to/cellranger/output/files/"

# List to store Seurat objects
seurat_objs <- list()

# Iterate over each sample ID
for (sample_id in sample_ids) {
  # File path for filtered_feature_bc_matrix
  file_path <- paste0(base_dir, sample_id, "/outs/filtered_feature_bc_matrix")
  
  # Read in the data
  data <- Read10X(file_path)
  
  # Create Seurat object
  seurat_obj <- CreateSeuratObject(counts = data, project = sample_id)
  
  # Add Seurat object to the list
  seurat_objs[[sample_id]] <- seurat_obj
}

# Merge all Seurat objects into one object
merged_obj <- merge(x = seurat_objs[[1]], y = seurat_objs[-1], 
                    add.cell.ids = sample_ids)

# Calculate percentage of mitochondrial genes
merged_obj <- PercentageFeatureSet(merged_obj, "^mt-", col.name = "percent_mito")

# Calculate percentage of ribosomal genes
merged_obj <- PercentageFeatureSet(merged_obj, "^Rp[sl]", col.name = "percent_ribo")

# Cell cycle scoring
merged_obj <- CellCycleScoring(merged_obj, g2m.features = str_to_title(cc.genes$g2m.genes), 
                               s.features = str_to_title(cc.genes$s.genes))

# Define a function to assign batch based on condition
assign_batch <- function(condition) {
  return(condition)
}

# Add batch to metadata using sample info
merged_obj$batch <- sapply(merged_obj$orig.ident, function(sample_id) {
  condition <- sample_info[sample_info$sample_ID == sample_id, "condition"]
  if (length(condition) > 0) {
    return(condition)
  } else {
    return("unknown")  # or any default value you prefer
  }
})



# Define a function to assign age based on age
assign_age <- function(age) {
  return(age)
}

# Add age to metadata using sample info
merged_obj$age <- sapply(merged_obj$orig.ident, function(sample_id) {
  age <- sample_info[sample_info$sample_ID == sample_id, "age"]
  if (length(age) > 0) {
    return(assign_age(age))
  } else {
    return("unknown")  # or any default value you prefer
  }
})


# Save merged Seurat object
saveRDS(merged_obj, file = "/path/to/dir/merged_seurat_obj.rds")

### cell / gene filtering 

selected_cells <- WhichCells(merged_obj, expression = nFeature_RNA > 500 & nFeature_RNA < 10000 & percent_mito < 5)

selected_genes <- rownames(merged_obj)[Matrix::rowSums(merged_obj) > 5]

sobj <- subset(merged_obj, features = selected_genes, cells = selected_cells)

# remove Malat1
sobj <- sobj[ ! grepl("Malat1", rownames(sobj)), ]

# dopaminergic nuclei filtering (mDA)
sobj <- subset(sobj, subset = Slc6a3 > 0 | Th > 0)

sobj <- SCTransform(sobj, variable.features.n = 1000)

sobj <- RunPCA(sobj, npcs = 100, verbose = FALSE)

ElbowPlot(sobj, reduction = "pca", ndims = 100)

sobj <- RunUMAP(sobj, dims = 1:30)

sobj <- RunTSNE(sobj, dims = 1:30)

DefaultAssay(sobj) <- "RNA"

sobj <- NormalizeData(sobj)

sobj <- FindVariableFeatures(sobj, selection.method = "vst", nfeatures = 1000)

all.genes <- rownames(sobj)

sobj <- ScaleData(sobj, features = all.genes)

#  PC 1:30   kmeans clustering  

pcmat <- Embeddings(sobj, reduction = 'pca')[,1:30]

set.seed(84)
km69s84 <- kmeans(pcmat, centers = 69, nstart = 50, iter.max = 1000, algorithm="MacQueen")

# add clusters to object
sobj@meta.data$km69s84 <- km69s84$cluster 

### identify and annotate clusters on this dataset based on markers' expression, dendrogram, etc, 

# create new metadata entry named "class" from annotated clusters 

sobj$class <- plyr::mapvalues(
  x = sobj$km69s84, 
  from = c('34', '9', '69', '6', '21', '13', '62', '39', '60', '7', '59', 
           '20', '63', '3', '12', '10', '19', '68', '1', '47', '31', '50', 
           '52', '23', '28', '53', '64', '65', '41', '2', '24', '27', '35', 
           '4', '42', '29', '44', '61', '55', '26', '37', '11', '8', '15', 
           '22', '32', '49', '14', '30', '51', '25', '46', '36', '67', '16', 
           '38', '56', '66', '18', '17', '40', '54', '58', '33', '5', '57', 
           '43', '45', '48'),
  to = rep(c('mDA', 'mODC', 'unassigned', 'Glut', 'Hy_DA', 'GABA', 'Glut'), 
           c(47, 1, 2, 4, 3, 10, 2)))

DimPlot(sobj, group.by = "class", cols = c("mODC" = "#B35806", "mDA" = "#2171B5", 
                                           "Glut" = "#DF65B0", "Hy_DA"="#00441B", 
                                           "GABA"="#FB9A99", "unassigned"="#969696"), 
        label = F, order = c('mODC', 'unassigned', 'GABA')) + coord_fixed() + 
  theme(text = element_text(size = 8, face = "bold"), legend.text=element_text(size = 14, face = 'bold'))

DimPlot(sobj, reduction = "umap", split.by = "batch", order = "untreated", 
        cols = c("intact"="#08519C", "untreated"="green")) + coord_fixed()

DefaultAssay(sobj) <- 'RNA'

markers <- c('Th', 'Slc6a3', 'Ddc', 'En1', 'Sox6', 'Calb1', 'Aldh1a1', 'Slc32a1', 'Gad2', 
             'Prlr', 'Satb2', 'Slc17a6', 'Nfib', 'Mog', 'Mag')

FeaturePlot(sobj, features = markers, slot = 'data', ncol = 5, coord.fixed = T)


sessionInfo()