# setwd( "PLOTTING" )

library( ggplot2 )
library( SuppDists )
library( MASS )
library( Haplin )
## Registered S3 methods overwritten by 'ffbase':
##   method   from
##   [.ff     ff  
##   [.ffdf   ff  
##   [<-.ff   ff  
##   [<-.ffdf ff
source( "../common_variables.R" )
load( "../GxM_PoOxM_new_betas_mult_CLorP.RData" )

gene <- "8q21.3"

# promoter x 8q21.3
cur.plot.obj <- haptable.prom.res.all[[ gene ]]
cur.plot.obj.list <- lapply( unique( cur.plot.obj$stratum ), function(x){
  return( dplyr::filter( cur.plot.obj, stratum == x ) )
} )
names( cur.plot.obj.list ) <- c( "all", "1", "2", "3" )

cur.signif.plot <- plot.haplinSlide( cur.plot.obj.list, title = "", y.lim = c( 1/2,30 ) )

intercept.hline <- ( as.numeric( cur.signif.plot$data$haplos[ 3 ] ) -
    as.numeric( cur.signif.plot$data$haplos[ 2 ] ) ) / 2

y.labels <- c( 1/2, 1, 2, 5, 10, 20, 30 )

cur.signif.plot +
  geom_vline( xintercept = 1.5, linetype = 2 ) +
    scale_y_continuous( breaks = y.labels, labels = as.character( y.labels ) ) +
  xlab( "stratum" ) +
    theme( axis.text.x = element_text( angle = 0, size = 12 ),
            axis.text.y = element_text( size = 12 ),
            axis.line.y = element_line( color = "black" ),
            axis.title = element_text( size = 14 ),
            panel.background = element_rect( fill = "white", colour = "grey50" ),
            panel.grid.major.y = element_line( colour = "grey70", linetype = 2 ) )

ggsave( "GxM_CLorP_mult_betas_promoter-8q213.png", width = 15, height = 12, units = "cm" )
load( "../GxM_PoOxM_new_betas_mult_cntrl.RData" )

gene <- "TPM1"

# promoter x TPM1
cur.plot.obj <- haptable.prom.res.all[[ gene ]]
cur.plot.obj.list <- lapply( unique( cur.plot.obj$stratum ), function(x){
  return( dplyr::filter( cur.plot.obj, stratum == x ) )
} )
names( cur.plot.obj.list ) <- c( "all", "1", "2", "3" )

cur.signif.plot <- plot.haplinSlide( cur.plot.obj.list, title = "", y.lim = c( 1/8,20 ) )

intercept.hline <- ( as.numeric( cur.signif.plot$data$haplos[ 3 ] ) -
    as.numeric( cur.signif.plot$data$haplos[ 2 ] ) ) / 2

y.labels <- c( 1/8, 1/4, 1/2, 1, 2, 5, 10, 15, 20 )

cur.signif.plot +
  geom_vline( xintercept = 1.5, linetype = 2 ) +
    scale_y_continuous( breaks = y.labels, labels = as.character( y.labels ) ) +
  xlab( "stratum" ) +
    theme( axis.text.x = element_text( angle = 0, size = 12 ),
            axis.text.y = element_text( size = 12 ),
            axis.line.y = element_line( color = "black" ),
            axis.title = element_text( size = 14 ),
            panel.background = element_rect( fill = "white", colour = "grey50" ),
            panel.grid.major.y = element_line( colour = "grey70", linetype = 2 ) )

ggsave( "GxM_control_mult_betas_promoter-TPM1.png", width = 15, height = 12, units = "cm" )


gene <- "FOXE1"

# gene x FOXE1
cur.plot.obj <- haptable.gene.res.all[[ gene ]]
cur.plot.obj.list <- lapply( unique( cur.plot.obj$stratum ), function(x){
  return( dplyr::filter( cur.plot.obj, stratum == x ) )
} )
names( cur.plot.obj.list ) <- c( "all", "1", "2", "3" )

cur.signif.plot <- plot.haplinSlide( cur.plot.obj.list, title = "", y.lim = c( 1/8,15 ) )

intercept.hline <- ( as.numeric( cur.signif.plot$data$haplos[ 3 ] ) -
    as.numeric( cur.signif.plot$data$haplos[ 2 ] ) ) / 2

y.labels <- c( 1/8, 1/4, 1/2, 1, 2.5, 5, 10, 15 )

cur.signif.plot +
  geom_vline( xintercept = 1.5, linetype = 2 ) +
    scale_y_continuous( breaks = y.labels, labels = as.character( y.labels ) ) +
  xlab( "stratum" ) +
    theme( axis.text.x = element_text( angle = 0, size = 12 ),
            axis.text.y = element_text( size = 12 ),
            axis.line.y = element_line( color = "black" ),
            axis.title = element_text( size = 14 ),
            panel.background = element_rect( fill = "white", colour = "grey50" ),
            panel.grid.major.y = element_line( colour = "grey70", linetype = 2 ) )

ggsave( "GxM_control_mult_betas_gene-FOXE1.png", width = 15, height = 12, units = "cm" )

# this one is just to show comparison with the significan plot for CL/P:
gene <- "8q21.3"

# promoter x 8q21.3
cur.plot.obj <- haptable.prom.res.all[[ gene ]]
cur.plot.obj.list <- lapply( unique( cur.plot.obj$stratum ), function(x){
  return( dplyr::filter( cur.plot.obj, stratum == x ) )
} )
names( cur.plot.obj.list ) <- c( "all", "1", "2", "3" )

cur.signif.plot <- plot.haplinSlide( cur.plot.obj.list, title = "", y.lim = c( 1/4,4 ) )

intercept.hline <- ( as.numeric( cur.signif.plot$data$haplos[ 3 ] ) -
    as.numeric( cur.signif.plot$data$haplos[ 2 ] ) ) / 2

y.labels <- c( 1/4, 1/2, 1, 2, 4 )

cur.signif.plot +
  geom_vline( xintercept = 1.5, linetype = 2 ) +
    scale_y_continuous( breaks = y.labels, labels = as.character( y.labels ) ) +
  xlab( "stratum" ) +
    theme( axis.text.x = element_text( angle = 0, size = 12 ),
            axis.text.y = element_text( size = 12 ),
            axis.line.y = element_line( color = "black" ),
            axis.title = element_text( size = 14 ),
            panel.background = element_rect( fill = "white", colour = "grey50" ),
            panel.grid.major.y = element_line( colour = "grey70", linetype = 2 ) )

ggsave( "GxM_control_mult_betas_promoter-8q213.png", width = 15, height = 12, units = "cm" )
load( "../GxM_PoOxM_new_betas_mult_CLP.RData" )

gene <- "FOXE1"

# gene x FOXE1
cur.plot.obj <- haptable.gene.res.all[[ gene ]]
cur.plot.obj.list <- lapply( unique( cur.plot.obj$stratum ), function(x){
  return( dplyr::filter( cur.plot.obj, stratum == x ) )
} )
names( cur.plot.obj.list ) <- c( "all", "1", "2", "3" )

cur.signif.plot <- plot.haplinSlide( cur.plot.obj.list, title = "", y.lim = c( 1/32,20 ) )

intercept.hline <- ( as.numeric( cur.signif.plot$data$haplos[ 3 ] ) -
    as.numeric( cur.signif.plot$data$haplos[ 2 ] ) ) / 2

y.labels <- c( 1/32, 1/16, 1/8, 1/4, 1/2, 1, 2.5, 5, 10, 20 )

cur.signif.plot +
  geom_vline( xintercept = 1.5, linetype = 2 ) +
    scale_y_continuous( breaks = y.labels, labels = as.character( y.labels ) ) +
  xlab( "stratum" ) +
    theme( axis.text.x = element_text( angle = 0, size = 12 ),
            axis.text.y = element_text( size = 12 ),
            axis.line.y = element_line( color = "black" ),
            axis.title = element_text( size = 14 ),
            panel.background = element_rect( fill = "white", colour = "grey50" ),
            panel.grid.major.y = element_line( colour = "grey70", linetype = 2 ) )

ggsave( "GxM_CLP_mult_betas_gene-FOXE1.png", width = 15, height = 12, units = "cm" )