library(ggplot2)
library(tibble)
library(dplyr)
library(tidyr)
library(readr)
library(glue)
library(gprofiler2)

setwd(paste0(Sys.getenv("HOME"), "/code/malawi_transposon/script/R_SCRIPT/"))
source("modules/custom_gene_ontology.R")

################################################################################
# SV in body only
gene_list <- readLines("./local/gene_ontology/sv_bodyOnly.txt")

go_result <- run_GO_analysis(gene_list, "acalliptera", 0.01)
go_result$result$query <- ""
highlighted_terms <- go_result$result %>% filter(highlighted == TRUE) %>% .$term_id

p <- gostplot(go_result, capped = TRUE, interactive = FALSE)
pp <- publish_gostplot_custom(p, highlight_terms = highlighted_terms)
ggsave(filename = "./local/gene_ontology_recode/gProfiler_withSV.pdf", pp, width = 9, height = 2.8)
publish_gosttable(go_result, highlight_terms = highlighted_terms, use_colors = TRUE, 
  filename = "./local/gene_ontology_recode/gProfiler_withSV-table.pdf")

################################################################################
# no SV
gene_list <- readLines("./local/gene_ontology/sv_bodyOnly-None.txt")

go_result <- run_GO_analysis(gene_list, "acalliptera")
go_result$result$query <- ""
go_result$result$term_name[go_result$result$term_id == "GO:0006617"] <- "SRP-dependent cotranslational protein targeting to membrane"
highlighted_terms <- go_result$result %>% filter(highlighted == TRUE) %>% .$term_id

p <- gostplot(go_result, capped = TRUE, interactive = FALSE)
pp <- publish_gostplot_custom(p, highlight_terms = highlighted_terms)
ggsave(filename = "./local/gene_ontology_recode/gProfiler_noSV.pdf", pp, width = 9, height = 2.8)
publish_gosttable(go_result, highlight_terms = highlighted_terms, use_colors = TRUE, 
                  filename = "./local/gene_ontology_recode/gProfiler_noSV-table.pdf")

################################################################################
# no SV in body and flanking regions
gene_list <- readLines("./local/gene_ontology/sv_bodyAndFlank-None.txt")

go_result <- run_GO_analysis(gene_list, "acalliptera")
go_result$result$query <- ""
highlighted_terms <- go_result$result %>% filter(highlighted == TRUE) %>% .$term_id

p <- gostplot(go_result, capped = TRUE, interactive = FALSE)
pp <- publish_gostplot_custom(p, highlight_terms = highlighted_terms)
ggsave(filename = "./local/gene_ontology_recode/gProfiler_noSVbodyAndFlank.pdf", pp, width = 9, height = 2.8)
publish_gosttable(go_result, highlight_terms = highlighted_terms, use_colors = TRUE, 
                  filename = "./local/gene_ontology_recode/gProfiler_noSVbodyAndFlank-table.pdf")

################################################################################
# SVs in upstream region
gene_list <- readLines("./local/gene_ontology/sv_flankUpOnly.txt")

go_result <- run_GO_analysis(gene_list, "acalliptera")
go_result$result$query <- ""
highlighted_terms <- go_result$result %>% filter(highlighted == TRUE) %>% .$term_id

p <- gostplot(go_result, capped = TRUE, interactive = FALSE)
pp <- publish_gostplot_custom(p, highlight_terms = highlighted_terms)
ggsave(filename = "./local/gene_ontology_recode/gProfiler_flankUpOnly.pdf", pp, width = 9, height = 2.8)
publish_gosttable(go_result, highlight_terms = highlighted_terms, use_colors = TRUE, 
                  filename = "./local/gene_ontology_recode/gProfiler_flankUpOnly-table.pdf")

################################################################################
# SVs in flanking region
gene_list <- readLines("./local/gene_ontology/sv_flankDownOnly.txt")

go_result <- run_GO_analysis(gene_list, "acalliptera")
go_result$result$query <- ""
highlighted_terms <- go_result$result %>% filter(highlighted == TRUE) %>% .$term_id

p <- gostplot(go_result, capped = TRUE, interactive = FALSE)
pp <- publish_gostplot_custom(p, highlight_terms = highlighted_terms)
ggsave(filename = "./local/gene_ontology_recode/gProfiler_flankDownOnly.pdf", pp, width = 9, height = 2.8)
publish_gosttable(go_result, highlight_terms = highlighted_terms, use_colors = TRUE, 
                  filename = "./local/gene_ontology_recode/gProfiler_flankDownOnly-table.pdf")

################################################################################
# differentially expressed in testis, with SVs in body and flanking regions
gene_list <- readLines("./local/gene_ontology/DiffExpTestis-svSimple_bodyAndFlank.txt")
go_result <- run_GO_analysis(gene_list, "acalliptera")
go_result$result$query <- ""
highlighted_terms <- go_result$result %>% filter(highlighted == TRUE) %>% .$term_id
publish_gosttable(go_result, highlight_terms = highlighted_terms, use_colors = TRUE, 
                  filename = "./local/gene_ontology_recode/gProfiler_diffExpTestes-SVbodyAndFlank-table.pdf")

################################################################################
# differentially expressed in testis (whole set)
gene_list <- readLines("./local/gene_ontology/diffExp-testes.txt")
go_result <- run_GO_analysis(gene_list, "acalliptera")
go_result$result$query <- ""
highlighted_terms <- go_result$result %>% filter(highlighted == TRUE) %>% .$term_id
publish_gosttable(go_result, highlight_terms = highlighted_terms, use_colors = TRUE, 
                  filename = "./local/gene_ontology_recode/gProfiler_diffExpTestes-noSV-table.pdf")

################################################################################
# HISTONE MARK: H3K4me1
gene_list <- readLines("./local/enhancer/H3K4me1-genes_polymorphic.txt")
go_result <- run_GO_analysis(gene_list, "acalliptera", 0.1)
go_result$result$query <- ""
highlighted_terms <- go_result$result %>% filter(highlighted == TRUE) %>% .$term_id
publish_gosttable(go_result, highlight_terms = highlighted_terms, use_colors = TRUE, 
                  filename = "./local/gene_ontology_recode/gprofiler-H3K4me1-table.pdf")

################################################################################
# HISTONE MARK: H3K27ac
gene_list <- readLines("./local/enhancer/H3K27ac-genes_polymorphic.txt")
go_result <- run_GO_analysis(gene_list, "acalliptera", 0.1)
go_result$result$query <- ""
highlighted_terms <- go_result$result %>% filter(highlighted == TRUE) %>% .$term_id
publish_gosttable(go_result, highlight_terms = highlighted_terms, use_colors = TRUE, 
                  filename = "./local/gene_ontology_recode/gprofiler-H3K27ac-table.pdf")

################################################################################
# HISTONE MARK: H3K4me3
gene_list <- readLines("./local/enhancer/H3K4me3-genes_polymorphic.txt")
go_result <- run_GO_analysis(gene_list, "acalliptera", 0.1)
go_result$result$query <- ""
highlighted_terms <- go_result$result %>% filter(highlighted == TRUE) %>% .$term_id
publish_gosttable(go_result, highlight_terms = highlighted_terms, use_colors = TRUE, 
                  filename = "./local/gene_ontology_recode/gprofiler-H3K4me3-table.pdf")


