#CpGs from position 136669244 to 13666449 (the entire fragment of the intron 4 putative transcription starting site), 
#CpGs from position 136670325 to 136670386 (5’ of the intron 7) and 
#CpGs from position 136670827 to 136670945 (3’ of the intron 7).

# load data:

source("MethylAnalysis.R")

MethylAnalysis(methcall.folder = "input/", 
               outfolder = "Adults_vs_Prog", 
               samplesheet = "SampleSheets_BisulphiteExp_Sets.xlsx",
               sheetnum = 2,
               do.subset = T,
               chr.subset = "chr9",
               start.subset = 136668000,
               end.subset = 136671000,
               idcol = "SampleID", 
               design.var = "TestVar", 
               case.var = "Case", 
               control.var = "Controls",
               mincoverage = 10,
               idstoremove = c("S150310"),
               assem = "hg38", 
               pipeline.meth = "bismarkCoverage", 
               proj.id = "Adults_vs_Prog_Case_vs_Control", 
               plot.covariates = c("Condition","MRD","Cytogenetics","TestVar"), 
               lowperc = 5, 
               lowcount = 10, 
               delta.meth = 10, 
               plot.height = 12, plot.width = 18,
               plot.categorical.vars = c("Condition","Cytogenetics"), 
               plot.continuous.vars = c("miR-126"),cellw = 14, cellh = 8, fontnumsize = 7, fontsize = 7
)

library(openxlsx)
library(pheatmap)
library(RColorBrewer)

# load data pct methylation

pctmeth_matrix <- readRDS(file = "Adults_vs_Prog/Adults_vs_Prog_Case_vs_Control_CpG_percent_methylation_matrix_pheatmap.rds")

# fixing colors
annotations <- readRDS("Adults_vs_Prog/Adults_vs_Prog_Case_vs_Control_annotations_matrix_pheatmap.rds")
annotations2 <- annotations

annotations2$anncolors$Condition <- c("#006600","#929292","#c9c9c9","#b30101") # #006600: verdescuro, #b30101: rosso scuro, grigio-scuro: #929292, grigiochiaro: #c9c9c9
names(annotations2$anncolors$Condition) <- names(annotations$anncolors$Condition)

annotations2$anncolors$Cytogenetics <- c("#0080cc","#f6a145","white") # https://www.color-hex.com/color-palette/96440
names(annotations2$anncolors$Cytogenetics) <- names(annotations$anncolors$Cytogenetics)

png(filename = "Adults_vs_Prog/Adults_vs_Prog_Case_vs_Control_CpG_percent_methylation_matrix_pheatmap_CUSTOM.png", 
    width = 15, height = 7, units = "in", res = 300)
pheatmap(mat = pctmeth_matrix, 
         width = 15, height = 7,
         na_col = "black", 
         cluster_cols = FALSE, 
         cluster_rows = TRUE, 
         annotation_row = annotations2$annrow, 
         cellwidth = 14, 
         cellheight = 13, 
         display_numbers = T, 
         fontsize = 10,
         fontsize_number = 7, 
         number_format = "%.0f",
         border_color = FALSE, 
         annotation_col = annotations2$anncol, 
         legend = F,
         annotation_legend = TRUE,
         annotation_names_col = T,
         #labels_row = lrow,
         #labels_col = lcol,
         # gaps_row = c(23),
         gaps_col = c(9,43),
         annotation_colors = annotations2$anncolors,
         color = c("#F5F5F5","#EEEEEE","#CCCCCC","#999999", "#666666","#333333","#000000"), breaks = c(0,10,20,30,50,70,90,100)
)
dev.off()

#################
# Disease 12-27 #
#################

source("MethylAnalysis_bulk.R")

MethylAnalysis_bulk(methcall.folder = "input/", 
                    outfolder = "S120167_Case_vs_Controls", 
                    samplesheet = "SampleSheets_BisulphiteExp_Sets.xlsx",
                    sheetnum = 4,
                    do.subset = T,
                    chr.subset = "chr9",
                    start.subset = 136668000,
                    end.subset = 136671000,
                    idcol = "SampleID", 
                    design.var = "TestVar", 
                    case.var = "POP", 
                    control.var = "DGN",
                    mincoverage = 10,
                    idstoremove = NULL,
                    assem = "hg38", 
                    pipeline.meth = "bismarkCoverage", 
                    proj.id = "S120167_POP_vs_DGN", 
                    plot.covariates = c("Condition","TestVar"), 
                    lowperc = 5, 
                    lowcount = 10, 
                    delta.meth = 10, 
                    plot.height = 12, plot.width = 18,
                    plot.categorical.vars = c("Condition"), 
                    plot.continuous.vars = c("miR-126"),cellw = 14, cellh = 8, fontnumsize = 7, fontsize = 7
)

# RUNNING GFPLOW VS GFPHIGH to get diff meth results and stats

source("MethylAnalysis_bulk_v2.R") 

MethylAnalysis_bulk_v2(methcall.folder = "input/", 
                       outfolder = "S120167_GFP_L_vs_GFP_H", 
                       samplesheet = "SampleSheets_BisulphiteExp_Sets.xlsx",
                       sheetnum = 4,
                       do.subset = T,
                       chr.subset = "chr9",
                       start.subset = 136668000,
                       end.subset = 136671000,
                       idcol = "SampleID", 
                       design.var = "Condition", 
                       case.var = "GFP_L", 
                       control.var = "GFP_H",
                       mincoverage = 10,
                       idstoremove = NULL,
                       assem = "hg38", 
                       pipeline.meth = "bismarkCoverage", 
                       proj.id = "S120167_GFP_L_vs_GFP_H", 
                       plot.covariates = c("Condition"), 
                       lowperc = 5, 
                       lowcount = 10, 
                       delta.meth = 10, 
                       plot.height = 12, plot.width = 18,
                       plot.categorical.vars = c("Condition"), 
                       plot.continuous.vars = c("miR-126"),
                       cellw = 14, 
                       cellh = 8, 
                       fontnumsize = 7, 
                       fontsize = 7
)

# FIXING THE IMAGE

pctmeth_matrix <- readRDS(file = "S120167_Case_vs_Controls/S120167_POP_vs_DGN_CpG_percent_methylation_matrix_pheatmap.rds")

annotations <- readRDS("S120167_Case_vs_Controls/S120167_POP_vs_DGN_annotations_matrix_pheatmap.rds")
annotations_GFP_L_vs_GFP_H <- readRDS("S120167_GFP_L_vs_GFP_H/S120167_GFP_L_vs_GFP_H_annotations_matrix_pheatmap.rds")

# replacing diff meth and GFP_L vs GFP_H qvalue_r in POP vs DGN chart

annotations$anncol$qvalue_r <- annotations_GFP_L_vs_GFP_H$anncol$qvalue_r
annotations$anncol$meth.diff <- annotations_GFP_L_vs_GFP_H$anncol$meth.diff

annotations2 <- annotations

annotations2$anncolors$Condition <- c("#656161","#0b64a1","#df4041") # #006600: verdescuro, #b30101: rosso scuro, grigio-scuro: #929292, grigiochiaro: #c9c9c9
names(annotations2$anncolors$Condition) <- names(annotations$anncolors$Condition)

annotations2$anncolors$qvalue_r <- annotations_GFP_L_vs_GFP_H$anncolors$qvalue_r
annotations2$anncolors$meth.diff <- annotations_GFP_L_vs_GFP_H$anncolors$meth.diff

annotations2$annrow$`miR-126` <- NULL

# plot heatmap custom

png(filename = "S120167_Case_vs_Controls/S120167_GFP_Lvs_GFP_H_CpG_percent_methylation_matrix_pheatmap_CUSTOM_CpG_main.png", 
    width = 15, height = 5, units = "in", res = 300)
pheatmap(mat = pctmeth_matrix[,pctmeth_matrix_ADULTs_CpG],
         #filename = "S10272_Case_vs_Controls/S10272_GFP_Lvs_GFP_H_CpG_percent_methylation_matrix_pheatmap_CUSTOM.png", 
         width = 15, height = 7,
         na_col = "black", 
         cluster_cols = FALSE, 
         cluster_rows = TRUE, 
         annotation_row = annotations2$annrow, 
         cellwidth = 14, 
         cellheight = 13, 
         display_numbers = T, 
         fontsize = 10,
         fontsize_number = 7, 
         number_format = "%.0f",
         border_color = FALSE, 
         annotation_col = annotations2$anncol, 
         legend = F,
         annotation_legend = TRUE,
         annotation_names_col = T,
         #labels_row = lrow,
         #labels_col = lcol,
         # gaps_row = c(23),
         gaps_col = c(9,43),
         annotation_colors = annotations2$anncolors,
         color = c("#F5F5F5","#EEEEEE","#CCCCCC","#999999", "#666666","#333333","#000000"), breaks = c(0,10,20,30,50,70,90,100)
)
dev.off()