---
title: "geing GSEA reactome plot divided by category"
author: ""
output:
  html_document:
  df_print: paged
word_document: default
pdf_document: default
---
  
  
```{r echo=FALSE,echo=FALSE, warning=FALSE,message=FALSE} 
library(ggplot2)
library(dplyr)
library(stringi)
library(tidyverse)
library(RColorBrewer)
library(circlize)
library(ComplexHeatmap)
library(reshape2)
library(scales)
```

```{r echo=FALSE,echo=FALSE, warning=FALSE,message=FALSE} 
folder = ('./data')

reactome <- read.csv(paste(folder,'ON_Old-ON_Young-gsea-Pathway.txt',sep = '/'),
                     sep='\t',   header = TRUE)


### plot the following path from reactome
dt <- reactome 
dt <- dt[dt$p.adjust <= 0.05,]
dt <- dt[order(dt$p.adjust, decreasing = FALSE),]

dt <- dt%>%  mutate(Color = ifelse(p.adjust >  mean(dt[,'p.adjust'] ), "red", "blue")) 
dt$logp.adjust <- -log10(dt$p.adjust)
dt$logp.adjust <- round(dt$logp.adjust,1)

dt <- dt[!dt$Description %in% c("Neuronal System","Neurotransmitter receptors and postsynaptic signal transmission"),]
dt <- dt[!dt$Description %in% c("Transmission across Chemical Synapses","Muscle contraction","Axon guidance","HIV Life Cycle" ,"HIV Infection","Late Phase of HIV Life Cycle","Host Interactions of HIV factors"),]
dt$name <- ifelse(dt$Description == "HDR through Homologous Recombination (HRR) or Single Strand Annealing (SSA)", 'HDR through HRR or SSA',dt$Description )
dt$Name <- dt$name #paste0(dt$ID, '-' ,dt$name)

l <- c('Hematopoiesis','Metabolic process','Gene expression regulation','DNA repair','Cell cycle',
       'Cell cycle','Chromatin conformation','DNA repair',
       'Cell cycle','Gene expression regulation','Chromatin conformation','Cellular organization and transport',
       'Protein Translation and localization', 'Cell cycle', 'Cell cycle','Metabolic process','Inflammation and Senescence','Gene expression regulation',
       'Cell cycle','Inflammation and Senescence','Gene expression regulation','Hematopoiesis',
       'Cell cycle','Cell cycle','Gene expression regulation','Inflammation and Senescence','Cell cycle','Inflammation and Senescence',
       'DNA repair','Metabolic process','Cellular organization and transport','Cell cycle',
       'Cell cycle','Protein Translation and localization',
       'Gene expression regulation')

dt$Category <- l

dt$Category <- factor(dt$Category, levels =c("Cell cycle",
                                             "DNA repair",
                                             "Inflammation and Senescence",
                                             "Gene expression regulation",
                                             "Chromatin conformation",
                                             "Protein Translation and localization",
                                             "Cellular organization and transport",
                                             "Metabolic process",
                                             "Hematopoiesis"))

m <- round(mean(seq(min(dt$logp.adjust),max(dt$logp.adjust),by = 0.1)),1)
p <- ggplot(dt, aes(x = NES, y = reorder(Name, NES), fill=logp.adjust)) + 
  geom_bar(stat = "identity",width=0.6)+ # aes(color= logp.adjust))+
  #geom_point(aes(size=Count, color= p.adjust))+
  #scale_colour_gradient(low="red", high="blue") +
  scale_fill_gradient(low="red",high="blue",
                      breaks = c(min(dt$logp.adjust),m,max(dt$logp.adjust)),
                      labels = c(min(dt$logp.adjust),m,max(dt$logp.adjust)),
                      limits = c(min(dt$logp.adjust),max(dt$logp.adjust)),
                      name = "-log10(FDR)") +
  cowplot::theme_cowplot() + labs(y= '',colour="-log10(FDR)")+
  #facet_grid(cols = vars(plot.0$Dataset), scales="free_y") +
  geom_vline(xintercept  = 0,linetype="dashed", color = "black",size=1)+
  theme(axis.title = element_text(size=6,color = 'black', face='bold'),
        axis.title.x  = element_text(size=12),
        axis.text.x = element_text(colour="black",size=12),#,face ='bold'),
        axis.text.y = element_text(colour="black",size=12),#,#face ='bold'),
        panel.grid.major.x = element_line(colour = "grey80", linetype = "dashed"),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_line(colour = "grey80"),
        panel.grid.minor.x = element_blank(), 
        legend.text=element_text(size=6, color ='black',face='bold'),
        legend.position="bottom",  legend.margin=margin(),#legend.box="vertical",
        strip.background = element_rect(colour="black",
                                        fill="white",size=0.5
        ),
        strip.text = element_text(color='black',face = 'bold', size = 8))+ ggtitle(paste0(''))

x <- p+ ggforce::facet_col(~Category,scales = "free_y",
                           space='free')
#strip.position = "right")

# svg(paste(folder, paste('ON_24',"gsea_Pathway_pos_HM_facet.svg", sep = "-"), sep = "/"), 
#     width=8,height=12,pointsize = 12)
# x #+ggtitle('REACTOME PATHWAY')
# dev.off()
```

```{r echo=FALSE,echo=FALSE, warning=FALSE,message=FALSE} 

###################################################################
#### check if genes in categories from core enrichment are signif
###################################################################
## signif Old vs Young ON (over night)
df_ss <- read.table('./data/DESeq-ON_Old-ON_Young-siggene.txt')
df_ss <- df_ss[df_ss$FDR < 0.05,]

real_sig <- list()
for (i in dt$Description){
  lg <- dt[dt$Description  == i,'core_enrichment']
  lg <- str_split(lg, "/")[[1]]
  real_sig[[i]] <- df_ss[rownames(df_ss) %in% lg,] 
}

## core enrichment not only signif genes
# df_ss <- read.table('/beegfs/scratch/ric.dimicco/ric.dimicco/DiMiccoR_1653_RNAseq/results/05-DGE-ON-no-Young2/tables/DESeq-ON_Old-ON_Young-allgene.txt')
# real <- list()
# for (i in hallmark.full_selected$ID){
#   lg <- hallmark.full_selected[hallmark.full_selected$ID == i,'core_enrichment']
#   lg <- str_split(lg, "/")[[1]]
#   real[[i]] <- df_ss[rownames(df_ss) %in% lg,] 
# }

## heatmap genes in core enrichment
doHeatmap2 <- function(df, sinfo, num_row_clu = 2, high_genes = NULL) {
  hm_matrix = as.matrix(df)
  cols = list()
  for (cc in colnames(sinfo)) {
    cols[[cc]] <- brewer.pal(9, "Set1")[1:length(unique(sinfo[[cc]]))]
    names(cols[[cc]]) <- rev(unique(sinfo[[cc]]))
  }
  ha_col <- columnAnnotation(df = sinfo,
                             col = cols,show_annotation_name = FALSE
                             #annotation_name_side = "right"
  )
  if (is.null(high_genes)) {
    row_names = TRUE
    ha_tre = NULL
  } else {
    row_names = FALSE
    gg <- intersect(rownames(hm_matrix), high_genes)
    gg_index <- which(rownames(hm_matrix) %in% gg)
    labels <- rownames(hm_matrix)[gg_index]
    ha_tre <- rowAnnotation(link = anno_mark(at = gg_index, labels = labels, labels_gp = gpar(fontsize = 7),
                                             link_width = unit(0, "mm"), 
                                             which = "row",side = "left")#,padding = unit(0.5, "mm"))
    )
  }
  hm = Heatmap(matrix = hm_matrix,
               #col = colorRampPalette(rev(brewer.pal(n = 9, name ="RdBu")))(100),
               #col = colorRamp2(c(-6, -4, -2, 0, 2, 4, 6), rev(brewer.pal(n = 7, name ="RdBu"))),
               #col = colorRamp2(c(-3, -2, -1, 0, 1, 2, 3), rev(brewer.pal(n = 7, name ="RdBu"))),
               col = colorRamp2(c(-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2), rev(brewer.pal(n = 9, name ="RdBu"))),
               cluster_columns = FALSE,#TRUE,
               cluster_rows = TRUE,
               rect_gp = gpar(col = "grey", lwd = 0),
               border = FALSE,
               show_row_names = TRUE,#FALSE,
               show_row_dend = FALSE,
               #show_row_names = FALSE, #row_names,
               row_names_gp = gpar(fontsize = 12),
               show_column_names = TRUE,
               column_title_gp = gpar(fontsize = 0),
               row_names_side = 'left',
               #row_split = as.numeric(levels(obj@meta.data[[clustering]])),
               #row_split = 
               #cluster_row_slices = TRUE,
               #column_split = sinfo$condition,
               #column_split = rep(0:(nrow(hm_matrix) - 1), each = ntop),
               #column_gap = unit(4, "mm"),
               #row_title_gp = gpar(fontsize = 12),
               #row_split = num_row_clu,
               #row_title_rot = 0,
               #left_annotation = ha_tre,
               top_annotation = ha_col,
               width = 10, height = 10, 
               heatmap_legend_param = list(
                 title = "Z-score"
               ))
  return(hm)
}

doHeatmap3 <- function(df, sinfo, num_row_clu = 2, high_genes = NULL, top_annot =NULL) {
  hm_matrix = as.matrix(df)
  rows = list()
  for (cc in colnames(sinfo)) {
    rows[[cc]] <- brewer.pal(9, "Set1")[1:length(unique(sinfo[[cc]]))]
    names(rows[[cc]]) <- rev(unique(sinfo[[cc]]))
  }
  ha_row <- rowAnnotation(df = sinfo,
                          col = rows,show_annotation_name = FALSE
                          #annotation_name_side = "right"
  )
  if (is.null(high_genes)) {
    col_names = TRUE
    ha_tre = NULL
  } else {
    row_names = FALSE
    gg <- intersect(colnames(hm_matrix), high_genes)
    gg_index <- which(colnames(hm_matrix) %in% gg)
    labels <- colnames(hm_matrix)[gg_index]
    ha_tre <- rowAnnotation(link = anno_mark(at = gg_index, labels = labels, labels_gp = gpar(fontsize = 7),
                                             link_width = unit(0, "mm"), 
                                             which = "row",side = "left")#,padding = unit(0.5, "mm"))
    )
  }
  hm = Heatmap(matrix = hm_matrix,
               #col = colorRampPalette(rev(brewer.pal(n = 9, name ="RdBu")))(100),
               #col = colorRamp2(c(-6, -4, -2, 0, 2, 4, 6), rev(brewer.pal(n = 7, name ="RdBu"))),
               #col = colorRamp2(c(-3, -2, -1, 0, 1, 2, 3), rev(brewer.pal(n = 7, name ="RdBu"))),
               col = colorRamp2(c(-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2), rev(brewer.pal(n = 9, name ="RdBu"))),
               cluster_columns = TRUE,#TRUE,
               cluster_rows =  FALSE,
               rect_gp = gpar(col = "black", lwd = 0),
               border = FALSE,
               show_row_names = FALSE,
               show_row_dend = FALSE,
               show_column_dend = FALSE,
               column_names_gp = gpar(fontsize = 12),
               show_column_names = TRUE,
               column_title_gp = gpar(fontsize = 0),
               #column_dend_side = 'top',
               column_names_rot = 75,
               #row_names_side = 'left',
               #row_split = as.numeric(levels(obj@meta.data[[clustering]])),
               #row_split = 
               #cluster_row_slices = TRUE,
               #column_split = sinfo$condition,
               #column_split = rep(0:(nrow(hm_matrix) - 1), each = ntop),
               #column_gap = unit(4, "mm"),
               #row_title_gp = gpar(fontsize = 12),
               #row_split = num_row_clu,
               #row_title_rot = 0,
               #left_annotation = ha_tre,
               top_annotation = top_annot,
               left_annotation = ha_row,
               width = 10, height = 10, 
               heatmap_legend_param = list(
                 title = "Z-score",
                 at = c(round(min(hm_matrix)),0,round(max(hm_matrix))),
                 legend_direction = "horizontal", 
                 legend_width = unit(4, "cm")),
  )
  return(hm)
}


folder <- "/beegfs/scratch/ric.dimicco/ric.dimicco/DiMiccoR_1653_RNAseq/bulk_rnaseq/data"
dimicco <- readRDS("data/ON_Old-ON_Young_20220331_175727_DGE.rds")

plot.dir <- folder
sinfo <-  dimicco$sample_info 
sinfo$condition <- factor(sinfo$condition, levels = c("ON_Young", "ON_Old" ))


## annotation list 
dt <- arrange(dt,Category)
cat_full <- data.frame(separate_rows(dt , core_enrichment,sep='/', convert = TRUE))
cat_full <- unique(cat_full[, c('Description','core_enrichment')])
table(unique(cat_full$core_enrichment))
length(unique(cat_full$core_enrichment)) 

annotcat <- dcast(cat_full, core_enrichment ~ Description)  ## condition
rownames(annotcat) <- annotcat$core_enrichment
annotcat <- annotcat[,-c(1)]
c <- colnames(annotcat)

for (val in c) {
  annotcat[[val]] <- ifelse(annotcat[[val]] != 'NA',val,NA)
  
}


##### all list together
gg <- c()
for (n in names(real_sig)) {
  gg <- append(gg,rownames(real_sig[[n]]))
}  
mat_gg <- as.data.frame(dimicco$deseq$rlog)[unique(gg), ]
annocat_signif <- annotcat[rownames(mat_gg),]

mat_gg_scaled <- t(scale(t(mat_gg)))
dimicco_hm <- doHeatmap2(df = t(mat_gg_scaled), sinfo = sinfo, high_genes =unique(gg))
rorder <- row_order(dimicco_hm)
mat_gg_scaled_names <- mat_gg_scaled[rorder,]

annocat_signif <- annocat_signif[,dt$Description]
annocat_signif  <- annocat_signif[rownames(mat_gg_scaled ),] # order automatically
table(rownames(annocat_signif ) == rownames(mat_gg_scaled))

n <- rep('black',dim(annocat_signif)[2])
names(n) <- colnames(annocat_signif)

my_cols <- hue_pal()(length(levels(dt$Category)))
names(my_cols) <- levels(dt$Category)
colors <- list()
for (i in dt$Description) {
  cn <- dt[dt$Description == i, 'Category']
  colors[i] <-  my_cols[[cn]]
}
ht3 = Heatmap(as.matrix(annocat_signif), name = " ",show_column_names = TRUE,
              show_row_names =FALSE,
              border = TRUE,row_names_gp = gpar(fontsize = 4),
              column_names_gp = gpar(fontsize = 6, rot = 90),
              rect_gp = gpar(col = "black", lwd = 0.5),
              col = colors ,na_col = "white",
              heatmap_legend_param = list(labels = dt$Category ))#c(rep('', dim(annocat_signif)[2]))) )

## hm of gene in core enrichment and also signif all categories
# svg(filename=paste(plot.dir, paste0("HM_full_categories", "_ON_Old_vs_Young.svg"), sep = "/"), 
#       width=12, 
#       height=25, 
#       pointsize=12)
#   draw(dimicco_hm+ht3, show_annotation_legend = FALSE, padding = unit(c(30, 2, 2, 60), "mm"))
#   dev.off()

####################################################### 
## hm of gene in core enrichment and also signif all categories reduced version
cat_full <- data.frame(separate_rows(dt , core_enrichment,sep='/', convert = TRUE))
cat_full <- unique(cat_full[, c('Category','core_enrichment')])
table(unique(cat_full$core_enrichment))
length(unique(cat_full$core_enrichment)) 
annotcat <- dcast(cat_full, core_enrichment ~ Category)  ## condition
rownames(annotcat) <- annotcat$core_enrichment
annotcat <- annotcat[,-c(1)]
c <- colnames(annotcat)

for (val in c) {
  annotcat[[val]] <- ifelse(annotcat[[val]] != 'NA',val,NA)
  
}

annocat_signif <- annotcat[rownames(mat_gg),]
#rorder <- row_order(dimicco_hm)
#mat_gg_scaled_names <- mat_gg_scaled[rorder,]

#annocat_signif <- annocat_signif[,dt$Category]
#annocat_signif  <- annocat_signif[rownames(mat_gg_scaled ),] # order automatically
table(rownames(annocat_signif ) == rownames(mat_gg_scaled))

my_cols <- c(hue_pal()(length(levels(dt$Category))))
names(my_cols) <- levels(dt$Category)
ht3 = Heatmap(as.matrix(annocat_signif), name = " ",show_column_names = TRUE,
              show_row_names =FALSE,
              border = TRUE,row_names_gp = gpar(fontsize = 12),
              column_names_gp = gpar(fontsize = 12, rot = 90),
              rect_gp = gpar(col = "black", lwd = 0.5),
              col = my_cols ,na_col = "white",
              heatmap_legend_param = list(labels = levels(dt$Category) ))#c(rep('', dim(annocat_signif)[2]))) )

# svg(filename=paste(plot.dir, paste0("HM_reduced_categories", "_ON_Old_vs_Young.svg"), sep = "/"), 
#     width=12, 
#     height=25, 
#     pointsize=12)
# draw(dimicco_hm+ht3, show_annotation_legend = FALSE, padding = unit(c(30, 2, 2, 60), "mm"))
# dev.off()

## rotate heatmap
for (i in dt$Category) {
  colors[i] <-  my_cols[[cn]]
}
ht3 = HeatmapAnnotation(df= annocat_signif,
                        gp = gpar(col = 'black',lwd=0.5),#which = 'col',
                        #row_names_gp = gpar(fontsize = 12, rot = 90),
                        #rect_gp = gpar(col = "black", lwd = 0.5),
                        border = TRUE,   na_col = "white",
                        col = list('colors' = c(colors)),#list('colors' = c(my_cols)), 
                        show_legend = FALSE)             

dimicco_hm <- doHeatmap3(df = t(mat_gg_scaled), sinfo = sinfo, high_genes =unique(gg), 
                         top_annot=ht3)
#draw(dimicco_hm)
# svg(filename=paste(plot.dir, paste0("HM_reduced_categories_reverse", "_ON_Old_vs_Young.svg"), sep = "/"),
#     width=26,#12,
#     height=6, #25,
#     pointsize=12)
# #bottom, left, top, right paddings.
# draw(dimicco_hm, show_annotation_legend = FALSE, padding = unit(c(30, 2, 2, 30), "mm"))
# dev.off()
svg(filename=paste(plot.dir, paste0("HM_reduced_categories_reverse", "_ON_Old_vs_Young_v1.svg"), sep = "/"),
    width=26,#12,
    height=6, #25,
    pointsize=12)
#bottom, left, top, right paddings.
draw(dimicco_hm, show_annotation_legend = FALSE, padding = unit(c(30, 2, 2, 30), "mm"),heatmap_legend_side = "top")
dev.off()
```
##### HM of genes in the core enrichment and also signif single HM
```{r echo=FALSE,echo=FALSE, warning=FALSE,message=FALSE} 
####################################################### 

for (y in unique(dt$Category)) { ## for each macro category
  gg <- c()
  for (n in names(real_sig)) {
    if (n %in% dt[dt$Category == y,'Description']){
      print(n)
      #gg <- rownames(real_sig[[n]])
      gg <- append(gg,rownames(real_sig[[n]]))
    }
  }
  gx <- unique(gg)
  mat_gg <- as.data.frame(dimicco$deseq$rlog)[gx, ]
  mat_gg_scaled <- t(scale(t(mat_gg)))
  dimicco_hm <- doHeatmap2(df = mat_gg_scaled, sinfo = sinfo)
  
  str(sinfo)
  n <- gsub(' ','_',gsub('/','_',n))
  # svg(filename=paste(plot.dir, paste0("HM_", y, "_ON_Old_vs_Young.svg"), sep = "/"), 
  #     width=6, 
  #     height=6, 
  #     pointsize=12)
  # draw(dimicco_hm)
  # dev.off()
  # 
  if (length(gx) > 39) {
    png(filename=paste(plot.dir, paste0("HM_", y, "_ON_Old_vs_Young.png"), sep = "/"), 
        units = "in", width = 8, height = 25, res = 300)
    draw(dimicco_hm, merge_legend = TRUE,heatmap_legend_side = 'top',column_title =y)
    dev.off()
  } else if (length(gx) > 22) { 
    png(filename=paste(plot.dir, paste0("HM_", y, "_ON_Old_vs_Young.png"), sep = "/"), 
        units = "in", width = 5, height = 8, res = 300)
    draw(dimicco_hm, merge_legend = TRUE,heatmap_legend_side = 'top',column_title =y)
    dev.off()
  } else  if (length(gx) > 10) {
    png(filename=paste(plot.dir, paste0("HM_", y, "_ON_Old_vs_Young.png"), sep = "/"), 
        units = "in", width = 5, height = 6, res = 300)
    draw(dimicco_hm, merge_legend = TRUE,heatmap_legend_side = 'top',column_title =y)
    dev.off()
  } else {
    png(filename=paste(plot.dir, paste0("HM_", y, "_ON_Old_vs_Young.png"), sep = "/"), 
        units = "in", width = 4, height = 4, res = 300)
    draw(dimicco_hm, merge_legend = TRUE,heatmap_legend_side = 'top',column_title =y)
    dev.off()
  }
}

```

