suppressPackageStartupMessages({
  library(Seurat)
  library(sctransform)
  library(stringr)
  library(ggplot2)
  library(ggpubr)
  library(future)
  require(scales)
  library(RColorBrewer)
  library("readxl")
  library(dplyr)
  library(dendextend)
  library('conover.test')
  library("devtools")
  library(ggpmisc)
  library(VennDiagram)
  library(randomcoloR)
})


### figure 7 + supplements 

# import the Dopaminergic nuclei (mDA) dataset.
sobj <- readRDS("/path/to/dir/mDA.rds")

# remove Louvain resolutions from the merged object. 
sobj$seurat_clusters <- NULL

to.remove <- sapply(grep('snn', colnames(sobj@meta.data), value = T),
                    function(x) c(paste(x, collapse = ",")))

for(i in to.remove) {
  sobj[[i]] <- NULL
}

# split the dataset into a list of two seurat objects based on condition 
list <- SplitObject(sobj, split.by = "condition")

# normalize and identify variable features for each dataset independently
list <- lapply(X = list, FUN = SCTransform, variable.features.n=1000)

features <- SelectIntegrationFeatures(object.list = list, nfeatures = 1000)

list <- PrepSCTIntegration(object.list = list, anchor.features = features)

anchors <- FindIntegrationAnchors(object.list = list, normalization.method = "SCT", anchor.features = features)

integObj <- IntegrateData(anchorset = anchors, normalization.method = "SCT")

DefaultAssay(integObj) <- "integrated"

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

ElbowPlot(integObj, reduction = "pca", ndims = 100) + ggtitle("mDA_hvg1k_integrated")

# non-linear dim reduction:

DefaultAssay(integObj) <- "integrated"

integObj <- RunUMAP(integObj, reduction = "pca", dims = 1:30, verbose = TRUE)

integObj <- RunTSNE(integObj, reduction = "pca", dims = 1:30, verbose = TRUE)

## find neighbors and find clusters 
integObj <- FindNeighbors(integObj, reduction = "pca", dims = 1:30)

# graph clustering Seurat 
integObj  <- FindClusters(integObj, resolution = seq(1, 4.5, by = 0.5), n.start = 100, n.iter = 100)


## kmeans clustering 
# get the pca embeddings matrix of the desired pca range: 1:30

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

set.seed(69)
km71 <- kmeans(pcmat, 71, nstart = 100, iter.max = 1000, algorithm="MacQueen")

## add kmeans$clsuter to Seurat Obj metadata 
integObj@meta.data$kmeans71_integ <- km71$cluster 

## get the centroid of the kmeans for hierarchical clustering plus show labels. 
hc71 <- hclust(dist(km71$centers), method = "ward.D2")


# Manual Annotation of integrated territories 
# based on membership (Frequency) of integrated cluster to merged territories. 
# see supplementary file 7 sheet = mergedTERs_integCLUSTERs
# ambiguous clusters were resolved using integrated dendrogram and markers expression. 

anno.df <- as.data.frame(table(integObj$territory, integObj$kmeans71_integ))

# annotation file uploaded, added to integrated object. 
anno.df <- read_excel("/path/supplementary_file7/integrated_territories.xlsx", sheet = 2)
anno.df <- as.data.frame(anno.df)

IntegClusters <- as.character(anno.df$integratedClusters) 

IntegTerritories <- as.character(anno.df$`assigned territories`)

integObj$integTerritory <- plyr::mapvalues(
  x = integObj$kmeans71_integ, 
  from = IntegClusters,
  to = IntegTerritories)

saveRDS(integObj, "/path/to/dir/mDA_integrated.rds")

# fig7 A & B 
Idents(sobj) <- "condition"

tiff(file = "/path/to/dir/fig7A_B.tiff", 
     units="cm", width=40, height=20, res=300)

DimPlot(sobj, split.by = "condition", 
        cols = c("intact"="#08519C", "lesion"="#A63603")) + coord_fixed()

dev.off()


# fig7C merged territories in integrated object
TER_cols <- c("costumized_colors")

tiff(file = "/path/to/dir/fig7c.tiff", 
     units="cm", width=25, height=25, res=300)

DimPlot(integObj, group.by = "territory", cols = TER_cols, 
        order = c("ML", "Hy_DA"), pt.size = 0.7) + 
  theme(legend.text = element_text(size = 12, face = "bold")) + 
  coord_fixed() + ggtitle("merged Territory in integrated obj")

dev.off()

# fig7 D & E Sox6 & Calb1 in merged object split by condition:

# same scale in feature plot 
sapply(grep("Sox6",rownames(sobj@assays$RNA@data),value = TRUE),
       function(x) max(sobj@assays$RNA@data[x,]))

sapply(grep("Calb1",rownames(sobj@assays$RNA@data),value = TRUE),
       function(x) max(sobj@assays$RNA@data[x,]))

sapply(grep("Sox6",rownames(sobj@assays$RNA@data),value = TRUE),
       function(x) min(sobj@assays$RNA@data[x,]))

sapply(grep("Calb1",rownames(sobj@assays$RNA@data),value = TRUE),
       function(x) min(sobj@assays$RNA@data[x,]))

# fixed scale range based on the max & min values of the two genes from above. 

fix.sc <- scale_color_gradientn( colours = c('grey90', 'blue'),  limits = c(0, 4))

tiff(file = "/path/to/dir/fig7D_E.tiff", 
     units="cm", width=25, height=25, res=300)

p1 <- FeaturePlot(sobj, features = c('Sox6', 'Calb1'), 
                  split.by = 'condition', order = T, combine = FALSE) 

p2 <- lapply(p1, function (x) x + fix.sc) 

CombinePlots(p2) 

dev.off()

# Fig7 F highlight ML clusters from merged object in the integrated data
Idents(integObj) <- "territory"

ML_mergedData <- WhichCells(integObj, idents = 'ML')

tiff(file = "/path/to/dir/fig7F.tiff", 
     units="cm", width=25, height=25, res=300)

DimPlot(integObj, cells.highlight = ML_mergedData) + 
  scale_color_manual(labels = c("others", "ML_mergedData"), 
                     values = c("grey", "#08306B")) + 
  theme(legend.text = element_text(size = 10, face = "bold")) + 
  coord_fixed() + ggtitle("merged ML in integrated")

dev.off()


# fig supplement 6-2 A related to integrated object 
# plot UMAP Integrated cluster labels BUT with Territory colors, based on their membership 

clusterCols <- c("costumized_TERcolors")

p <- DimPlot(integObj, group.by = "kmeans71_integ", cols = clusterCols) + coord_fixed() + NoLegend()

tiff(file = "/path/to/dir/figSuppl6-2A.tiff", 
     units="cm", width=25, height=25, res=300)

LabelClusters(plot = p, id = 'kmeans71_integ', color = 'white', size=4, fontface = 'bold', 
              box = T, max.overlapp = Inf, repel=T) 

dev.off()

# fig supplement 6-2 C  VUL module in integrated clusters 
tiff(file = "/path/to/dir/figSuppl6-2C.tiff", 
     units="cm", width=40, height=10, res=300)

VlnPlot(integObj, features = "newVulmodule1", cols = clusterCols, 
        group.by = "kmeans71_integ", pt.size = 0, combine = TRUE) + 
  stat_summary(fun = mean, geom='point', size = 6, colour = "black", shape=95) + NoLegend() +
  labs(y="avg.Exp.gene.set", title = "vulnerable module integ clusters")

dev.off()


# fig supplement 6-2 D  RES module in integrated clusters 

tiff(file = "/path/to/dir/figSuppl6-2D.tiff", 
     units="cm", width=40, height=10, res=300)

VlnPlot(integObj, features = "newResmodule1", cols = clusterCols, 
        group.by = "kmeans71_integ", pt.size = 0, combine = TRUE) + 
  stat_summary(fun = mean, geom='point', size = 6, colour = "black", shape=95) + NoLegend() +
  labs(y="avg.Exp.gene.set", title = "resilient module integ clusters")

dev.off()

# fig supplement 7-1 B 
# plot ML-enriched markers across territories (supplementary file 8)

markers <- c("Atf3", "Creb5", "Xirp2", "Cd9", "Ecel1", "Clic4", "Sprr1a", 
             "Mmp12", "Hrk", "Akain1", "Cd44", "Ighm", "Pim1", "Tnfrsf12a", 
             "Nms", "Igf2bp2", "Arid5a", "Mustn1", "Cd109", "Lifr", "Qrfpr",
             "Pde8a", "Ell2",  "D5Ertd615e", "Bcl2l11", "Dusp10", "Ngf", 
             "Spata13", "Lamb3", "Fbln5")

tiff(file = "/path/to/dir/figSuppl7-1B.tiff", 
     units="cm", width=20, height=20, res=300)

DotPlot(sobj, assay = "RNA", features = markers, group.by = "territory") + coord_flip() + 
  theme(axis.text.x = element_text(size = 14, face = "bold", angle = 90, hjust = 1, vjust = 0.5), 
        axis.text.y = element_text(size = 14, face = "bold"))

dev.off()


# fig supplement 7-1 C
# visualize the overlap of gene sets with VennDiagram package

# import the markers:

# lesion vs intact DE 
lvi <- read.csv("/path/to/dir/lesion_intact_markers.csv", 
                        header = TRUE, sep = ",", quote = "\"")

# ML clusters markers 
ml_up <- read.csv("/path/to/dir/ML_up.csv", 
                  header = TRUE, sep = ",", quote = "\"")

ml_down <- read.csv("/path/to/dir/ML_down.csv", 
                    header = TRUE, sep = ",", quote = "\"")

# get the gene vectors:
lvi_up <- lvi[lvi$p_val_adj < 0.05 & lvi$avg_log2FC > 0, ]$X
lvi_down <- lvi[lvi$p_val_adj < 0.05 & lvi$avg_log2FC < 0, ]$X
ml_up <- ml_up[ml_up$p_val_adj < 0.05, ]$X
ml_down <- ml_down[ml_down$p_val_adj < 0.05, ]$X

# set the colors:
brewer.pal(3, 'Paired')
# "#A6CEE3" "#1F78B4" "#B2DF8A"
show_col(brewer.pal(3, 'Paired'))
myCol <- c("#A6CEE3", "#B2DF8A")

venn.diagram(
  x = list(ml_up, lvi_up),
  category.names = c("ML upregulated" , "lesioned upregulated"),
  filename = '/path/to/dir/figSupp7-1C-up.tiff',
  output=TRUE,
  imagetype="tiff" ,
  height = 700 , 
  width = 700 , 
  resolution = 2000,
  compression = "lzw",
  main = 'up-regulayed genes',
  main.fontface = 'bold', 
  main.cex = 0.05,
  lwd = 0,
  lty = 'blank', 
  fill = myCol,
  cex = 0.25,
  fontface = "bold",
  ext.text = TRUE,
  ext.line.lwd = 0.1, 
  ext.dist = -0.09,
  ext.length = 0.8,
  cat.cex = 0.05,
  cat.fontface = "bold", 
  cat.default.pos = "outer",
  cat.pos = c(0, 0),
  cat.dist = c(0.05, 0.2), 
  margin = 0.05)


venn.diagram(
  x = list(ml_down, lvi_down),
  category.names = c("ML down-regulated" , "lesion down-regulated"),
  filename = '/path/to/dir/figSupp7-1C-down.tiff',
  output=TRUE,
  imagetype="tiff" ,
  height = 700 , 
  width = 700 , 
  resolution = 2000,
  compression = "lzw",
  main = 'down-regulayed genes',
  main.fontface = 'bold', 
  main.cex = 0.05,
  lwd = 0,
  lty = 'blank', 
  fill = myCol,
  cex = 0.25,
  fontface = "bold",
  ext.text = 1,
  ext.line.lwd = 0.1, 
  ext.dist = 0.1,
  ext.length = 0.8,
  cat.cex = 0.05,
  cat.fontface = "bold", 
  cat.default.pos = "outer",
  cat.pos = c(0, 0),
  cat.dist = c(0.05, 0.2), 
  margin = 0.05)

# fig supplement 7-1 D & E 
# mostly_lesioned (ML) clusters enriched markers GSE analysis
# only genes unique to ML clusters and NOT the ones also enriched in lesion_intact DE: 

ml_up_unique <- read.csv("/path/to/dir/ML_up_unique.csv", 
                         header = TRUE, sep = ",", quote = "\"")

ml_down_unique <- read.csv("/path/to/dir/ml_down_unique.csv", 
                           header = TRUE, sep = ",", quote = "\"")

# get the gene names column:  
ml_up_unique <- ml_up_unique$x
ml_down_unique <- ml_down_unique$x

# next,  Enrichr for analysis 
install.packages("enrichR")

library(enrichR)

# list Enrichr sites:

listEnrichrSites()

dbs <- listEnrichrDbs()

websiteLive <- TRUE

if (is.null(dbs)) websiteLive <- FALSE
if (websiteLive) head(dbs)

# select the desired gene_set_libraries from the Enrichr databases:  

dbs <- c("GO_Molecular_Function_2021", "GO_Cellular_Component_2021", "GO_Biological_Process_2021")

if (websiteLive) {
  enriched_up <- enrichr(ml_up_unique, dbs)
}

if (websiteLive) {
  enriched_down <- enrichr(ml_down_unique, dbs)
}

# you may make and write the results tables.

Up_mol <- if (websiteLive) enriched_up[["GO_Molecular_Function_2021"]]

Up_cell <- if (websiteLive) enriched_up[["GO_Cellular_Component_2021"]]

Up_bio <- if (websiteLive) enriched_up[["GO_Biological_Process_2021"]]

write.csv(Up_mol, '/path/to/dir/up_GO_Molecular_Function_2021.csv')

write.csv(Up_cell, '/path/to/dir/up_GO_Cellular_Component_2021.csv')

write.csv(Up_bio, '/path/to/dir/up_GO_Biological_Process_2021.csv')

down_mol <- if (websiteLive) enriched_down[["GO_Molecular_Function_2021"]]

down_cell <- if (websiteLive) enriched_down[["GO_Cellular_Component_2021"]]

down_bio <- if (websiteLive) enriched_down[["GO_Biological_Process_2021"]]

write.csv(down_mol, '/path/to/dir/down_GO_Molecular_Function_2021.csv')
write.csv(down_cell, '/path/to/dir/down_GO_Cellular_Component_2021.csv')
write.csv(down_bio, '/path/to/dir/down_GO_Biological_Process_2021.csv')

## Plot Enrichr GO-BP output. (Plotting function contributed by I-Hsuan Lin)

p1 <- if (websiteLive) plotEnrich(enriched_up[['GO_Biological_Process_2021']], 
                                  showTerms = 20, numChar = 90, y = "Count", orderBy = "P.value",
                                  title = 'up-regulated GO_Biological_Process_2021')

p1 + theme(plot.title = element_text(size = 13, face = "bold"), 
           axis.text.x = element_text(size = 10, face = "bold"), 
           axis.text.y = element_text(size = 10, face = "bold"))


p2 <- if (websiteLive) plotEnrich(enriched_down[['GO_Biological_Process_2021']], 
                                  showTerms = 20, numChar = 90, y = "Count", orderBy = "P.value",
                                  title = 'down-regulated GO_Biological_Process_2021')

p2 + theme(plot.title = element_text(size = 13, face = "bold"), 
           axis.text.x = element_text(size = 10, face = "bold"), 
           axis.text.y = element_text(size = 10, face = "bold"))