vox1 <- PercentageFeatureSet(vox1, pattern = "^MT-", col.name = 'percent.mt')
vox1 <- subset(vox1, subset = 
                 nFeature_RNA > 200 & nFeature_RNA < 5000 & 
                 nCount_RNA > 500 & 
                 percent.mt < 11)
counts <- GetAssayData(vox1, assay = "RNA")
mt <- counts[-grep("^MT-", rownames(counts)),]
rp <- counts[-grep("^RP[SL]", rownames(counts)),]
vox1 <- subset(vox1, features = rownames(rp))

vox2 <- PercentageFeatureSet(vox2, pattern = "^MT-", col.name = 'percent.mt')
vox2 <- subset(vox2, subset = 
                 nFeature_RNA > 200 & nFeature_RNA < 5000 & 
                 nCount_RNA > 500 & 
                 percent.mt < 11)
counts <- GetAssayData(vox2, assay = "RNA")
mt <- counts[-grep("^MT-", rownames(counts)),]
rp <- counts[-grep("^RP[SL]", rownames(counts)),]
vox2 <- subset(vox2, features = rownames(rp))

vox3 <- PercentageFeatureSet(vox3, pattern = "^MT-", col.name = 'percent.mt')
vox3 <- subset(vox3, subset = 
                 nFeature_RNA > 200 & nFeature_RNA < 5000 & 
                 nCount_RNA > 500 &
                 percent.mt < 11)
counts <- GetAssayData(vox3, assay = "RNA")
mt <- counts[-grep("^MT-", rownames(counts)),]
rp <- counts[-grep("^RP[SL]", rownames(counts)),]
vox3 <- subset(vox3, features = rownames(rp))

vox4 <- PercentageFeatureSet(vox4, pattern = "^MT-", col.name = 'percent.mt')
vox4 <- subset(vox4, subset = 
                 nFeature_RNA > 200 & nFeature_RNA < 5000 & 
                 nCount_RNA > 500 & 
                 percent.mt < 11)
counts <- GetAssayData(vox4, assay = "RNA")
mt <- counts[-grep("^MT-", rownames(counts)),]
rp <- counts[-grep("^RP[SL]", rownames(counts)),]
vox4 <- subset(vox4, features = rownames(rp))

data <- merge(vox1, y = c(vox2, vox3, vox4), add.cell.ids = c("MTG", "HC", "PT", "WM"), project = "MitoBrainMap")
rm(vox1, vox2, vox3, vox4, dataloc, counts, mt, rp)

data$Brainregion <- NA
data$Brainregion[data$orig.ident == "voxel1"] <- 'Cortex'
data$Brainregion[data$orig.ident == "voxel2"] <- 'Hippocampus'
data$Brainregion[data$orig.ident == "voxel3"] <- 'Putamen'
data$Brainregion[data$orig.ident == "voxel4"] <- 'WhiteMatter'

data <- NormalizeData(data, normalization.method = "LogNormalize", scale.factor = 1e4)
data <- FindVariableFeatures(data, selection.method = 'vst', nfeatures = 3000)
data <- ScaleData(data) 

bcmvn %>%
  ggplot(aes(x = pK, y = MeanBC)) +
  geom_point() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggsave(here::here(Sys.Date(), "Plots", "DF_bcmvn.png"), width = 6, height =4.5)
nExp_poi <- round(0.075*nrow(data_df@meta.data))  ## Assuming 7.5% doublet formation rate - tailor for your dataset
data_df <- doubletFinder_v3(data_df, PCs = 1:pcs, pN = 0.25, pK = 5e-04, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)
colnames(data_df@meta.data) <- c("orig.ident",
                                      "nCount_RNA",
                                      "nFeature_RNA",
                                      "percent.mt",
                                      "Brainregion",
                                      "pANN",
                                      "DF_class")

tmp <- RunUMAP(data_df, dims = 1:pcs, seed.use = 4321)
DimPlot(tmp,label=F, label.size = 5, reduction = "umap",
        group.by ="DF_class", cols = c("coral3","darkcyan")) +
  guides(color = guide_legend(override.aes = list(size = 6))) +
  theme_umap

ggsave(here::here(Sys.Date(), "Plots", "UMAP2_DoubledFinder.png"), width = 5, height =4.5)
data_singlets <-data_df %>% subset(subset = DF_class == "Singlet")
rm(tmp)

set.seed(2)
data_harmony <- harmony::RunHarmony(data_singlets,"orig.ident", plot_convergence = TRUE, num.cores = detectCores() - 1)

rm(data, data_singlets)

data <- data_harmony
Idents(data) <- "Brainregion"

data <- data %>%
  RunUMAP(reduction = "harmony", dims = 1:pcs, verbose = F, seed.use = 123) %>% #123 #1
  FindNeighbors(reduction = "harmony", dims = 1:pcs) %>%
  FindClusters(resolution =0.2, random.seed = 2)

DimPlot(data,label=T, reduction = "umap")+
  NoLegend() +
  theme_umap

ggsave(here::here(Sys.Date(), "Plots", "UMAP3_clusters.png"), width = 5, height =4.5)

reference <- Azimuth:::LoadReference(path = "Azimuth", seconds = 10L)
query <- data
if (!all(c("nCount_RNA", "nFeature_RNA") %in% c(colnames(x = query[[]])))) {
    calcn <- as.data.frame(x = Seurat:::CalcN(object = query))
    colnames(x = calcn) <- paste(
      colnames(x = calcn),
      "RNA",
      sep = '_'
    )
    query <- AddMetaData(
      object = query,
      metadata = calcn
    )
    rm(calcn)
}

# Calculate percent mitochondrial genes if the query contains genes
# matching the regular expression "^MT-"
if (any(grepl(pattern = '^MT-', x = rownames(x = query)))) {
  query <- PercentageFeatureSet(
    object = query,
    pattern = '^MT-',
    col.name = 'percent.mt',
    assay = "RNA"
  )
}

# Filter cells based on the thresholds for nCount_RNA and nFeature_RNA
# you set in the app
cells.use <- query[["nCount_RNA", drop = TRUE]] <= 26235 &
  query[["nCount_RNA", drop = TRUE]] >= 485 &
  query[["nFeature_RNA", drop = TRUE]] <= 4984 &
  query[["nFeature_RNA", drop = TRUE]] >= 207

# If the query contains mitochondrial genes, filter cells based on the
# thresholds for percent.mt you set in the app
if ("percent.mt" %in% c(colnames(x = query[[]]))) {
  cells.use <- cells.use & (query[["percent.mt", drop = TRUE]] <= 12 &
    query[["percent.mt", drop = TRUE]] >= 0)
}

# Remove filtered cells from the query
query <- query[, cells.use]

# Preprocess with SCTransform
query <- SCTransform(
  object = query,
  assay = "RNA",
  new.assay.name = "refAssay",
  residual.features = rownames(x = reference$map),
  reference.SCT.model = reference$map[["refAssay"]]@SCTModel.list$refmodel,
  method = 'glmGamPoi',
  ncells = 2000,
  n_genes = 2000,
  do.correct.umi = FALSE,
  do.scale = FALSE,
  do.center = TRUE, verbose = F
)

# Find anchors between query and reference
anchors <- FindTransferAnchors(
  reference = reference$map,
  query = query,
  k.filter = NA,
  reference.neighbors = "refdr.annoy.neighbors",
  reference.assay = "refAssay",
  query.assay = "refAssay",
  reference.reduction = "refDR",
  normalization.method = "SCT",
  features = intersect(rownames(x = reference$map), VariableFeatures(object = query)),
  dims = 1:50,
  n.trees = 20,
  mapping.score.k = 100
)

# Transfer cell type labels and impute protein expression
#
# Transferred labels are in metadata columns named "predicted.*"
# The maximum prediction score is in a metadata column named "predicted.*.score"
# The prediction scores for each class are in an assay named "prediction.score.*"
# The imputed assay is named "impADT" if computed

refdata <- lapply(X = c("subclass", "class", "cluster"), function(x) {
  reference$map[[x, drop = TRUE]]
})
names(x = refdata) <- c("subclass", "class", "cluster")
if (FALSE) {
  refdata[["impADT"]] <- GetAssayData(
    object = reference$map[['ADT']],
    slot = 'data'
  )
}
query <- TransferData(
  reference = reference$map,
  query = query,
  dims = 1:50,
  anchorset = anchors,
  refdata = refdata,
  n.trees = 20,
  store.weights = TRUE, verbose = F
)

# Calculate the embeddings of the query data on the reference SPCA
query <- IntegrateEmbeddings(
  anchorset = anchors,
  reference = reference$map,
  query = query,
  reductions = "pcaproject",
  reuse.weights.matrix = TRUE,
  verbose = F
)

# Calculate the query neighbors in the reference
# with respect to the integrated embeddings
query[["query_ref.nn"]] <- FindNeighbors(
  object = Embeddings(reference$map[["refDR"]]),
  query = Embeddings(query[["integrated_dr"]]),
  return.neighbor = TRUE,
  l2.norm = TRUE
)

# The reference used in the app is downsampled compared to the reference on which
# the UMAP model was computed. This step, using the helper function NNTransform,
# corrects the Neighbors to account for the downsampling.
query <- Azimuth:::NNTransform(
  object = query,
  meta.data = reference$map[[]]
)

# Project the query to the reference UMAP.
query[["proj.umap"]] <- RunUMAP(
  object = query[["query_ref.nn"]],
  reduction.model = reference$map[["refUMAP"]],
  reduction.key = 'UMAP_',verbose = F
)


# Calculate mapping score and add to metadata
query <- AddMetaData(
  object = query,
  metadata = MappingScore(anchors = anchors),
  col.name = "mapping.score"
)

Idents(query) <- "orig.ident"
query_sub <- subset(query, idents = c("voxel4"), invert = T)

## Going back to data (without predictions but with voxel 4)
Idents(data) <- "seurat_clusters"
data <- subset(data, idents = c("0"), invert = T)
DimPlot(object = data, reduction = "umap", group.by = "seurat_clusters", label = TRUE) +
  labs(title = "UMAP after cleaning")

data$CellType <- NA
data$CellType[data$seurat_clusters == 0] <- "unspecific"
data$CellType[data$seurat_clusters == 1] <- "Oligodendrocytes"
data$CellType[data$seurat_clusters == 2] <- "Oligodendrocytes"
data$CellType[data$seurat_clusters == 3] <- "Excitatory"
data$CellType[data$seurat_clusters == 4] <- "Astrocytes"
data$CellType[data$seurat_clusters == 5] <- "Microglia"
data$CellType[data$seurat_clusters == 6] <- "OPC"
data$CellType[data$seurat_clusters == 7] <- "Endothelial"
data$CellType[data$seurat_clusters == 8] <- "Inhibitory"
data$CellType[data$seurat_clusters == 9] <- "Excitatory"
data$CellType[data$seurat_clusters == 10] <- "Pericytes"
data$CellType[data$seurat_clusters == 11] <- "VLMCs"
data$CellType_cluster <- paste0(data$CellType, '-',data$seurat_clusters)
Idents(data) <- "CellType"

DimPlot(data,label=T, reduction = "umap", group.by = "CellType",
        cols = unlist(col_cells))+
  NoLegend() +
  guides(color = guide_legend(override.aes = list(size = 6))) +
  theme_umap

data <- data %>%
  RunUMAP(reduction = "harmony", dims = 1:pcs, verbose = F, seed.use = 3) 

DimPlot(data,label=T, reduction = "umap", group.by = "CellType",
        cols = unlist(col_cells))+
  NoLegend() +
  guides(color = guide_legend(override.aes = list(size = 6))) +
  theme_umap

ggsave(here::here(Sys.Date(), "Plots", "Figure3B_UMAP_annotated.png"), width = 5, height =4.5)


p <- DimPlot(data,label=F, reduction = "umap", group.by = "CellType",
        cols = unlist(col_cells))+
  NoLegend() +
  guides(color = guide_legend(override.aes = list(size = 6))) +
  theme_umap
ggsave(here::here(Sys.Date(), "Plots", "Figure3B_UMAP_blank.png"),p, width = 5, height =4.5)

## Add cells per cluster annotation 
cells_per_cluster <- data@meta.data %>%
  filter(!seurat_clusters %in% 0) %>%
  group_by(Brainregion, CellType) %>%
  mutate(count = n()) %>%
  select(Brainregion, CellType, count) %>%
  unique()

cells_per_cluster  %>%
    mutate(Brainregion = case_when(
    Brainregion == "Cortex" ~ "Cort",
    Brainregion == "Hippocampus" ~ "Hipp",
    Brainregion == "Putamen" ~ "Put",
    Brainregion == "WhiteMatter" ~ "WM",
  )) %>%
  ggplot(aes(x =Brainregion,
             y= count, fill = fct_relevel(CellType, "Pericytes", "Endothelial", "VLMCs", "Microglia", 
                                          "OPC", "Oligodendrocytes", "Astrocytes", "Inhibitory", "Excitatory"))) +
  geom_bar(position = "fill", stat="identity",
           size = 0.2, color = "black") +
  scale_fill_manual(values = unlist(col_cells)) +
  labs(y="Fraction of total") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_blank(),
        axis.title.x = element_blank())

ggsave(here::here(Sys.Date(), "Plots", "Figure3C_Proportions.png"), width = 4, height =3)#width = 7.5, height = 6.3)

mito_data <- NormalizeData(data,normalization.method = "LogNormalize", scale.factor = 10000)
mito_data <- subset(mito_data, features = mitogenes)
pseudobulk_mito <- AverageExpression(mito_data, group.by = c("orig.ident"), layer = "data", return.seurat=F)
## Since slot="data", the data is assumed to be log transformed. Return.seurat = F, average non log values are returned 
pseudobulk_mito  <- pseudobulk_mito[["RNA"]] %>%
  as.data.frame() %>%
  rownames_to_column("Gene") %>%
  pivot_longer(cols = -Gene) %>%
  mutate(name = case_when(
    name == "voxel1" ~ "Cort",
    name == "voxel2" ~ "Hipp",
    name == "voxel3" ~ "Put",
    name == "voxel4" ~ "WM"
  )) %>%
  rename(Voxel = name) %>%
  group_by(Gene) %>%
  mutate(sum = sum(value)) %>%
  filter(!sum %in% 0) %>%
  select(-sum) %>%
  pivot_wider(names_from = "Gene", values_from = "value") 

pseudobulk_mitopathway_scores <- pseudobulk_mito %>%
  pivot_longer(cols = -c(Voxel), names_to = "Gene") %>%
  full_join(gene_to_pathway, by = "Gene") %>%
  filter(!is.na(Pathway)) %>%
  group_by(Voxel,Pathway) %>%
  mutate(score = mean(value, na.rm= T)) %>%
  select(-c(Gene, value)) %>%
  unique() %>%
  na.omit() 

pseudobulk_oxphos_scores <- pseudobulk_mitopathway_scores %>%
  filter(Pathway %in% oxphos_pathways) %>%
  mutate(Pathway =case_when(
    Pathway == "Complex I" ~ "CI exprs",
    Pathway == "Complex II" ~ "CII exprs",
    Pathway == "Complex III" ~ "CIII exprs",
    Pathway == "Complex IV" ~ "CIV exprs",
    Pathway == "Complex V" ~ "CV exprs",
    Pathway == "Cytochrome C"~"Cyt C",
    Pathway == "OXPHOS"~"OxPhos",
    Pathway == "OXPHOS assembly factors"~"OxPhos assembly factors",
    Pathway == "OXPHOS subunits"~"OxPhos subunits",
    TRUE ~ Pathway))  %>%
## Zscore
  group_by(Pathway) %>%
  mutate(mean = mean(score, na.rm = T)) %>%
  mutate(sd = sd(score, na.rm = T)) %>%
  mutate(score = (score - mean) / sd) %>%
  select(-c(mean, sd)) %>%
  pivot_wider(names_from = "Pathway", values_from = "score") 

pseudobulk_mito_content <- pseudobulk_mito %>%
  pivot_longer(cols = -Voxel) %>%
  group_by(Voxel) %>%
  mutate(estim_content = mean(value, na.rm = T)) %>%
  select(-c(name, value)) %>%
  unique() 

additional_data <- read.csv(here::here("MRC", "voxel_MRC_data.csv")) %>%
  select(Voxel, MRC, TRC, CI, CII, CIV, MitoD) %>%
   mutate(Voxel = case_when(
    Voxel == "voxel1" ~ "Cort",
    Voxel == "voxel2" ~ "Hipp",
    Voxel == "voxel3" ~ "Put",
    Voxel == "voxel4" ~ "WM"
  )) %>%
  full_join(pseudobulk_mito_content) %>%
  pivot_longer(cols = -Voxel) %>%
  mutate(name = case_when(
    name == "CI" ~ "CI-activity",
    name == "CII" ~ "CII-activity",
    name == "CIV" ~ "CIV-activity",
    TRUE ~ name
  )) %>%
  group_by(name) %>%
  mutate(mean = mean(value, na.rm = T)) %>%
  mutate(sd = sd(value, na.rm = T)) %>%
  mutate(value = (value - mean) / sd) %>%
  select(-c(mean, sd)) %>%
  pivot_wider(names_from = "name", values_from = "value") 

pseudobulk_oxphos_scores <- pseudobulk_oxphos_scores %>%
  full_join(additional_data) %>%
  select(Voxel, 
  OxPhos, `OxPhos assembly factors`, `OxPhos subunits`,
  `CI exprs`, `CI subunits`,`CI assembly factors`,`CI-activity`,
  `CII exprs`,`CII subunits`, `CII assembly factors`,`CII-activity`,
  `CIII exprs`,`CIII subunits`, `CIII assembly factors`,
  `CIV exprs`,`CIV subunits`, `CIV assembly factors`,`CIV-activity`,
  `CV exprs`,`CV subunits`, `CV assembly factors`, 
  `Cyt C`, estim_content, 
  MitoD, MRC, TRC  )



## Transpose data
exprs <- t(pseudobulk_oxphos_scores%>%column_to_rownames("Voxel"))

## Heatmap annotation
anno_df <- as.data.frame(pseudobulk_oxphos_scores) %>%
  pivot_longer(cols = -c(Voxel), names_to= "Gene") %>%
  na.omit() %>%
  pivot_wider(names_from = "Gene", values_from = "value") %>%
  dplyr::select(Voxel)
col_anno = columnAnnotation(
  Voxel = anno_df$Voxel, 
  col=list(
  Voxel = c(
    "Cort" = "#4F5088", 
    "Hipp" = "#A5458D", 
    "Put" = "#5D9E4F", 
    "WM" = "#4E77BC"
  )
  ),
  show_annotation_name = F,
  show_legend =  T,
  annotation_legend_param = list(nrow=12)) 

## Build the heatmap
hm <- Heatmap(exprs, 
        name = "Rel. expression",
        show_column_dend=F,
        show_column_names =T,
        show_row_names =T,
        row_title="",
        cluster_columns =F,
        cluster_rows =F,
        row_names_gp = grid::gpar(fontsize = 8),
        row_title_gp = grid::gpar(fontsize = 10),
        column_names_rot = 45,
        column_names_gp = grid::gpar(fontsize = 12),
        column_title_gp = grid::gpar(fontsize = 10),
        top_annotation=col_anno,
        width = unit(5, "cm")
           )
pdf(here::here(Sys.Date(), "Plots", "Figure3D_Heatmap_oxphos_pseudobulk_per_voxel.pdf"), width = 8, height = 4)
draw(hm)
invisible(dev.off())
draw(hm)

rm(pseudobulk_mito, pseudobulk_mitopathway_scores, pseudobulk_oxphos_scores, pseudobulk_mito_content, additional_data, exprs, anno_df, col_anno, hm)

mito_data <- NormalizeData(data,normalization.method = "LogNormalize", scale.factor = 10000)
mito_data <- subset(mito_data, features = mitogenes)
pseudobulk_mito_vox_cell <- AverageExpression(mito_data, group.by = c("orig.ident", "CellType"),
                                              slot = "data", return.seurat=F)

pseudobulk_mito_vox_cell  <-pseudobulk_mito_vox_cell[["RNA"]] %>%
  as.data.frame() %>%
  rownames_to_column("Gene") %>%
  pivot_longer(cols = -Gene) %>%
  separate(name, into = c("Voxel", "CellType"), remove = T) %>%
  mutate(Voxel = case_when(
    Voxel == "voxel1" ~ "Cort",
    Voxel == "voxel2" ~ "Hipp",
    Voxel == "voxel3" ~ "Put",
    Voxel == "voxel4" ~ "WM"
  )) %>%
  group_by(Gene) %>%
  mutate(sum = sum(value)) %>%
  filter(!sum %in% 0) %>%
  select(-sum) 

pseudobulk_mito_vox_cell_pathways = pseudobulk_mito_vox_cell %>%
  full_join(gene_to_pathway, by = "Gene") %>%
  filter(!is.na(Pathway)) %>%
  group_by(Voxel,CellType, Pathway) %>%
  mutate(score = mean(value, na.rm= T)) %>%
  na.omit() %>%
  select(-c(Gene, value)) %>%
  unique() 

mito_raw <- pseudobulk_mito_vox_cell_pathways %>%
  pivot_wider(names_from = "Pathway", values_from = "score") %>%
  unite(name, Voxel, CellType, remove =F) %>%
  column_to_rownames("name")

col_vox <- list(  
    "Cort" = "#4F5088", 
    "Hipp" = "#A5458D", 
    "Put" = "#5D9E4F", 
    "WM" = "#4E77BC")

theme <-   theme(
    legend.title = element_blank(),
    axis.title = element_text(size = 10), 
    axis.text = element_text(size =8)) 

p <- mito_raw %>%
  ggplot(aes(x = `Complex II`, y = `Complex I`)) +
  geom_point(aes(fill = Voxel, color = Voxel), alpha = 0.6, size = 4) +
  scale_color_manual(values = col_vox) +
  scale_fill_manual(values = col_vox) +
  labs(title = "CI + CII") + 
  theme_bw() +
  theme(legend.position = "none")
p

p <- mito_raw %>%
  ggplot(aes(x = `Complex III`, y = `Complex I`)) +
  geom_point(aes(fill = Voxel, color = Voxel), alpha = 0.6, size = 4) +
  scale_color_manual(values = col_vox) +
  scale_fill_manual(values = col_vox) +
  labs(title = "CI + CIII") + 
  theme_bw() +
  theme(legend.position = "none")
p

p <- mito_raw %>%
  ggplot(aes(x = `Complex IV`, y = `Complex I`)) +
  geom_point(aes(fill = Voxel, color = Voxel), alpha = 0.6, size = 4) +
  scale_color_manual(values = col_vox) +
  scale_fill_manual(values = col_vox) +
  labs(title = "CI + CIV (Figure 3e)") + 
  theme_bw() +
  theme(legend.position = "none")
p

ggsave(here::here(Sys.Date(), "Plots", "Figure3E_Bivariate_CI_CIV_Raw.png" ), height = 2.8, width = 2.55)

p <- mito_raw %>%
  ggplot(aes(x = `Complex V`, y = `Complex I`)) +
  geom_point(aes(fill = Voxel, color = Voxel), alpha = 0.6, size = 4) +
  scale_color_manual(values = col_vox) +
  scale_fill_manual(values = col_vox) +
  labs(title = "CI + CV") + 
  theme_bw() +
  theme(legend.position = "none")
p

pseudobulk_mito_vox_cell_pathways = pseudobulk_mito_vox_cell %>%
  full_join(gene_to_pathway, by = "Gene") %>%
  filter(!is.na(Pathway)) %>%
  group_by(Voxel,CellType, Pathway) %>%
  mutate(score = mean(value, na.rm= T)) %>%
  na.omit() %>%
  select(-c(Gene, value)) %>%
  unique() 

mito_raw <- pseudobulk_mito_vox_cell_pathways %>%
  pivot_wider(names_from = "Pathway", values_from = "score") %>%
  unite(name, Voxel, CellType, remove =F) %>%
  column_to_rownames("name")

## Transpose and scale matrix
exprs <- t(scale(mito_raw[,-c(1:2)]))

## Heatmap annotation
anno_df <- as.data.frame(mito_raw) %>%
  pivot_longer(cols = -c(Voxel, CellType), names_to= "Pathway") %>%
  na.omit() %>%
  pivot_wider(names_from = "Pathway", values_from = "value") %>%
  dplyr::select(Voxel, CellType)
col_anno = columnAnnotation(
  Celltype=anno_df$CellType,
  Voxel = anno_df$Voxel, 
  col=list(
  Celltype  =unlist(col_cells),
  Voxel = c(
    "Cort" = "#4F5088", 
    "Hipp" = "#A5458D", 
    "Put" = "#5D9E4F", 
    "WM" = "#4E77BC")
  ),
  show_annotation_name = F,
  show_legend =  T,
  annotation_legend_param = list(nrow=12)) 

col_dist    = stats::dist(t(as.matrix(exprs)), method="canberra")
coldend     = hclust(col_dist, method="ward.D2")
row_dist    = stats::dist(as.matrix(exprs), method="canberra")
rowdend     = hclust(row_dist, method="ward.D2")


hm <- Heatmap(exprs, 
        name = "Rel. expression",
        show_column_dend=T,
        show_column_names =T,
        show_row_names =F,
        show_row_dend=F,
        row_title="",
        cluster_columns =coldend,
        cluster_rows =rowdend,
        column_names_rot = 45,
        row_names_gp = grid::gpar(fontsize = 8),
        row_title_gp = grid::gpar(fontsize = 10),
        column_names_gp = grid::gpar(fontsize = 8),
        column_title_gp = grid::gpar(fontsize = 10),
        top_annotation=col_anno,
        width = unit(12, "cm"))
pdf(here::here(Sys.Date(), "Plots", "Figure3F_Heatmap_pathways_scaled.pdf"), width = 25, height = 5)
draw(hm)
invisible(dev.off())
Heatmap(exprs,
        name = "Rel. expression",
        show_column_dend=T,
        show_column_names =T,
        show_row_names =T,
        show_row_dend=T,
        row_title="Mitopathways",
        cluster_columns =coldend,
        cluster_rows =rowdend,
        column_names_rot = 45,
        row_names_gp = grid::gpar(fontsize = 2),
        row_title_gp = grid::gpar(fontsize = 10),
        column_names_gp = grid::gpar(fontsize = 6),
        column_title_gp = grid::gpar(fontsize = 10),
        top_annotation=col_anno,
        width = unit(8, "cm")
           )

rm(pseudobulk_mito_vox_cell_zscore, exprs, anno_df, col_anno, hm)

mito_zscore <- pseudobulk_mito_vox_cell_pathways %>%
  group_by(Voxel, Pathway) %>%
  mutate(mean = mean(score, na.rm = T)) %>%
  mutate(sd = sd(score, na.rm =T)) %>%
  mutate(score = (score-mean)/sd) %>%
  select(-c(mean, sd)) %>%
  pivot_wider(names_from = "Pathway", values_from = "score") %>%
  unite(name, Voxel, CellType, remove =F) %>%
  column_to_rownames("name")

p <- mito_zscore %>%
  ggplot(aes(x = `Complex II`, y = `Complex I`)) +
  geom_point(aes(fill = CellType, color = CellType), alpha = 0.6, size = 4) +
  scale_color_manual(values = col_cells) +
  scale_fill_manual(values = col_cells) +
  labs(title = "CI + CII") + 
  theme_bw() +
  theme(legend.position = "none")
p

p <- mito_zscore %>%
  ggplot(aes(x = `Complex III`, y = `Complex I`)) +
  geom_point(aes(fill = CellType, color = CellType), alpha = 0.6, size = 4) +
  scale_color_manual(values = col_cells) +
  scale_fill_manual(values = col_cells) +
  labs(title = "CI + CIII") + 
  theme_bw() +
  theme(legend.position = "none")
p

p <- mito_zscore %>%
  ggplot(aes(x = `Complex IV`, y = `Complex I`)) +
  geom_point(aes(fill = CellType, color = CellType), alpha = 0.6, size = 4) +
  scale_color_manual(values = col_cells) +
  scale_fill_manual(values = col_cells) +
  labs(title = "CI + CIV (Figure 3g)") + 
  theme_bw() +
  theme(legend.position = "none")
p

 ggsave(here::here(Sys.Date(), "Plots", "Figure3G_Bivariate_CI_CIV_Zscore.png" ), height = 2.8, width = 2.55)

p <- mito_zscore %>%
  ggplot(aes(x = `Complex V`, y = `Complex I`)) +
  geom_point(aes(fill = CellType, color = CellType), alpha = 0.6, size = 4) +
  scale_color_manual(values = col_cells) +
  scale_fill_manual(values = col_cells) +
  labs(title = "CI + CV") + 
  theme_bw() +
  theme(legend.position = "none")
p

pseudobulk_mito_vox_cell_zscore <- pseudobulk_mito_vox_cell %>%
  full_join(gene_to_pathway, by = "Gene") %>%
  filter(!is.na(Pathway)) %>%
  group_by(Voxel,CellType, Pathway) %>%
  mutate(score = mean(value, na.rm= T)) %>%
  select(-c(Gene, value)) %>%
  unique() %>%
  na.omit() %>%
  ## zscore
  group_by(Pathway, Voxel) %>%
  mutate(mean = mean(score, na.rm = T)) %>%
  mutate(sd = sd(score, na.rm =T)) %>%
  mutate(score = (score-mean)/sd) %>%
  select(-c(mean, sd)) %>%
  pivot_wider(names_from = "Pathway", values_from = "score") %>%
  unite(name, Voxel, CellType, remove =F) %>%
  column_to_rownames("name")

## Transpose matrix
exprs <- t(pseudobulk_mito_vox_cell_zscore[,-c(1:2)])

## Heatmap annotation
anno_df <- as.data.frame(pseudobulk_mito_vox_cell_zscore) %>%
  pivot_longer(cols = -c(Voxel, CellType), names_to= "Pathway") %>%
  na.omit() %>%
  pivot_wider(names_from = "Pathway", values_from = "value") %>%
  dplyr::select(Voxel, CellType)
col_anno = columnAnnotation(
  Celltype=anno_df$CellType,
  Voxel = anno_df$Voxel, 
  col=list(
  Celltype  =unlist(col_cells),
  Voxel = c(
    "Cort" = "#4F5088", 
    "Hipp" = "#A5458D", 
    "Put" = "#5D9E4F", 
    "WM" = "#4E77BC")
  ),
  show_annotation_name = F,
  show_legend =  T,
  annotation_legend_param = list(nrow=12)) 
col_dist    = dist(t(as.matrix(exprs)), method="canberra" )
coldend     = hclust(col_dist, method="ward.D2")
row_dist    = dist(as.matrix(exprs), method="canberra" )
rowdend     = hclust(row_dist, method="ward.D2")

hm <- Heatmap(exprs, 
        name = "Rel. expression",
        show_column_dend=T,
        show_column_names =T,
        show_row_names =F,
        show_row_dend=F,
        row_title="",
        cluster_columns =coldend,
        cluster_rows =rowdend,
        column_names_rot = 45,
        row_names_gp = grid::gpar(fontsize = 8),
        row_title_gp = grid::gpar(fontsize = 10),
        column_names_gp = grid::gpar(fontsize = 8),
        column_title_gp = grid::gpar(fontsize = 10),
        top_annotation=col_anno,
        width = unit(12, "cm"))
pdf(here::here(Sys.Date(), "Plots", "Figure3H_Heatmap_pathways_zscore.pdf"), width = 25, height = 5)
draw(hm)
invisible(dev.off())
Heatmap(exprs,
        name = "Rel. expression",
        show_column_dend=T,
        show_column_names =T,
        show_row_names =T,
        show_row_dend=T,
        row_title="Mitopathways",
        cluster_columns =coldend,
        cluster_rows =rowdend,
        column_names_rot = 45,
        row_names_gp = grid::gpar(fontsize = 2),
        row_title_gp = grid::gpar(fontsize = 10),
        column_names_gp = grid::gpar(fontsize = 6),
        column_title_gp = grid::gpar(fontsize = 10),
        top_annotation=col_anno,
        width = unit(8, "cm")
           )

rm(pseudobulk_mito_vox_cell_zscore, exprs, anno_df, col_anno, hm)