##################################################
##### integration of phQTL, mQTL, and eQTL #######
##################################################

# Set working directory and load libraries ####
  remove(list = ls())
  gc()
  
# load library
  library(dplyr)
  library(grid)
  library(scales)
  library(RColorBrewer)
  library(devtools)
  library(ggplot2)
  
# set directory
  work.dir <- "C:/Users/harta005/Projects/seed-germination-qtl"
  setwd(work.dir)
  
# gathering phQTL profile
  phqtl.peak <- readRDS('qtl-peaks/peak_phqtl_without-outliers-not-transformed.RDS') %>%
    mutate(qtl_level = 'phQTL') %>%
    mutate(stage = 'phQTL')

# gathering mQTL profile
  seed.stage <- c('pd', 'ar', 'im', 'rp')
  mqtl.peak <- data.frame(matrix(nrow = 0, ncol = 12))
  colnames(mqtl.peak) <- colnames(phqtl.peak)
  for (i in seed.stage) {
    table.tmp <- readRDS(paste0('qtl-peaks/peak_mqtl_', i, '-without-outliers-log-transformed.RDS')) %>%
      mutate(qtl_level = 'mQTL', stage = i)
      
    mqtl.peak <- rbind(mqtl.peak, table.tmp)
  }
  
# gathering eQTL profiles
  
  eqtl.peak <- data.frame(matrix(nrow = 0, ncol = 13))
  colnames(eqtl.peak) <- colnames(phqtl.peak)
  for (i in seed.stage) {
    table.tmp <- readRDS(paste0('qtl-tables/table_single-stage-eqtl_', i, '.rds')) %>%
      filter(qtl_type == 'trans') %>%
      mutate(qtl_level = 'eQTL', stage = i) %>%
      dplyr::select(colnames(phqtl.peak))
    eqtl.peak <- rbind(eqtl.peak, table.tmp)
  }

# creating the histogram plot ####
  
  # plot style ####
  presentation <- theme(axis.text.x = element_text(size=10, face="bold", color="black"),
                        axis.text.y = element_text(size=12, face="bold", color="black"),
                        axis.title.x = element_text(size=15, face="bold", color="black"),
                        axis.title.y = element_text(size=15, face="bold", color="black"),
                        strip.text.x = element_text(size=15, face="bold", color="black"),
                        strip.text.y = element_text(size=15, face="bold", color="black"),
                        plot.title = element_text(size=15, face="bold"),
                        panel.background = element_rect(fill = "white",color="black"),
                        panel.grid.major = element_line(colour = "grey80"),
                        panel.grid.minor = element_blank(),
                        legend.position = "right",
                        legend.title = element_text(size = 12, face = 'bold'),
                        legend.text = element_text(size = 12))
  
  myColors <- brewer.pal(9,"Set1")[c(3,5,9,6,7)] 
  myColors <- c("black",brewer.pal(9,"Set1")[c(2,9)])
  plot.colour <- c("#000000", "#E69F00", "#009E73", "#0072B2", "#D55E00")
  names(plot.colour) <- c("cis", "DryFresh", "DryAR",  "6H", "RP")
  hist.colour <-  c("#E69F00", "#009E73", "#0072B2", "#D55E00")
  names(hist.colour) <- c("DryFresh", "DryAR",  "6H", "RP")
  names(myColors) <- c("cis","trans","none")
  Snoekcols <- scale_colour_manual(name = "eQTL type",values = myColors)
  Snoekfill <- scale_fill_manual(name = "eQTL type",values = myColors)
  
# plotting ####
  
  overlap <- function(table, stage) {
    for (i in 1:length(stage)) {
      if (nrow(table) != 0) {
        table <- subset(table, table[stage[i]] == T)
      } else {
        break
      }
    }
    nrow(table)
  }
  
  all.peak <- rbind(phqtl.peak, mqtl.peak, eqtl.peak)
  all.peak$stage <- factor(all.peak$stage, levels = c('phQTL', 'pd', 'ar', 'im', 'rp'))
  qtl.hist <- ggplot(all.peak, aes(x=qtl_bp, fill = stage)) +
    geom_histogram(binwidth = 2000000, right = T, origin = 0, alpha = 1) +
    #geom_density(alpha = 0.5) +
    facet_grid(factor(qtl_level, level = c('phQTL', 'mQTL', 'eQTL')) ~ qtl_chromosome, scales="free") +
    
    presentation +
    scale_fill_manual(values = c('#000000', '#ccbb44', '#228833', '#4477aa', '#cc3311')) +
    theme(legend.position = "top") +
    labs(x="QTL peak position (Mb)",y="QTL counts") +
    scale_x_continuous(breaks=c(5, 10, 15, 20, 25, 30, 35, 40)*10^6,labels=c(5, 10, 15, 20, 25, 30, 35, 40))
  
  tiff(file = paste0("figures/all-qtl-hist.tiff"), 
       width = 2250, 
       height = 1600, 
       units = 'px',
       res = 300,
       compression = 'lzw')
  qtl.hist
  dev.off()
  
  # venn diagram ####
  
  qtl.table.summary <- data.frame(matrix(ncol = 18, nrow = 0))
  window.nu <- 2e6
  maxsize <- 100e6
  chr.num <- 5
  
  all.peak2 <- all.peak %>%
    mutate(interval = findInterval(qtl_bp, seq(1, maxsize, by = window.nu))) %>%
    select(trait, qtl_chromosome, interval, qtl_significance, stage, qtl_level) %>%
    #group_by(trait, qtl_chromosome, interval, qtl_significance, stage, qtl_level) %>%
    count(qtl_chromosome, interval, qtl_level) %>%
    #ungroup() %>%
    spread(qtl_level, n) %>%
    #select(phQTL, mQTL, eQTL) %>%
    replace(is.na(.), 0)
  
  overlap <- function(table, stage) {
    for (i in 1:length(stage)) {
      if (nrow(table) != 0) {
        table <- subset(table, table[stage[i]] > 0)
      } else {
        break
      }
    }
    nrow(table)
  }
  overlap(all.peak2, c('eQTL', 'phQTL', 'mQTL'))
  draw.triple.venn(area1 = overlap(all.peak2, 'phQTL'), 
                area3 = overlap(all.peak2, 'mQTL'), 
                area2 = overlap(all.peak2, 'eQTL'), 
                n13 = overlap(all.peak2, c('phQTL', 'mQTL')), 
                n12 = overlap(all.peak2, c('phQTL', 'eQTL')), 
                n23 = overlap(all.peak2, c('mQTL', 'eQTL')), 
                n123 = overlap(all.peak2, c('phQTL', 'mQTL', 'eQTL')), 
                category = c('phQTL', 'eQTL', 'mQTL'),
                lty = 'blank', 
                fill = c('#0d98ba', '#c71585', '#f8d568'))
                #alpha = 0.9)