DEG-GO_analysis.Rmd

---
title: "GO analysis"
author: "ERM"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE)
```

# GO Analysis

I have created several files from the RNA analysis that contain the significant genes(determined by adj.P.val \< 0.05) from each Time and Condition. The names of the files are in the following format: 'sigV'+Drug(2 letters)+time.

example: 'sigVDA3.txt' means this file contains the significant DE genes from the Daunorubicin 3 hour compared to Vehicle Control 3 hour analysis

```{r Import libraries, message=FALSE, warning=FALSE}
library(tidyverse)
library(gprofiler2)
library(readr)
library(BiocGenerics)
library(gridExtra)
library(kableExtra)
library(scales)
library(ggVennDiagram)
library(cowplot)
library(ggpubr)
library(rstatix)
library(zoo)
library(ggsignif)
```

```{r Import data, echo=FALSE, message=FALSE, warning=FALSE}
toplistall <- read.csv("output/toplistall.csv")
  #get sig files made with 0.05 data this way(data created on run)
siglist <- readRDS("data/siglist.RDS")
list2env(siglist,envir=.GlobalEnv)
##readin the filenameonly.txt file from data to use for naming and filtering.
# 
# 
# filenameonly <- read.csv("data/filenameonly.txt")
# #loop through the list of files and make a separate dataframe for each file under the 'real'  name of the data set
# for (k in 1:length(file.names)){
#     assign(paste0(filenameonly$x[k]) , read.csv(file.names[k]))
# }
# 
# ##rename the columns to the previous names
# 
# colnames(sigVDA24)<- c("ENTREZID","SYMBOL","logFC","AveExpr","t","P.Value","adj.P.Val","B")
# colnames(sigVDX24)<- c("ENTREZID","SYMBOL","logFC","AveExpr","t","P.Value","adj.P.Val","B")
# colnames(sigVEP24)<- c("ENTREZID","SYMBOL","logFC","AveExpr","t","P.Value","adj.P.Val","B")
# colnames(sigVMT24)<- c("ENTREZID","SYMBOL","logFC","AveExpr","t","P.Value","adj.P.Val","B")
# colnames(sigVTR24)<- c("ENTREZID","SYMBOL","logFC","AveExpr","t","P.Value","adj.P.Val","B")
# colnames(sigVDA3)<- c("ENTREZID","SYMBOL","logFC","AveExpr","t","P.Value","adj.P.Val","B")
# colnames(sigVDX3)<- c("ENTREZID","SYMBOL","logFC","AveExpr","t","P.Value","adj.P.Val","B")
# colnames(sigVEP3)<- c("ENTREZID","SYMBOL","logFC","AveExpr","t","P.Value","adj.P.Val","B")
# colnames(sigVMT3)<- c("ENTREZID","SYMBOL","logFC","AveExpr","t","P.Value","adj.P.Val","B")
# colnames(sigVTR3)<- c("ENTREZID","SYMBOL","logFC","AveExpr","t","P.Value","adj.P.Val","B")
```

The analysis is based on all genes that passed the rowMeans\>0 from the previous page which are about 14084 genes expressed in my RNA-seq data [link](https://reneeisnowhere.github.io/Cardiotoxicity/run_all_analysis.html)

```{r uploading the background genes, echo=FALSE,  message=FALSE, warning=FALSE, paged.print=FALSE}
#backGL <- efit2$genes  # ***making the list

#write.csv(backGL, "data/backGL.txt", row.names = FALSE)
##read the  .csv then render into txt in the data file!
 backGL <- read_csv("data/backGL.txt", 
    col_types = cols(...1 = col_skip()))
```

All analysis is completed with an adjusted p value of 0.05
This is updated from the 0.1 value I had been using.

2023-04-19 RM

```{r tras had no DEGs}

drug_palNoVeh <- c("#8B006D" ,"#DF707E", "#F1B72B" ,"#3386DD", "#707031")


toplistall %>%
  mutate(id = as.factor(id)) %>%
  mutate(time=factor(time, levels=c("3_hours","24_hours"))) %>%
  group_by(time, id) %>% 
  mutate(sigcount = if_else(adj.P.Val < 0.05,'sig','notsig'))%>%
  count(sigcount) %>% 
  pivot_wider(id_cols = c(time,id), names_from=sigcount, values_from=n) %>% 
  mutate(prop = sig/(sig+notsig)*100) %>% 
  mutate(prop=if_else(is.na(prop),0,prop)) %>% 
  ggplot(., aes(x=id, y= prop))+
  geom_col(aes(fill=id))+
  geom_text(aes(label = sprintf("%.2f",prop)),
            position=position_dodge(0.9),vjust=-.2 )+
  scale_fill_manual(values =drug_palNoVeh)+
  guides(fill=guide_legend(title = "Treatment"))+
  facet_wrap(id~time)+#labeller = (time = facettimelabel) )+
  theme_bw()+
  xlab("")+
  ylab("Percentage of expressed genes")+
  theme_bw()+
  facet_wrap(~time)+
  ggtitle("Percent DEGs (adj. P value <0.05)")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        # axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        strip.background = element_rect(fill = "transparent"),
        axis.text.x = element_text(size = 8, color = "white", angle = 0),
        axis.text.y = element_text(size = 8, color = "black", angle = 0),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))


```
Knowels, et al. mentioned a relationship between transcriptional change level indicating cardiotoxicity.  This may be an observation that supports that claim as well.
## Daunorubicin

This is analysis on the significantly differentially expressed genes for Daunorubicin by 3 and 24 hours.
```{r intial analysis DA3 and DA24, echo=FALSE, message=FALSE, warning=FALSE}
# DAgostres3 <- gost(query = siglist$sigVDA3$ENTREZID,
#                     organism = "hsapiens",
#                     ordered_query = TRUE,
#                     domain_scope = "custom",
#                     measure_underrepresentation = FALSE,
#                     evcodes = FALSE,
#                     user_threshold = 0.2,
#                     correction_method = c("fdr"),
#                     custom_bg = backGL$ENTREZID,
#                     sources=c("GO:BP","KEGG"))
# saveRDS(DAgostres3,"data/DEG-GO/DAgostres3.RDS")
DAgostres3 <- readRDS("data/DEG-GO/DAgostres3.RDS")
DAp3 <- gostplot(DAgostres3, capped = FALSE, interactive = TRUE)
DAp3

DAtable3 <- DAgostres3$result %>%
  dplyr::select(c(source, term_id,term_name,intersection_size, 
                   term_size, p_value))
  
DAtable3 %>% 
    dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Daunorubicin 3 hour GO:BP terms') +
    xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 10, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))
DAtable3 %>% mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")
```

###  Daunorubicin specific genes at 3 hours 
    How they may relate to 24 hour genes using GO analysis
 3 hours specifically (see venn below too):

```{r daun 3 hour sig DEG,}

total24 <-list(sigVDA24$ENTREZID,sigVDX24$ENTREZID,sigVEP24$ENTREZID,sigVMT24$ENTREZID)
total3 <- list(sigVDA3$ENTREZID,sigVDX3$ENTREZID,sigVEP3$ENTREZID,sigVMT3$ENTREZID)

list3totvenn3 <- VennDiagram::get.venn.partitions(total3)
list24totvenn24 <- VennDiagram::get.venn.partitions(total24)

Daun3sp <- list3totvenn3$..values..[[15]]
Daun24sp <- list24totvenn24$..values..[[15]]


# DAspgostres3 <- gost(query = Daun3sp, 
#                     organism = "hsapiens",
#                     ordered_query = TRUE,
#                     domain_scope = "custom",
#                     measure_underrepresentation = FALSE,
#                     evcodes = FALSE,
#                     user_threshold = 0.05,
#                     correction_method = c("fdr"),
#                     custom_bg = backGL$ENTREZID,
#                     sources=c("GO:BP","KEGG"))
# saveRDS(DAspgostres3,"data/DEG-GO/DAspgostres3.RDS")
DAspgostres3 <- readRDS("data/DEG-GO/DAspgostres3.RDS")

DAspp3 <- gostplot(DAspgostres3, capped = FALSE, interactive = TRUE)
DAspp3

DAsptable3 <- DAspgostres3$result %>% 
  dplyr::select(c(source, term_id,
                  term_name,intersection_size, 
                   term_size, p_value))# %>% 
  
#---------------
```

###  Daunorubicin specific genes at 24 hours 

```{r DA24sp }
###no enrichment
# Daun24spgost <- gost(query = Daun24sp, 
#                    organism = "hsapiens",
#                     ordered_query = TRUE,
#                      domain_scope = "custom",
#                      measure_underrepresentation = FALSE,
#                      evcodes = FALSE,
#                      user_threshold = 0.05,
#                      correction_method = c("fdr"),
#                      custom_bg = backGL$ENTREZID,
#                      sources=c("GO:BP","KEGG"))
# saveRDS(DAspgostres3,"data/DEG-GO/DAspgostres3.RDS")
# DAspgostres3 <- readRDS("data/DEG-GO/DAspgostres3.RDS")
# 
# DAspp3 <- gostplot(DAspgostres3, capped = FALSE, interactive = TRUE)
# DAspp3
# 
# DAsptable3 <- DAspgostres3$result %>% 
#   dplyr::select(c(source, term_id,
#                   term_name,intersection_size, 
#                    term_size, p_value))# %>% 
#   
# 

```




### All Daun DEGs adj. p value < 0.05

```{r daun specific Venn genes}
# 3 hour daun specific genes (venn)
 


DAsptable3 <- DAspgostres3$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value))
  DAsptable3 %>% 
    dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Daunorubicin specific 3 hour GO:BP \nterms n = 335') +
    xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 10, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))
  
  
DAsptable3 %>% mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")



length(union(Daun3sp,sigVDX3$ENTREZID))
ggVennDiagram::ggVennDiagram(list(Daun3sp, sigVDA24$ENTREZID),
              category.names = c("3 hour\n Daun- specific","24 hours\n DEG p < 0.05"),
              show_intersect = FALSE,
              set_color = "black",
              label = "both",
              label_percent_digit = 1,
              label_size = 3.5,
              label_alpha = 0,
              label_color = "black",
              edge_lty = "solid") +
  scale_x_continuous(expand = expansion(mult = .2))+
  scale_fill_gradient(low = "red4", high = "tan1")+
  labs(title = "3 hour venn with 24 hour Daun", caption = "n = 351 genes")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5))


```
Interestingly, Over half of the daunorubicin specific genes are in the Daunorucin 24 hour group.

When I look at those 187 genes, I do not see a significant enrichment in GO:BP terms, only the KEGG pathway connection to Herpes simplex virus

```{r 3hour24hour daun overlap}
list324venn <- get.venn.partitions(list(Daun3sp, sigVDA24$ENTREZID))

over3_24daun <- list324venn$..values..[[1]]

# over3_24gost <- gost(query = over3_24daun,
#                   organism = "hsapiens",
#                   ordered_query = TRUE,
#                   domain_scope = "custom",
#                   measure_underrepresentation = FALSE,
#                   evcodes = FALSE,
#                   user_threshold = 0.05,
#                   correction_method = c("fdr"),
#                   custom_bg = backGL$ENTREZID,
#                   sources=c("GO:BP","KEGG"))
#saveRDS(over3_24gost, "data/DEG-GO/over3_24gost.RDS")
over3_24gost <- readRDS("data/DEG-GO/over3_24gost.RDS")
DA3_24table <- gostplot(over3_24gost, capped = FALSE, interactive = TRUE)
DA3_24table
plDA3_24table <- over3_24gost$result%>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value)) %>% 
  mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")
```
###Daunorubicin 24 sigDEG

Looking at the Daunorubicin 24 hour geneset using adj. P value of 0.05, I have enrichment  in the terms cell cycle and chromosome segregation. 
```{r DA sig deg 24}
# DAgostres24 <- gost(query = sigVDA24$ENTREZID,
#                   organism = "hsapiens",
#                   ordered_query = TRUE,
#                   domain_scope = "custom",
#                   measure_underrepresentation = FALSE,
#                   significant = FALSE,
#                   evcodes = FALSE,
#                   user_threshold = 0.05,
#                   correction_method = c("fdr"),
#                   custom_bg = backGL$ENTREZID,
#                   sources=c("GO:BP","KEGG"))
# saveRDS(DAgostres24, "data/DAgostres24.RDS")
DAgostres24 <- readRDS("data/DAgostres24.RDS")
DAp24 <- gostplot(DAgostres24, capped = FALSE, interactive = TRUE)
DAp24

DAtable124 <- DAgostres24$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value))  
  

  DAtable124 %>% 
  dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Daunorubicin 24 hour GO:BP terms') +
   xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+ 
    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 10, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))
  
  DAtable124 %>% mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")
  
#------------------------------------------------------
```
  
###Daunorubicin 3 v 24 intersection

All 187 3 hour Daun-specific genes are in the 304 genes between the two full lists of daunorubicin sigDEGs for 3 hours and 24 hours. This could possibly support the idea that Daunorubicin's more toxic effect is related to earlier impacts on cell function.  Note the enriched go terms for the 304 genes that they share.  I need to see if there is a trend in EP and DX.  I do see in Mito 3 hour (looked a little later).

```{r DA 3v24 int, echo=FALSE, message=FALSE, warning=FALSE}
daun3v24int <- intersect(sigVDA3$ENTREZID,sigVDA24$ENTREZID)
length(union(sigVDA3$ENTREZID,sigVDA24$ENTREZID))
ggVennDiagram::ggVennDiagram(list(sigVDA3$ENTREZID, sigVDA24$ENTREZID),
              category.names = c("3 hour Daunrubicin\nn=555","24 hours Daunorubicin\n   n=6864"),
              show_intersect = FALSE,
              set_color = "black",
              label = "both",
              label_percent_digit = 1,
              label_size = 3.5,
              label_alpha = 0,
              label_color = "black",
              edge_lty = "solid") +
  scale_x_continuous(expand = expansion(mult = .2))+
  scale_fill_gradient(low = "red4", high = "tan1")+
  labs(title = "DEG p < 0.05, All daunorubicin 3hour and 24 hours sets", caption = "n = 7115 genes")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5))


length(intersect(over3_24daun, intersect(sigVDA3$ENTREZID, sigVDA24$ENTREZID)))
# DAgostresint <- gost(query = daun3v24int,
#                     organism = "hsapiens",
#                     ordered_query = TRUE,
#                     domain_scope = "custom",
#                     measure_underrepresentation = FALSE,
#                     evcodes = FALSE,
#                     user_threshold = 0.05,
#                     correction_method = c("fdr"),
#                     custom_bg = backGL$ENTREZID,
#                     sources=c("GO:BP","KEGG"))

#saveRDS(DAgostresint,"data/DEG-GO/DAgostresint.RDS")
DAgostresint <- readRDS("data/DEG-GO/DAgostresint.RDS")
DApint <- gostplot(DAgostresint, capped = FALSE, interactive = TRUE)
DApint

DAtable_int <- DAgostresint$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value))  
 
  DAtable_int %>% 
    dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
    ggplot(., aes(x = log_val, y =reorder(term_name,p_value), 
                  col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", 
           size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Daunorubicin intersect 3 v 24 hour GO:BP terms') +
     xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
    theme_bw()+
    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
          axis.title = element_text(size = 15, color = "black"),
          axis.ticks = element_line(linewidth = 1.5),
          axis.line = element_line(linewidth = 1.5),
          axis.text = element_text(size = 10, color = "black", angle = 0),
          strip.text.x = element_text(size = 15, color = "black", face = "bold"))
  
  DAtable_int %>% 
    mutate_at(.vars = 6, .funs= scientific_format()) %>% 
    kable(.,) %>% 
    kable_paper("striped", full_width = FALSE) %>%
    kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>%
    scroll_box(width = "100%", height = "400px") 
  
  
  
  
```

## Doxorubicin 3 hour and 24 hour
```{r intial analysis DX3 and DX24, echo=FALSE, message=FALSE, warning=FALSE}
# DXgostres <- gost(query = sigVDX3$ENTREZID,  
#                   organism = "hsapiens",
#                   ordered_query = TRUE,
#                   domain_scope = "custom",
#                   measure_underrepresentation = FALSE,
#                   evcodes = FALSE,
#                   user_threshold = 0.05,
#                   correction_method = c("fdr"),
#                   custom_bg = backGL$ENTREZID,
#                   sources=c("GO:BP","KEGG"))
# saveRDS(DXgostres,"data/DEG-GO/DXgostres.RDS")

DXgostres <- readRDS("data/DEG-GO/DXgostres.RDS")
DXp <- gostplot(DXgostres, capped = FALSE, interactive = TRUE)
DXp

DXtable1 <- DXgostres$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value)) 


  DXtable1 %>% 
    dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Doxorubicin 3 hour GO:BP terms') +
   xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 9, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))

  DXtable1%>% mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")
  
 # ---------------------------------------------------------
  # 
  # DXgostres24 <- gost(query = sigVDX24$ENTREZID,  
  #                   organism = "hsapiens",
  #                   ordered_query = TRUE,
  #                   domain_scope = "custom",
  #                   measure_underrepresentation = FALSE,
  #                   evcodes = FALSE,
  #                   user_threshold = 0.05,
  #                   correction_method = c("fdr"),
  #                   custom_bg = backGL$ENTREZID,
  #                   sources=c("GO:BP","KEGG"))
#saveRDS(DXgostres24,"data/DEG-GO/gostres24.RDS") 
  DXgostres24 <- readRDS("data/DEG-GO/gostres24.RDS")
gostres24 <- readRDS("data/DEG-GO/gostres24.RDS")  
  p24 <- gostplot(gostres24, capped = FALSE, interactive = TRUE)
p24

DXtable124 <- gostres24$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value)) 


  DXtable124 %>% 
    dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected\n terms"))+
    ggtitle('Doxorubicin 24 hour GO:BP terms') +
   xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 10, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))
DXtable124 %>% 
  mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")
```
 
## Epirubicin  
All Epi DEGs adj. p value < 0.05


```{r intial analysis EP3 and EP24, echo=FALSE, message=FALSE, warning=FALSE}



# EPgostres <- gost(query = sigVEP3$ENTREZID,  
#                   organism = "hsapiens",
#                   ordered_query = TRUE,
#                   domain_scope = "custom",
#                   measure_underrepresentation = FALSE,
#                   evcodes = FALSE,
#                   user_threshold = 0.05,
#                   correction_method = c("fdr"),
#                   custom_bg = backGL$ENTREZID,
#                   sources=c("GO:BP","KEGG"))
# saveRDS(EPgostres,"data/DEG-GO/EPgostres.RDS")
EPgostres <- readRDS("data/DEG-GO/EPgostres.RDS")
EPp <- gostplot(EPgostres, capped = FALSE, interactive = TRUE)
EPp 

EPtable1 <- EPgostres$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value))  




  EPtable1 %>% 
    dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Epirubicin 3 hour top ten GO:BP terms') +
   xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 10, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))
  EPtable1%>%
    mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")
  
 # ---------------------------------------------------------
  
  # EPgostres24 <- gost(query = sigVEP24$ENTREZID,  
  #                   organism = "hsapiens",
  #                   ordered_query = TRUE,
  #                   domain_scope = "custom",
  #                   measure_underrepresentation = FALSE,
  #                   evcodes = FALSE,
  #                   user_threshold = 0.05,
  #                   correction_method = c("fdr"),
  #                   custom_bg = backGL$ENTREZID,
  #                   sources=c("GO:BP","KEGG"))
#saveRDS(EPgostres24,"data/DEG-GO/EPgostres24.RDS")
EPgostres24 <- readRDS("data/DEG-GO/EPgostres24.RDS")
  
  
  EPp24 <- gostplot(EPgostres24, capped = FALSE, interactive = TRUE)
EPp24

EPtable124 <- EPgostres24$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value)) 
EPtable124 %>% 
    dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Epirubicin 24 hour top 10 GO:BP terms') +
    xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 10, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))
EPtable124%>% 
  mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")



##epi 24 specifi n = 533
epsp24 <- list24totvenn24$..values..[[12]]





# EPgostres24sp <- gost(query = epsp24,
#                     organism = "hsapiens",
#                     ordered_query = TRUE,
#                     domain_scope = "custom",
#                     measure_underrepresentation = FALSE,
#                     evcodes = FALSE,
#                     user_threshold = 0.05,
#                     correction_method = c("fdr"),
#                     custom_bg = backGL$ENTREZID,
#                     sources=c("GO:BP","KEGG"))

#saveRDS(EPgostres24sp,"data/DEG-GO/EPgostres24sp.RDS")
EPgostres24sp <- readRDS("data/DEG-GO/EPgostres24sp.RDS")

EPsptable124 <- EPgostres24sp$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value)) 
EPsptable124%>% 
    dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle("Epirubicin specific 24 hour top GO:BP terms") +
     xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 10, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))
```
Epirubicin specific genes at 24 hours n= 489 genes are starting to show some apoptosis enrichment for myeloid cells. 


```{r mito specific at 24 hours}

mtsp24 <- list24totvenn24$..values..[[8]]
# MTgostres24sp <- gost(query = mtsp24,  
#                     organism = "hsapiens",
#                     ordered_query = TRUE,
#                     domain_scope = "custom",
#                     measure_underrepresentation = FALSE,
#                     evcodes = FALSE,
#                     user_threshold = 0.05,
#                     correction_method = c("fdr"),
#                     custom_bg = backGL$ENTREZID,
#                     sources=c("GO:BP","KEGG"))
# 
# saveRDS(MTgostres24sp,"data/DEG-GO/MTgostres24sp.RDS")
MTgostres24sp <- readRDS("data/DEG-GO/MTgostres24sp.RDS")


MTsptable124 <- MTgostres24sp$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value)) 


MTsptable124%>% 
    dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle("Mitoxantrone specific 24 hour top GO:BP terms") +
     xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 10, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))



```
Mitoxantrone-specific 24 hour genes show an enrichment for DNA DSB repair. (n=188)

## Mitoxantrone DEGs (adj. p value < 0.05)
  3 hour  = 58
  24 hour = 1327

```{r intial analysis MT3 and MT24, echo=FALSE, message=FALSE, warning=FALSE}
# MTgostres <- gost(query = sigVMT3$ENTREZID,
#                   organism = "hsapiens",
#                   ordered_query = TRUE,
#                   domain_scope = "custom",
#                   measure_underrepresentation = FALSE,
#                   evcodes = FALSE,
#                   user_threshold = 0.05,
#                   correction_method = c("fdr"),
#                   custom_bg = backGL$ENTREZID,
#                   sources=c("GO:BP","KEGG"))

# saveRDS(MTgostres,"data/DEG-GO/MTgostre.RDS")
MTgostres <- readRDS("data/DEG-GO/MTgostre.RDS")


MTp <- gostplot(MTgostres, capped = FALSE, interactive = TRUE)
MTp

MTtable1 <- MTgostres$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value))
MTtable1 %>% 
 dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Mitoxatrone 3 hour GO:BP terms') +
    xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 10, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold")) 
  MTtable1 %>% 
    mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")
 # ---------------------------------------------------------
#   
#   MTgostres24 <- gost(query = sigVMT24$ENTREZID,
#                     organism = "hsapiens",
#                     ordered_query = TRUE,
#                     domain_scope = "custom",
#                     measure_underrepresentation = FALSE,
#                     evcodes = FALSE,
#                     user_threshold = 0.05,
#                     correction_method = c("fdr"),
#                     custom_bg = backGL$ENTREZID,
#                     sources=c("GO:BP","KEGG"))
# saveRDS(MTgostres24,"data/DEG-GO/MTgostres24.RDS")
MTgostres24 <- readRDS("data/DEG-GO/MTgostres24.RDS")
  
  
  MTp24 <- gostplot(MTgostres24, capped = FALSE, interactive = TRUE)
MTp24

MTtable124 <- MTgostres24$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value)) 

MTtable124 %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Mitoxantrone 24 hour GO:BP terms') +
    xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 10, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))

  
  MTtable124 %>% 
    mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")
  

  
  
  
  

```
AT 3 hours, the Mitoxantrone DEGs, adj. p value < 0.01, look similar to daunorubicin gene sets a 3 hours.  hmmm



```{r up down trial data3, eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
#results_sig = subset(sigVDA3,adj.P.Val < 0.05)
# # get the significant up-regulated genes
#   up = subset(sigVDA3, logFC > 0)
# # get the significant down-regulated genes
# down = subset(sigVDA3, logFC < 0)
# gp_up = gost(query = up$ENTREZID,  organism = "hsapiens",
#   ordered_query = TRUE,
#   domain_scope = "custom",
#   measure_underrepresentation = FALSE,
#   evcodes = FALSE,
#   user_threshold = 0.001,
#   correction_method = c("fdr"),
#   custom_bg = backGL$ENTREZID,
#   sources=c("GO:BP","GO:MF", "GO:CC"))
#   
#   
# #gp_down = gost(query = down$ENTREZID,  organism = "hsapiens",
#   # ordered_query = TRUE,
#   # domain_scope = "custom",
#   # measure_underrepresentation = FALSE,
#   # evcodes = FALSE,
#   # user_threshold = 0.05,
#   # correction_method = c("bonferroni"),
#   # custom_bg = backGL$ENTREZID,
#   # sources=c("GO:BP","GO:MF", "GO:CC"))
# p2_up <- gostplot(gp_up, capped = FALSE, interactive = TRUE)
# 
# p2_up #+ ggtitle("Daunorubicin up regulated gene enrichment at 3 hours")

```

```{r showing it all, eval = FALSE, echo=FALSE, message=FALSE, warning=FALSE}

#p2_down <- gostplot(gp_down, capped = FALSE, interactive = TRUE)

#p2_down #+ ggtitle("Daunorubicin down regulated gene enrichment at 3 hours")


```

## Venn Diagrams 
### All DEGs adj. p value < 0.05, 24 hour
```{r venndiagram, echo=FALSE, message=FALSE, warning=FALSE}
library(paletteer)
total24 <-list(sigVDA24$ENTREZID,sigVDX24$ENTREZID,sigVEP24$ENTREZID,sigVMT24$ENTREZID)
in_common24 <-c(sigVDA24$ENTREZID,sigVDX24$ENTREZID,sigVEP24$ENTREZID,sigVMT24$ENTREZID)

length(unique(in_common24))

ggVennDiagram(total24,
              category.names = c("Daunorubicin\n n = 7017",
                                 "Doxorubicin\n n = 6645","Epirubicin\n n = 6328","Mitoxantrone\n n = 1115"),
              show_intersect = FALSE,
              set_color = "black",
              catagory_size = c(4,4,4,4),
              label = "count",
              label_percent_digit = 1,
              label_size = 3,
              label_alpha = 0,
              label_color = "black",
              edge_lty = "solid", set_size = 4.5)+
  scale_x_continuous(expand = expansion(mult = .3))+
   scale_y_continuous(expand = expansion(mult = .2))+
  
  scale_color_paletteer_d(palette = "fishualize::Bodianus_pulchellus")+
scale_fill_distiller(palette="Spectral", direction = -1)+
  #scale_hue_discrete(values=drug_palNoVeh)+
  labs(title = "24 hour comparison of DE genes p.adj. <0.05", caption = "n = 8188 genes")+
  theme(plot.title = element_text(size = rel(1.6), hjust = 0.5, vjust =1))



#paletteer_d("fishualize::Bodianus_pulchellus")
```

### All DEGs adj. p value < 0.05, 3 hour

```{r 3 hours, echo=FALSE, message=FALSE, warning=FALSE}
# total 3 ----------------------------------------------------------------
drug_palNoVeh <- c("#8B006D","#DF707E",
                   "#F1B72B","#3386DD","#707031","#41B333")
total3 <- list(sigVDA3$ENTREZID,sigVDX3$ENTREZID, sigVEP3$ENTREZID,sigVMT3$ENTREZID)
totalin_common3 <- c(sigVDA3$SYMBOL,sigVDX3$SYMBOL, sigVEP3$SYMBOL,sigVMT3$SYMBOL)
length(unique(totalin_common3))

ggVennDiagram(total3,
              category.names = c("Daunorubicin\nn = 555",
                                 "Doxorubicin\nn = 16",
                                 "Epirubicin\nn = 220",
                                 "Mitoxantrone\nn = 58"),
              show_intersect = FALSE,
              set_color = "black",
              catagory_size = c(6,6,6,6),
              label = "both",
              label_percent_digit = 1,
              label_size = 3,
              label_alpha = 0,
              label_color = "black",
              edge_lty = "solid", set_size = )+
  scale_x_continuous(expand = expansion(mult = .3))+
   scale_y_continuous(expand = expansion(mult = .2))+
scale_fill_distiller(palette="Spectral")+
  labs(title = "3 hour comparison of significant DE genes", caption = "n = 579 genes")+
  theme(plot.title = element_text(size = rel(1.6), hjust = 0.5, vjust =1))
```
Because the 3 hour DEGs are low (expecially in Doxorubicin),  we have chosen to focus on analyzing the 24 hour DEGs for GO: BP analysis
### Venn 3 and 24 hour sigDEGs by treatment (adj. p value <0.01) similarities

```{r 3V 24 BY DRUG, echo=FALSE}
# Dauno comp --------------------------------------------------------------

Dauncomp <- list(sigVDA24$ENTREZID,sigVDA3$ENTREZID)
in_commonDa <- c(sigVDA24$ENTREZID,sigVDA3$ENTREZID)
length(unique(in_commonDa))
ggVennDiagram(Dauncomp,
              category.names = c("Daunorubicin-24","Daunorubicin-3"),
              show_intersect = FALSE,
              set_color = "black",
              label = "both",
              label_percent_digit = 1,
              label_size = 3,
              label_alpha = 0,
              label_color = "black",
              edge_lty = "solid", set_size = )+
  scale_x_continuous(expand = expansion(mult = .2))+
  scale_fill_gradient(low = "light blue", high = "yellow")+
  labs(title = "Comparision of Dauno 3h v 24h", caption = "n = 8203 genes")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5))
#Davenlist <- intersect(Dauncomp[[1]],Dauncomp[[2]])


# Doxocomp ----------------------------------------------------------------
Doxcomp <- list(sigVDX24$ENTREZID,sigVDX3$ENTREZID)
in_commonDx <- c(sigVDX24$ENTREZID,sigVDX3$ENTREZID)
length(unique(in_commonDx))
ggVennDiagram(Doxcomp,
              category.names = c("Doxorubicin-24","Doxorubicin-3"),
              show_intersect = FALSE,
              set_color = "black",
              label = "both",
              label_percent_digit = 1,
              label_size = 3,
              label_alpha = 0,
              label_color = "black",
              edge_lty = "solid", set_size = )+
  scale_x_continuous(expand = expansion(mult = .2))+
  scale_fill_gradient(low = "light blue", high = "yellow")+
  labs(title = "Comparision of Doxo 3h v 24h", caption = "n = 7573 genes")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5))
#Dxvenlist <- intersect(Doxcomp[[1]],Doxcomp[[2]])
#length(intersect(Dxvenlist,Davenlist))## 7 of DX are in DA
# Epi Comp ----------------------------------------------------------------

Epicomp <- list(sigVEP24$ENTREZID,sigVEP3$ENTREZID)
in_commonEp <- c(sigVEP24$ENTREZID,sigVEP3$ENTREZID)
length(unique(in_commonEp))

ggVennDiagram(Epicomp,
              category.names = c("Epirubicin-24","Epirubicin-3"),
              show_intersect = FALSE,
              set_color = "black",
              label = "both",
              label_percent_digit = 1,
              label_size = 3,
              label_alpha = 0,
              label_color = "black",
              edge_lty = "solid", set_size = )+
  scale_x_continuous(expand = expansion(mult = .2))+
  scale_fill_gradient(low = "light blue", high = "yellow")+
  labs(title = "Comparision of Epi 3h v 24h", caption = "n = 7394 genes")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5))

#
# Mito comp ---------------------------------------------------------------
Mitocomp <- list(sigVMT24$ENTREZID,sigVMT3$ENTREZID)
in_commonMt <-c(sigVMT24$ENTREZID,sigVMT3$ENTREZID)
length(unique(in_commonMt))

ggVennDiagram(Mitocomp,
              category.names = c("Mitoxantrone-24","Mitoxantrone-3"),
              show_intersect = FALSE,
              set_color = "black",
              label = "both",
              label_percent_digit = 1,
              label_size = 3,
              label_alpha = 0,
              label_color = "black",
              edge_lty = "solid", set_size = )+
  scale_x_continuous(expand = expansion(mult = .21))+
  scale_fill_gradient(low = "light blue", high = "yellow")+
  labs(title = "Comparision of Mito 3h v 24h", caption = "n = 1642 genes")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5))



```

## GO analysis on 24 hour Venn, adj. p value < 0.05

```{r Making Go lists from SigDE genes, message=FALSE, warning=FALSE, paged.print=TRUE}

# > list24totvenn$..count..
#  [1]  977   15   17   14   58   24   34  188 4178  225  287  489  697  342  616
# > list24totvenn$..values..[?]

##pairwise
DDEresp<- list24totvenn24$..values..[[9]]
#DDEresp1<- read.csv("data/DDEresp_list.csv",row.names = 1)
# write.csv(DDEresp,"data/DDEresp_list.csv")
# DDEMresp<- list24totvenn24$..values..[[1]]
DDEMresp<-read.csv("data/DDEMresp_list.csv")
# write.csv(DDEMresp,"data/DDEMresp_list.csv")
# DDEall <-(intersect(sigVEP24$ENTREZID,intersect(sigVDA24$ENTREZID,sigVDX24$ENTREZID)))
# length(intersect(DDEresp,sigVMT24))
# 
# print("Number of AC drug DEG genes in common") 
# length(intersect(sigVEP24$ENTREZID,intersect(sigVDA24$ENTREZID,sigVDX24$ENTREZID)))
# #5317
DXsprespon<- list24totvenn24$..values..[[14]]
## (342)
print("Number of Doxorubicin specific DEGs, adj. P value < 0.05") 
length(DXsprespon)
complete <- c(sigVDA24$ENTREZID,sigVDX24$ENTREZID,sigVEP24$ENTREZID,sigVMT24$ENTREZID, sigVDA3$ENTREZID, sigVDX3$ENTREZID, sigVEP3$ENTREZID, sigVMT3$ENTREZID)
complete <- as.data.frame(unique(complete))
colnames(complete) <- "ENTREZID"  #9047 genes
print("Number of All genes that respond to some condition in DEG set, adj. P value  < 0.05")
length(complete$ENTREZID)
#8368
      
NoRespDEG <- setdiff(backGL$ENTREZID ,complete$ENTREZID) #
print("Number of No response genes in DEG set, adj. P value  < 0.05")
length(NoRespDEG) 

 #5716
```
## DDE response set

adj. P value  < 0.05, intersection of DNE-DOX-EPI, excluding those in common with MTX
n= 4435

no significant enrichment found!   I tried going for only the 
I did the intersection of all AC drugs and found a list that was 6129 genes long. 
Below is that analysis.

```{r GO DDEresp genelist of SDEGs,  message=FALSE, warning=FALSE, include=FALSE}

DDEresp <- list(DDEresp$x)#(ac +mtx)
DDEgostres24 <- gost(query = DDEresp,
                    organism = "hsapiens",
                    ordered_query = TRUE,
                    domain_scope = "custom",
                    measure_underrepresentation = FALSE,
                    evcodes = FALSE,
                    user_threshold = 0.1,
                    correction_method = c("fdr"),
                    custom_bg = backGL$ENTREZID,
                    sources=c("GO:BP","KEGG"))
# saveRDS(DDEgostres24,"data/DEG-GO/DDEgostres24-5155.RDS")
# DDEgostres24 <- readRDS("data/DEG-GO/DDEgostres24-5155.RDS")
# saveRDS(DDEgostres24only,"data/DEG-GO/DDEgostres24only.RDS")
DDEgostres24 <- readRDS("data/DEG-GO/DDEgostres24.RDS")
# saveRDS(DDEgostres24kegg,"data/DEG-GO/DDEgostres24kegg.RDS")
DDE24kegg <- readRDS("data/DEG-GO/DDEgostres24kegg.RDS")
DDEp24 <- gostplot(DDEgostres24only, capped = FALSE, interactive = TRUE)
DDEp24

kegg_term <- c("p53 signaling pathway","DNA replication","Cell cycle" ,"Base excision repair")         

DDEtable24kegg <- DDE24kegg$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value)) 
DDEtable24<- DDEgostres24only$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value)) 


DDEtable24 %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    dplyr::slice_min(., n= 15 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Anthracycline DEGs 24 hr GO:BP terms') +
    xlab(expression(" -"~log[10]~"(adj. p-value)"))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.2), hjust = 0.5),
        axis.title = element_text(size = 12, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 8, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))

DDEtable24kegg %>% dplyr::filter(source=="KEGG") %>%
    dplyr::select(p_value,term_name,intersection_size) %>%
  filter(term_name %in% kegg_term) %>% 
    slice_min(., n= 10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
     geom_point(aes(size = intersection_size, col = "red")) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Anthracycline DEGs 24 hr KEGG pathway terms') +
    xlab(expression(" -"~log[10]~"(adj. p-value)"))+
    ylab("KEGG term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.2), hjust = 0.5),
        axis.title = element_text(size = 12, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 8, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))


  
  DDEtable24 %>% 
    mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")

``` 

of the AC genes that overlapped without Mito (n=4178)
and all AC gens plus those which overlapped with Mito (n=5155), the significantly  enriched GO: BP terms are the same.  This may lend the idea that the specific AC set can tell us more about mechanism in AC induced cardiotoxicity.



## DDEM response of 977
```{r GO DDEM 977 genelist of SDEGs, echo=FALSE, message=FALSE, warning=FALSE}


# DDEMgostres24 <- gost(query = DDEMresp,  
#                     organism = "hsapiens",
#                     ordered_query = TRUE,
#                     domain_scope = "custom",
#                     measure_underrepresentation = FALSE,
#                     evcodes = FALSE,
#                     user_threshold = 0.05,
#                     correction_method = c("fdr"),
#                     custom_bg = backGL$ENTREZID,
#                     sources=c("GO:BP","KEGG"))

#saveRDS(DDEMgostres24,"data/DEG-GO/DDEMgostres24.RDS")
DDEMgostres24 <- readRDS("data/DEG-GO/DDEMgostres24.RDS")

DDEMp24 <- gostplot(DDEMgostres24, capped = FALSE, interactive = TRUE)
DDEMp24

DDEMtable24 <- DDEMgostres24$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value)) 

DDEMtable24 %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n= 10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Top2Bi DEGs (n=977) 24 hr GO:BP terms') +
    xlab(expression(" -"~log[10]~"(adj. p-value)"))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.2), hjust = 0.5),
        axis.title = element_text(size = 12, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 8, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))


DDEMtable24 %>% dplyr::filter(source=="KEGG") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n= 10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size, col = "red")) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Top2B inhibitor DEGs (n=977) 24 hr KEGG terms') +
    xlab(expression(" -"~log[10]~"(adj. p-value)"))+
    ylab("KEGG term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.2), hjust = 0.5),
        axis.title = element_text(size = 12, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 8, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))
  
  DDEMtable24 %>% 
    mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")
  
```

What I will do is take the terms from the tope three KEGG pathways and look at count changes or at FC values.  Updated later.
 

## DXspecific response (only the subset unique to Dox)

Whoa!  look at the calcium things
```{r GO DXsprespon$ENTREZID genelist of SDEGs, echo=FALSE, message=FALSE, warning=FALSE}

### no enrichment
# DXspgostres24 <- gost(query = DXsprespon, 
#                     organism = "hsapiens",
#                     ordered_query = TRUE,
#                     domain_scope = "custom",
#                     measure_underrepresentation = FALSE,
#                     evcodes = FALSE,
#                     user_threshold = 0.05,
#                     correction_method = c("fdr"),
#                     custom_bg = backGL$ENTREZID,
#                     sources=c("GO:BP","KEGG"))
# 
# saveRDS(DXspgostres24,"data/DEG-GO/DXspgostres24.RDS")
DXspgostres24 <- readRDS("data/DEG-GO/DXspgostres24.RDS")


DXspp24 <- gostplot(DXspgostres24, capped = FALSE, interactive = TRUE)
DXspp24

DXsptable24 <- DXspgostres24$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value)) 

DXsptable24 %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n= 10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Doxorubicin specific DEGs (n=342) 24 hr GO:BP terms') +
   xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.2), hjust = 0.5),
        axis.title = element_text(size = 12, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 8, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))

# 
# DXsptable24 %>% dplyr::filter(source=="KEGG") %>% 
#     dplyr::select(p_value,term_name,intersection_size) %>%
#     slice_min(., n= 10 ,order_by=p_value) %>%
#     mutate(log_val = -log10(p_value)) %>%
#    # slice_max(., n=10,order_by = p_value) %>%
#    ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
#     geom_point(aes(size = intersection_size, col = "red")) +
#     scale_y_discrete(labels = wrap_format(30))+
#     guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
#     ggtitle('Doxorubicin specific DEGs (n=400) 24 hr KEGG terms') +
#     xlab(expression(" -"~log[10]~("adj. p-value")))+
#     ylab("KEGG term")+
#       theme_bw()+
#     theme(plot.title = element_text(size = rel(1.2), hjust = 0.5),
#         axis.title = element_text(size = 12, color = "black"),
#         axis.ticks = element_line(linewidth = 1.5),
#         axis.line = element_line(linewidth = 1.5),
#         axis.text = element_text(size = 8, color = "black", angle = 0),
#         strip.text.x = element_text(size = 15, color = "black", face = "bold"))
  
  DXsptable24 %>% 
    mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")
  



``` 
I feel a little dumb for not remembering this results on Wednesday....   now that I ran this using the new adj. p value of < 0.05, I am more confident they are consistent.   as above, I will pull the specific genes and look at logCPM of these genes soon to examine if the change is seen in the other AC genes at 24 hours.


## No Response genes n=5727
 

```{r no response cluster, echo=FALSE, message=FALSE, warning=FALSE}
# 
# gostresNoResp <- gost(query = NoRespDEG,  
#                     organism = "hsapiens",
#                     ordered_query = FALSE,
#                     domain_scope = "custom",
#                     measure_underrepresentation = FALSE,
#                     evcodes = FALSE,
#                     user_threshold = 0.05,
#                     correction_method = c("fdr"),
#                     custom_bg = backGL$ENTREZID,
#                     sources=c("GO:BP","KEGG"))
# 
#saveRDS(gostresNoResp,"data/DEG-GO/gostresNoResp.RDS")
gostresNoResp <- readRDS("data/DEG-GO/gostresNoResp.RDS")

Noresp24 <- gostplot(gostresNoResp, capped = FALSE, interactive = TRUE)
Noresp24

Noresp24table <- gostresNoResp$result %>% 
  dplyr::select(c(source, term_id, term_name,intersection_size, term_size, p_value)) 

Noresp24table %>% dplyr::filter(source=="GO:BP") %>%
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n= 10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size)) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('No response genes(n=4593) 24 hr GO:BP terms') +
    xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("GO: BP term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.2), hjust = 0.5),
        axis.title = element_text(size = 12, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 8, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))


Noresp24table %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n= 10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
   # slice_max(., n=10,order_by = p_value) %>%
   ggplot(., aes(x = log_val, y =reorder(term_name,p_value), col= intersection_size)) +
    geom_point(aes(size = intersection_size, col = "red")) +
    scale_y_discrete(labels = wrap_format(30))+
    guides(col="none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('No response genes(n=4593) 24 hr KEGG terms') +
     xlab(expression(" -"~log[10]~("adj. p-value")))+
    ylab("KEGG term")+
      theme_bw()+
    theme(plot.title = element_text(size = rel(1.2), hjust = 0.5),
        axis.title = element_text(size = 12, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 8, color = "black", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))
  
  Noresp24table %>% 
    mutate_at(.vars = 6, .funs= scientific_format()) %>% 
  kable(.,) %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped","hover")) %>% 
  scroll_box(width = "100%", height = "400px")
  

  
```
##GO:BP full figure 24 DEG analysis


```{r GOBP dEGlong, echo=FALSE, fig.height=12, fig.width=8,message=FALSE, warning=FALSE}

DDEdeg <- DDEtable24 %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.)))) 

DDEMdeg <-  DDEMtable24 %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.))))

DAsp3deg <- DAsptable3%>% 
  dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.))))

DXsp24deg <- DXsptable24 %>% 
  dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=5 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.))))

EPsp24deg <- EPsptable124 %>% 
  dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.))))

MTsp24deg <-   MTsptable124 %>% 
    dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.))))
  

DAdeg <-   DAtable124 %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.))))

DXdeg <-  DXtable124 %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.))))

  EPdeg <-   EPtable124 %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
    mutate(order_val=rev(as.numeric(rownames(.))))
  
  MTdeg <-   MTtable124 %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
    mutate(order_val=rev(as.numeric(rownames(.))))

#filter<- dplyr::filter 
#select <- dplyr::select
  
GODEGlong <-list("DAdeg"=DAdeg,
                 "DXdeg"=DXdeg,
                 "DXsp24deg"=DXsp24deg,
                 "EPdeg"=EPdeg,
                 "EPsp24deg"=EPsp24deg,
                 "MTdeg"=MTdeg,
                 "MTsp24deg"= MTsp24deg,
                 "DDEdeg"=DDEdeg,
                 "DDEMdeg"=DDEMdeg)

GODEGlongtable <- data.table::rbindlist(GODEGlong, idcol = "deg")
 p <- GODEGlongtable%>% 
   mutate(deg = fct_inorder(deg)) %>% 
ggplot(., aes(x = log_val, y=reorder(term_name, order_val, desc=TRUE),col=deg)) +
   geom_point(aes(size = intersection_size, col=deg)) +
   scale_y_discrete(labels = wrap_format(40))+
   scale_color_manual(values=c("#8B006D","#DF707E","#e38893","#F1B72B","#f2c350","#3386DD","#33a1dd","red", "blue", alpha= 0.8))+
  facet_grid(deg ~.,scales="free_y" , labeller=label_bquote(.(as.expression(
               eval(parse(text = paste0('facet_deg', '$`', deg, '`')))
               )))) +
  
    guides(col=guide_legend(title="DEG Treatment group"), size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Top enriched GO:BP terms for each 24 hour DEG list') +
    xlab("-log 10 (adj. p-value)")+
    
    ylab("GO:BP term")+
      theme_bw()+
 theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 8, color = "black", angle = 0),
        #strip.text.x = element_text(size = 15, color = "black", face = "bold"),
        strip.text.y = element_text(color = "white"))

 
 
 facet_deg=(
   c(DAdeg = expression(atop("Daunorubicin 24 hr",~italic("n = 6864"))),
     DXdeg = expression(atop("Doxorubicin 24 hr",~italic("n = 6516"))),
     DXsp24deg =expression(atop("Dox-specific 24 hr",~italic("n = 342"))),
     EPdeg = expression(atop("Epirubicin 24 hr",~italic("n = 6202"))),
     EPsp24deg=expression(atop("Epi-specific 24 hr",~italic("n = 489"))),
     MTdeg = expression(atop( "Mitoxantrone 24 hr",~italic("n = 1327"))),
     MTsp24deg =expression(atop( "Mito-specific 24 hr",~italic("n = 188"))),
     DDEdeg = expression(atop("Anthracycline 24 hr",~italic("n = 4915,"))),
     DDEMdeg = expression(atop( "Top2Bi 24 hr",~italic("n = 977")))))

#  title_labeller_function <- function(nuisance_parameter) {
#   return(facet_degnames() )
# }
 
 
#https://stackoverflow.com/questions/41631806/change-facet-label-text-and-background-colour/60046113#60046113
 
# drug_palNoVeh <- c("#8B006D","#DF707E","#F1B72B", "#3386DD","#707031","#41B333")
 
 
g <- ggplot_gtable(ggplot_build(p))
stripr <- which(grepl('strip-r', g$layout$name))
fills <- c("#8B006D","#DF707E","#e38893","#F1B72B","#f2c350","#3386DD","#33a1dd","red", "blue", alpha= 0.8)
k <- 1
for (i in stripr) {
  j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
  g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
  k <- k+1
}
print(grid::grid.draw(g))
grid::grid.draw(g)
#hex <- hue_pal()(4)

#hex
#> hex (red/green/blue/purple default in ggplot)
#[1] "#F8766D" "#7CAE00" "#00BFC4" "#C77CFF"
#p



```



## CVRI poster figures not run

## Visualization of FC in cormotif sets
```{r Correlation of FC values in all values and in DEG}
library(limma)
library(corrplot)
library(pheatmap)

efit2<- readRDS(file ="data/efit2results.RDS")

# saveRDS(efit2, file ="data/efit2results.RDS")

V.DA.top= topTable(efit2, coef=1, adjust="BH", number=Inf, sort.by="p")
V.DX.top= topTable(efit2, coef=2, adjust="BH", number=Inf, sort.by="p")
V.EP.top= topTable(efit2, coef=3, adjust="BH", number=Inf, sort.by="p")
V.MT.top= topTable(efit2, coef=4, adjust="BH", number=Inf, sort.by="p")
V.TR.top= topTable(efit2, coef=5, adjust="BH", number=Inf, sort.by="p")
V.DA24.top= topTable(efit2, coef=6, adjust="BH", number=Inf, sort.by="p")
V.DX24.top= topTable(efit2, coef=7, adjust="BH", number=Inf, sort.by="p")
V.EP24.top= topTable(efit2, coef=8, adjust="BH", number=Inf, sort.by="p")
V.MT24.top= topTable(efit2, coef=9, adjust="BH", number=Inf, sort.by="p")
V.TR24.top= topTable(efit2, coef=10, adjust="BH", number=Inf, sort.by="p")

toplist3hours <-list(V.DA.top[,c(1,3)],
                              V.DX.top[,c(1,3)],
                              V.EP.top[,c(1,3)],
                              V.MT.top[,c(1,3)],
                              V.TR.top[,c(1,3)])
names(toplist3hours) <- c("Daunorubicin-3h","Doxorubicin-3h", "Epirubicin-3h","Mitoxantrone-3h", "Trastuzumab-3h")
library(RColorBrewer)
#toplist3hours <- map_df(toplist3hours, ~as.data.frame(.x), .id="id")

FCmatrix <- subset(efit2$coefficients)
is.matrix(subset(FCmatrix))
colnames(FCmatrix) <- c("Daunorubicin-3h","Doxorubicin-3h", "Epirubicin-3h","Mitoxantrone-3h", "Trastuzumab-3h","Daunorubicin-24h","Doxorubicin-24h", "Epirubicin-24h","Mitoxantrone-24h", "Trastuzumab-24h")


mat_col <- data.frame(time= c(rep("3 hours",5), rep("24 hours",5)))
rownames(mat_col) <- colnames(FCmatrix)
# list(group=c("#C25757", "#A77E69")
# names(mat_colors$group) <-c("3 hours", "24 hours")
mat_colors <- list(time=c("#C25757", "#A77E69"))
names(mat_colors$time) <- unique(mat_col$time)
corrFC <- cor(FCmatrix)
# corrplot(corrFC, method = 'square', order = 'AOE', addCoef.col = 'black')
# corrplot(corrFC,type = 'full', order = 'hclust', tl.col = 'black',
#          cl.ratio = 0.2, tl.srt = 45, col = COL2('PuOr', 20), col.lim=c(-0.1,1), addrect=3)


pheatmap::pheatmap(corrFC, display_numbers = TRUE, number_format = "%.2f",main=("Correlation of FC values across all expressed genes, n=14084"),annotation_col = mat_col, annotation_colors = mat_colors)



#making of the side columns for class
```


```{r Poster DEG GOBP, eval=FALSE, include=FALSE}
DDEdeg <- DDEtable24 %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=6 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%mutate(order_val=rev(as.numeric(rownames(.)))) %>% slice_max(., n=10, order_by=order_val)

DDEMdeg <-  DDEMtable24 %>% dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%mutate(order_val=rev(as.numeric(rownames(.))))
DAsp3deg <- DAsptable3%>% 
  dplyr::filter(source=="GO:BP") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    slice_min(., n=10 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.))))
MTsp24 <- MTsptable124 %>% 
  dplyr::filter(source=="GO:BP") %>% 
  dplyr::select(p_value,term_name,intersection_size) %>%
  slice_min(., n=10 ,order_by=p_value) %>%
  mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.))))





# GODEGposter <-list("DAsp3deg"=DAsp3deg,"DDEdeg"=DDEdeg,"DDEMdeg"=DDEMdeg)
GODEGposter <-list("MTsp24"=MTsp24,"DDEdeg"=DDEdeg,"DDEMdeg"=DDEMdeg)
GODEGpostertable <- data.table::rbindlist(GODEGposter, idcol = "deg")
facet_degp=(
   c(MTsp24 = expression(atop("Mitoxantrone-specific 24 hours",~italic("n = 126"))),
    DDEdeg=expression(atop("Shared Anthracycline genes 24 hours",~italic("n = 4400"))),
     DDEMdeg=expression(atop("Shared Top2"*beta*"i genes 24 hours",~italic("n = 952")))))

DEG_posterp <- GODEGpostertable%>% mutate(deg = fct_inorder(deg)) %>% 
  ggplot(., aes(x = log_val, y=reorder(term_name, order_val, desc=TRUE),col=deg)) +
    geom_point(aes(size = intersection_size, col=deg)) +
  scale_y_discrete(labels = wrap_format(40))+
   scale_color_manual(values=c("#3386DD","red3", "orange"))+
   facet_grid(deg ~.,scales="free_y" , labeller=label_bquote(.(as.expression(
               eval(parse(text = paste0('facet_degp', '$`', deg, '`')))
               )))) +
  guides(color= "none", size=guide_legend(title = "# of intersected \n terms"))+
    ggtitle('Top 10 enriched GO:BP terms\n in specific overlaping sets') +
     xlab(expression(" -"~log[10]~"(adj. p-value)"))+
      ylab("GO:BP term")+
      theme_bw()+
 theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text = element_text(size = 8, color = "black", angle = 0),
        #strip.text.x = element_text(size = 15, color = "black", face = "bold"),
        strip.text.y = element_text(color = "black"))

 
 g <- ggplot_gtable(ggplot_build(DEG_posterp))
stripr <- which(grepl('strip-r', g$layout$name))
fills <- c("#3386DD","red3", "orange", alpha= 0.8)
k <- 1
for (i in stripr) {
  j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
  g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
  k <- k+1
}
 

plot_grid(g)

```


# Genes of Interest log2 cpm
Genes of Interest:

Top2B  

Top2A  

CDKN1a  

ATM  

ATR  

mTOR  

Rictor  

Raptor  

RARG  



```{r GOI-top2b, echo=FALSE}

cpmcounts <- read.csv("data/filtered_cpm_counts.csv", row.names = 1)
drug_pal_vehend <- c("#8B006D","#DF707E","#F1B72B", "#3386DD","#707031","#41B333")
time <- rep((rep(c("3h", "24h"), c(6,6))), 6)
time <- ordered(time, levels =c("3h", "24h"))
drug <- rep(c("Daunorubicin","Doxorubicin","Epirubicin","Mitoxantrone","Trastuzumab", "Vehicle"),12)
indv <- as.factor(rep(c(1,2,3,4,5,6), c(12,12,12,12,12,12)))
labeltop <- (interaction(substring(drug, 0, 2), indv, time))
names(cpmcounts) <- labeltop

### TOP2B
cpmcounts %>% dplyr::filter(rownames(.)==7155) %>%
  pivot_longer(everything(), names_to = "treatment",values_to = "counts") %>%

  separate(treatment,c("drug","indv","time")) %>%

  mutate(time=factor(time, levels =c("3h", "24h"))) %>%
  mutate(indv=factor(indv, levels = c(1,2,3,4,5,6))) %>%
  mutate(drug =case_match(drug, "Da"~"Daunorubicin",
                          "Do"~"Doxorubicin",
                          "Ep"~"Epirubicin",
                          "Mi"~"Mitoxantrone",
                          "Tr"~"Trastuzumab",
                          "Ve"~"Vehicle", .default = drug))  %>%
  ggplot(., aes(x=drug, y=counts))+
  geom_boxplot(position="identity",aes(fill=drug))+
  geom_point(aes(col=indv, size=2, alpha=0.5))+
  guides(alpha= "none", size= "none")+
  scale_color_brewer(palette = "Dark2",guide = "none")+
  scale_fill_manual(values=drug_pal_vehend)+
  # geom_signif(comparisons =list(c("Daunorubicin","Vehicle"),
  #                               c("Doxorubicin","Vehicle"),
  #                               c("Epirubicin","Vehicle"),
  #                               c("Mitoxantrone","Vehicle"),
  #                               c("Trastuzumab","Vehicle")),
  #             test= "t.test",
  #             map_signif_level=FALSE,
  #             textsize =6,
  #             tip_length = .1,
  #             vjust = 0.2,
  #             step_increase = 0.1)+
  facet_wrap("time",  nrow=1, ncol=2)+
  theme_bw()+
  ylab(expression(atop(" ",italic("Top2"*beta)~log[2]~"cpm ")))+
  xlab("")+
  xlab("")+
  theme(strip.background = element_rect(fill = "#C77CFF"),
        plot.title = element_text(size=18,hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text.x = element_text(size = 12, color = "white", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))
```

```{r GOI CDKN1a, echo=FALSE}
### CDKN1a
cpmcounts %>% dplyr::filter(rownames(.)==1026) %>%
  pivot_longer(everything(), names_to = "treatment",values_to = "counts") %>%

  separate(treatment,c("drug","indv","time")) %>%

  mutate(time=factor(time, levels =c("3h", "24h"))) %>%
  mutate(indv=factor(indv, levels = c(1,2,3,4,5,6))) %>%
  mutate(drug =case_match(drug, "Da"~"Daunorubicin",
                          "Do"~"Doxorubicin",
                          "Ep"~"Epirubicin",
                          "Mi"~"Mitoxantrone",
                          "Tr"~"Trastuzumab",
                          "Ve"~"Vehicle", .default = drug))  %>%

   ggplot(., aes(x=drug, y=counts))+
  geom_boxplot(position="identity",aes(fill=drug))+
  geom_point(aes(col=indv, size=2, alpha=0.5))+
  guides(alpha= "none", size= "none")+
  scale_color_brewer(palette = "Dark2",guide = "none")+
  scale_fill_manual(values=drug_pal_vehend)+
  # geom_signif(comparisons =list(c("Daunorubicin","Vehicle"),
  #                               c("Doxorubicin","Vehicle"),
  #                               c("Epirubicin","Vehicle"),
  #                               c("Mitoxantrone","Vehicle"),
  #                               c("Trastuzumab","Vehicle")),
  #             test= "t.test",
  #             map_signif_level=FALSE,
  #             textsize =6,
  #             tip_length = .1,
  #             vjust = 0.2,
  #             step_increase = 0.1)+
  facet_wrap("time", nrow=1, ncol=2)+
  theme_bw()+
  ylab(expression(atop(" ",italic("CDKN1"*alpha)~log[2]~"cpm ")))+
  xlab("")+
  xlab("")+
  theme(strip.background = element_rect(fill = "#7CAE00"),
        plot.title = element_text(size=18,hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.line = element_line(linewidth = 1.5),
        axis.text.x = element_text(size = 12, color = "white", angle = 0),
        strip.text.x = element_text(size = 15, color = "black", face = "bold"))

```

```{r GOI-top2a, echo=FALSE}
### Top2a

cpmcounts %>% dplyr::filter(rownames(.)==7153) %>%
  pivot_longer(everything(), names_to = "treatment",values_to = "counts") %>%
  separate(treatment,c("drug","indv","time")) %>%
  mutate(time=factor(time, levels =c("3h", "24h"))) %>%
  mutate(indv=factor(indv, levels = c(1,2,3,4,5,6))) %>%
  mutate(drug =case_match(drug, "Da"~"Daunorubicin",
                          "Do"~"Doxorubicin",
                          "Ep"~"Epirubicin",
                          "Mi"~"Mitoxantrone",
                          "Tr"~"Trastuzumab",
                          "Ve"~"Vehicle", .default = drug))  %>%
 
  ggplot(., aes(x=drug, y=counts))+
  geom_boxplot(position="identity",aes(fill=drug))+
  geom_point(aes(col=indv, size=2, alpha=0.5))+
  guides(alpha= "none", size= "none")+
  scale_color_brewer(palette = "Dark2",guide = "none")+
  scale_fill_manual(values=drug_pal_vehend)+
  # geom_signif(comparisons =list(c("Daunorubicin","Vehicle"),
  #                               c("Doxorubicin","Vehicle"),
  #                               c("Epirubicin","Vehicle"),
  #                               c("Mitoxantrone","Vehicle"),
  #                               c("Trastuzumab","Vehicle")),
  #             test= "t.test",
  #             map_signif_level=FALSE,
  #             textsize =6,
  #             tip_length = .1,
  #             vjust = 0.2,
  #             step_increase = 0.1)+
  facet_wrap("time", nrow=1, ncol=2)+
  theme_bw()+
  ylab(expression(atop(" ",italic("Top2"*alpha)~log[2]~"cpm ")))+
  xlab("")+
  xlab("")+
  theme(strip.background = element_rect(fill = "white"),
    plot.title = element_text(size=18,hjust = 0.5),
    axis.title = element_text(size = 15, color = "black"),
    axis.ticks = element_line(linewidth = 1.5),
    axis.line = element_line(linewidth = 1.5),
    axis.text.x = element_text(size = 12, color = "white", angle = 0),
    strip.text.x = element_text(size = 15, color = "black", face = "bold"))
```

```{r GOI-atm, echo=FALSE}
### ATM

cpmcounts %>% dplyr::filter(rownames(.)==472) %>%
  pivot_longer(everything(), names_to = "treatment",values_to = "counts") %>%

  separate(treatment,c("drug","indv","time")) %>%

  mutate(time=factor(time, levels =c("3h", "24h"))) %>%
  mutate(indv=factor(indv, levels = c(1,2,3,4,5,6))) %>%
  mutate(drug =case_match(drug, "Da"~"Daunorubicin",
                          "Do"~"Doxorubicin",
                          "Ep"~"Epirubicin",
                          "Mi"~"Mitoxantrone",
                          "Tr"~"Trastuzumab",
                          "Ve"~"Vehicle", .default = drug))  %>%

  
  ggplot(., aes(x=drug, y=counts))+
  geom_boxplot(position="identity",aes(fill=drug))+
  geom_point(aes(col=indv, size=2, alpha=0.5))+
  guides(alpha= "none", size= "none")+
  scale_color_brewer(palette = "Dark2",guide = "none")+
  scale_fill_manual(values=drug_pal_vehend)+
  # geom_signif(comparisons =list(c("Daunorubicin","Vehicle"),
  #                               c("Doxorubicin","Vehicle"),
  #                               c("Epirubicin","Vehicle"),
  #                               c("Mitoxantrone","Vehicle"),
  #                               c("Trastuzumab","Vehicle")),
  #             test= "t.test",
  #             map_signif_level=FALSE,
  #             textsize =6,
  #             tip_length = .1,
  #             vjust = 0.2,
  #             step_increase = 0.1)+
  facet_wrap("time",  nrow=1, ncol=2)+
  theme_bw()+
  ylab(expression(atop(" ",italic("ATM")~log[2]~"cpm ")))+
  xlab("")+
  xlab("")+
  theme(#strip.background = element_rect(fill = "#C77CFF"),
    plot.title = element_text(size=18,hjust = 0.5),
    axis.title = element_text(size = 15, color = "black"),
    axis.ticks = element_line(linewidth = 1.5),
    axis.line = element_line(linewidth = 1.5),
    axis.text.x = element_text(size = 12, color = "white", angle = 0),
    strip.text.x = element_text(size = 15, color = "black", face = "bold"))
```

```{r GOI-atr, echo=FALSE}
###  ATR

cpmcounts %>% dplyr::filter(rownames(.)==545) %>%
  pivot_longer(everything(), names_to = "treatment",values_to = "counts") %>%

  separate(treatment,c("drug","indv","time")) %>%

  mutate(time=factor(time, levels =c("3h", "24h"))) %>%
  mutate(indv=factor(indv, levels = c(1,2,3,4,5,6))) %>%
  mutate(drug =case_match(drug, "Da"~"Daunorubicin",
                          "Do"~"Doxorubicin",
                          "Ep"~"Epirubicin",
                          "Mi"~"Mitoxantrone",
                          "Tr"~"Trastuzumab",
                          "Ve"~"Vehicle", .default = drug))  %>%

 
  ggplot(., aes(x=drug, y=counts))+
  geom_boxplot(position="identity",aes(fill=drug))+
  geom_point(aes(col=indv, size=2, alpha=0.5))+
  guides(alpha= "none", size= "none")+
  scale_color_brewer(palette = "Dark2",guide = "none")+
  scale_fill_manual(values=drug_pal_vehend)+
  # geom_signif(comparisons =list(c("Daunorubicin","Vehicle"),
  #                               c("Doxorubicin","Vehicle"),
  #                               c("Epirubicin","Vehicle"),
  #                               c("Mitoxantrone","Vehicle"),
  #                               c("Trastuzumab","Vehicle")),
  #             test= "t.test",
  #             map_signif_level=FALSE,
  #             textsize =6,
  #             tip_length = .1,
  #             vjust = 0.2,
  #             step_increase = 0.1)+
  facet_wrap("time", nrow=1, ncol=2)+
  theme_bw()+
  ylab(expression(atop(" ",italic("ATR")~log[2]~"cpm ")))+
  xlab("")+
  xlab("")+
  theme(#strip.background = element_rect(fill = "#C77CFF"),
    plot.title = element_text(size=18,hjust = 0.5),
    axis.title = element_text(size = 15, color = "black"),
    axis.ticks = element_line(linewidth = 1.5),
    axis.line = element_line(linewidth = 1.5),
    axis.text.x = element_text(size = 12, color = "white", angle = 0),
    strip.text.x = element_text(size = 15, color = "black", face = "bold"))
```

```{r GOI-rictor, echo=FALSE}
##Rictor


cpmcounts %>% dplyr::filter(rownames(.)==545) %>%
  pivot_longer(everything(), names_to = "treatment",values_to = "counts") %>%

  separate(treatment,c("drug","indv","time")) %>%

  mutate(time=factor(time, levels =c("3h", "24h"))) %>%
  mutate(indv=factor(indv, levels = c(1,2,3,4,5,6))) %>%
  mutate(drug =case_match(drug, "Da"~"Daunorubicin",
                          "Do"~"Doxorubicin",
                          "Ep"~"Epirubicin",
                          "Mi"~"Mitoxantrone",
                          "Tr"~"Trastuzumab",
                          "Ve"~"Vehicle", .default = drug))  %>%

 
  ggplot(., aes(x=drug, y=counts))+
  geom_boxplot(position="identity",aes(fill=drug))+
  geom_point(aes(col=indv, size=2, alpha=0.5))+
  guides(alpha= "none", size= "none")+
  scale_color_brewer(palette = "Dark2",guide = "none")+
  scale_fill_manual(values=drug_pal_vehend)+
  # geom_signif(comparisons =list(c("Daunorubicin","Vehicle"),
  #                               c("Doxorubicin","Vehicle"),
  #                               c("Epirubicin","Vehicle"),
  #                               c("Mitoxantrone","Vehicle"),
  #                               c("Trastuzumab","Vehicle")),
  #             test= "t.test",
  #             map_signif_level=FALSE,
  #             textsize =6,
  #             tip_length = .1,
  #             vjust = 0.2,
  #             step_increase = 0.1)+
  facet_wrap("time", nrow=1, ncol=2)+
  theme_bw()+
  ylab(expression(atop(" ",italic("RICTOR")~log[2]~"cpm ")))+
  xlab("")+
  xlab("")+
  theme(#strip.background = element_rect(fill = "#C77CFF"),
    plot.title = element_text(size=18,hjust = 0.5),
    axis.title = element_text(size = 15, color = "black"),
    axis.ticks = element_line(linewidth = 1.5),
    axis.line = element_line(linewidth = 1.5),
    axis.text.x = element_text(size = 12, color = "white", angle = 0),
    strip.text.x = element_text(size = 15, color = "black", face = "bold"))
```

```{r GOImtor, echo=FALSE}
### mTOR


cpmcounts %>% dplyr::filter(rownames(.)==2475) %>%
  pivot_longer(everything(), names_to = "treatment",values_to = "counts") %>%

  separate(treatment,c("drug","indv","time")) %>%

  mutate(time=factor(time, levels =c("3h", "24h"))) %>%
  mutate(indv=factor(indv, levels = c(1,2,3,4,5,6))) %>%
  mutate(drug =case_match(drug, "Da"~"Daunorubicin",
                          "Do"~"Doxorubicin",
                          "Ep"~"Epirubicin",
                          "Mi"~"Mitoxantrone",
                          "Tr"~"Trastuzumab",
                          "Ve"~"Vehicle", .default = drug))  %>%

 
  ggplot(., aes(x=drug, y=counts))+
  geom_boxplot(position="identity",aes(fill=drug))+
  geom_point(aes(col=indv, size=2, alpha=0.5))+
  guides(alpha= "none", size= "none")+
  scale_color_brewer(palette = "Dark2",guide = "none")+
  scale_fill_manual(values=drug_pal_vehend)+
  # geom_signif(comparisons =list(c("Daunorubicin","Vehicle"),
  #                               c("Doxorubicin","Vehicle"),
  #                               c("Epirubicin","Vehicle"),
  #                               c("Mitoxantrone","Vehicle"),
  #                               c("Trastuzumab","Vehicle")),
  #             test= "t.test",
  #             map_signif_level=FALSE,
  #             textsize =6,
  #             tip_length = .1,
  #             vjust = 0.2,
  #             step_increase = 0.1)+
  facet_wrap("time",  nrow=1, ncol=2)+
  theme_bw()+
  ylab(expression(atop(" ",italic("mTOR")~log[2]~"cpm ")))+
  xlab("")+
  xlab("")+
  theme(#strip.background = element_rect(fill = "#C77CFF"),
    plot.title = element_text(size=18,hjust = 0.5),
    axis.title = element_text(size = 15, color = "black"),
    axis.ticks = element_line(linewidth = 1.5),
    axis.line = element_line(linewidth = 1.5),
    axis.text.x = element_text(size = 12, color = "white", angle = 0),
    strip.text.x = element_text(size = 15, color = "black", face = "bold"))
```

```{r GOI RARG, echo=FALSE}
##RARG

cpmcounts %>% dplyr::filter(rownames(.)==5916) %>%
  pivot_longer(everything(), names_to = "treatment",values_to = "counts") %>%

  separate(treatment,c("drug","indv","time")) %>%

  mutate(time=factor(time, levels =c("3h", "24h"))) %>%
  mutate(indv=factor(indv, levels = c(1,2,3,4,5,6))) %>%
  mutate(drug =case_match(drug, "Da"~"Daunorubicin",
                          "Do"~"Doxorubicin",
                          "Ep"~"Epirubicin",
                          "Mi"~"Mitoxantrone",
                          "Tr"~"Trastuzumab",
                          "Ve"~"Vehicle", .default = drug))  %>%

 
  ggplot(., aes(x=drug, y=counts))+
  geom_boxplot(position="identity",aes(fill=drug))+
  geom_point(aes(col=indv, size=2, alpha=0.5))+
  guides(alpha= "none", size= "none")+
  scale_color_brewer(palette = "Dark2",guide = "none")+
  scale_fill_manual(values=drug_pal_vehend)+
  # geom_signif(comparisons =list(c("Daunorubicin","Vehicle"),
  #                               c("Doxorubicin","Vehicle"),
  #                               c("Epirubicin","Vehicle"),
  #                               c("Mitoxantrone","Vehicle"),
  #                               c("Trastuzumab","Vehicle")),
  #             test= "t.test",
  #             map_signif_level=FALSE,
  #             textsize =6,
  #             tip_length = .1,
  #             vjust = 0.2,
  #             step_increase = 0.1)+
  facet_wrap("time",  nrow=1, ncol=2)+
  theme_bw()+
  ylab(expression(atop(" ",italic("RARG")~log[2]~"cpm ")))+
  xlab("")+
  xlab("")+
  theme(#strip.background = element_rect(fill = "#C77CFF"),
    plot.title = element_text(size=18,hjust = 0.5),
    axis.title = element_text(size = 15, color = "black"),
    axis.ticks = element_line(linewidth = 1.5),
    axis.line = element_line(linewidth = 1.5),
    axis.text.x = element_text(size = 12, color = "white", angle = 0),
    strip.text.x = element_text(size = 15, color = "black", face = "bold"))


```


## expression correlation:


```{r expression correlation, echo=FALSE}
level_order2 <- c('75','87','77','79','78','71')
drug_pal_fact <- c("#8B006D","#DF707E","#F1B72B", "#3386DD","#707031","#41B333")

GOIcorr <-(genes = c("RARG", "ZNF740", "KAZN", "SLC28A3", "FRS2" ,"HDDC2","UTGT1A6","UGT1A10","KDM4B","KAT6B","GNT1"))
# ensembl <- useMart("ensembl", dataset="hsapiens_gene_ensembl")
# my_chr <- c(1:22, 'M', 'X', 'Y')  ## creates a filter for each
# 
# my_attributes <- c('entrezgene_id', 'ensembl_gene_id', 'hgnc_symbol')
# GOI_genelist<- getBM(attributes=my_attributes,filters ='hgnc_symbol',
#          values = GOIcorr, mart = ensembl)
# GOI_genelist <- GOI_genelist %>% distinct(entrezgene_id, .keep_all = TRUE)
# write.csv(GOI_genelist,"output/GOI_genelist.txt")
GOI_genelist <- read.csv("output/GOI_genelist.txt",row.names = 1)
RNAnormlist <- read.csv("output/TNNI_LDH_RNAnormlist.txt", row.names = 1)

GOI_join <- RNAnormlist %>% 
  mutate(indv=factor(indv,levels = level_order2)) %>% 
  mutate(indv=as.numeric(indv)) %>%
  mutate(indv=factor(indv)) %>% 
  mutate(Drug = factor(Drug, levels = c("Daunorubicin", 
                                       "Doxorubicin",
                                       "Epirubicin",
                                       "Mitoxantrone",
                                       "Trastuzumab",
                                       "Vehicle"))) %>% 
  dplyr::select(indv, Drug,rldh,rtnni)

#make a list of GOIs then individual info:

gene_corr_df <- cpmcounts %>% 
  rownames_to_column(var ="entrezgene_id") %>% 
  # filter(entrezgene_id %in% GOI_genelist$entrezgene_id) %>% ##this works to only fileter
  # inner_join(., GOI_genelist,by ="entrezgene_id") %>%  view()## this works to filter and append. just depends on the"when" I want to filter and join dataframes
  mutate(entrezgene_id=as.numeric(entrezgene_id)) %>% 
  pivot_longer(!entrezgene_id, names_to="treatment", values_to = "counts") %>%
  separate(treatment,c("Drug","indv","time")) %>%
  filter(time =="24h") %>%
  mutate(indv=factor(indv, levels = c(1,2,3,4,5,6))) %>%
  mutate(Drug =case_match(Drug, "Da"~"Daunorubicin",
                          "Do"~"Doxorubicin",
                          "Ep"~"Epirubicin",
                          "Mi"~"Mitoxantrone",
                          "Tr"~"Trastuzumab",
                          "Ve"~"Vehicle", .default = Drug)) 

  
# now I have the long dataframe, time to combine with other data snf limit to just the GOIs
  gene_corr_frame <- gene_corr_df %>% 
    filter(entrezgene_id %in% GOI_genelist$entrezgene_id) %>%  
    left_join(., GOI_join, by=c("Drug","indv")) %>% 
    dplyr::select(indv, Drug,entrezgene_id,counts, rldh,rtnni) %>% 
    mutate(Drug=factor(Drug)) %>% 
    full_join(.,GOI_genelist, by="entrezgene_id")
  
# saveRDS(gene_corr_frame,"data/gene_corr_frame.RDS")

   
   library(Hmisc)
  rcorr(as.matrix(gene_corr_frame[,c(1,4,5,6)]), type = "pearson")
   
for (gene in GOI_genelist$entrezgene_id){
    gene_plot <- gene_corr_frame %>% 
      dplyr::filter(entrezgene_id == gene) %>%
      ggplot(., aes(x=rtnni, y=counts))+
      geom_point(aes(col=indv))+
      geom_smooth(method="lm")+
      facet_wrap(hgnc_symbol~Drug, scales="free")+
      theme_classic()+
      xlab("troponin I expression") +
      ylab("Gene counts in log2 cpm") +
      ggtitle(expression(paste("Correlation between counts and troponin I by drug")))+
      scale_color_brewer(palette = "Dark2",name = "Individual", label = c("1","2","3","4","5","6"))+
      stat_cor(method="pearson",
               aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")),
               color = "red",
               label.x.npc = 0.01,
               label.y.npc=0.01, 
               size = 3)+
      theme(plot.title = element_text(size = rel(1.5), hjust = 0.5,face = "bold"),
            axis.title = element_text(size = 15, color = "black"),
            axis.ticks = element_line(size = 1.5),
            axis.text = element_text(size = 8, color = "black", angle = 20),
            strip.text.x = element_text(size = 12, color = "black", face = "italic"))
   print(gene_plot)
   
  }
  print(GOI_genelist)
## ggplot by indiv for tnni
  
  
  
 for (gene in GOI_genelist$entrezgene_id){
    gene_plot <- gene_corr_frame %>% 
      dplyr::filter(entrezgene_id == gene) %>%
      ggplot(., aes(x=rtnni, y=counts))+
      geom_point(aes(col=Drug))+
      geom_smooth(method="lm")+
      facet_wrap(hgnc_symbol~Drug, scales="free")+
      theme_classic()+
      xlab("troponin I expression") +
      ylab("Gene counts in log2 cpm") +
      ggtitle(expression(paste("Correlation between counts and troponin I by individual")))+
      scale_color_brewer(palette = "Dark2",name = "Individual", label = c("1","2","3","4","5","6"))+
      stat_cor(method="pearson",
               aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")),
               color = "red",
               label.x.npc = 0.01,
               label.y.npc=0.01, 
               size = 3)+
      theme(plot.title = element_text(size = rel(1.5), hjust = 0.5,face = "bold"),
            axis.title = element_text(size = 15, color = "black"),
            axis.ticks = element_line(size = 1.5),
            axis.text = element_text(size = 8, color = "black", angle = 20),
            strip.text.x = element_text(size = 12, color = "black", face = "italic"))
   print(gene_plot)
   
  }

print(GOI_genelist)

## ggplot by drug for ldh

  
  for (gene in GOI_genelist$entrezgene_id){
    gene_plot <- gene_corr_frame %>% 
      dplyr::filter(entrezgene_id == gene) %>%
      ggplot(., aes(x=rldh, y=counts))+
      geom_point(aes(col=indv))+
      geom_smooth(method="lm")+
      facet_wrap(hgnc_symbol~Drug, scales="free")+
      xlab("troponin I expression") +
      ylab("Gene counts in log2 cpm") +
      ggtitle(expression(paste("Correlation between counts and LDH by drug")))+
      scale_color_brewer(palette = "Dark2",name = "Individual", label = c("1","2","3","4","5","6"))+
      stat_cor(method="pearson",
               aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")),
               color = "red",
               label.x.npc = 0.01,
               label.y.npc=0.01, 
               size = 3)+
      theme(plot.title = element_text(size = rel(1.5), hjust = 0.5,face = "bold"),
            axis.title = element_text(size = 15, color = "black"),
            axis.ticks = element_line(size = 1.5),
            axis.text = element_text(size = 8, color = "black", angle = 20),
            strip.text.x = element_text(size = 12, color = "black", face = "italic"))
   print(gene_plot)
   
  }
  print(GOI_genelist)
## ggplot by indiv for ldh
  
  
  
 for (gene in GOI_genelist$entrezgene_id){
    gene_plot <- gene_corr_frame %>% 
      dplyr::filter(entrezgene_id == gene) %>%
      ggplot(., aes(x=rldh, y=counts))+
      geom_point(aes(col=Drug))+
      geom_smooth(method="lm")+
      facet_wrap(hgnc_symbol~Drug, scales="free")+
      theme_classic()+
      xlab("troponin I expression") +
      ylab("Gene counts in log2 cpm") +
      ggtitle(expression(paste("Correlation between counts and LDH by individual")))+
      scale_color_brewer(palette = "Dark2",name = "Individual", label = c("1","2","3","4","5","6"))+
      stat_cor(method="pearson",
               aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")),
               color = "red",
               label.x.npc = 0.01,
               label.y.npc=0.01, 
               size = 3)+
      theme(plot.title = element_text(size = rel(1.5), hjust = 0.5,face = "bold"),
            axis.title = element_text(size = 15, color = "black"),
            axis.ticks = element_line(size = 1.5),
            axis.text = element_text(size = 8, color = "black", angle = 20),
            strip.text.x = element_text(size = 12, color = "black", face = "italic"))
   print(gene_plot)
   
  }
GOIsig <- toplistall %>%
  filter(ENTREZID %in% GOI_genelist$entrezgene_id) %>%
  filter(adj.P.Val < 0.05) %>% 
  select(time,id,ENTREZID,SYMBOL,adj.P.Val) 
GOIsig%>% 
  kable(., caption= "Significant (adj. P value of <0.05) Genes of interest by treatment") %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped"),font_size = 18) %>% 
  scroll_box(width = "60%", height = "400px")
# write.csv(GOIsig,"data/GOIsig.csv")

print(GOI_genelist)
```
### Kegg 24 hour heatmap

```{r kegg 24 hour DEG heat}


DDEdegk <- DDEtable24 %>% dplyr::filter(source=="KEGG") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
  filter(term_name %in% kegg_term) %>% 
    slice_min(., n=4 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.)))) 

DDEMdegk <-  DDEMtable24 %>% dplyr::filter(source=="KEGG") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
  filter(term_name %in% kegg_term) %>% 
    slice_min(., n=3 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.))))

DAdegk <-   DAtable124 %>% dplyr::filter(source=="KEGG") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
  filter(term_name %in% kegg_term) %>% 
    slice_min(., n=3 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.))))

DXdegk <-  DXtable124 %>% dplyr::filter(source=="KEGG") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
  filter(term_name %in% kegg_term) %>% 
    slice_min(., n=3 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
  mutate(order_val=rev(as.numeric(rownames(.))))

  EPdegk <-   EPtable124 %>% dplyr::filter(source=="KEGG") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    filter(term_name %in% kegg_term) %>% 
    slice_min(., n=3 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
    mutate(order_val=rev(as.numeric(rownames(.))))
  
  MTdegk <-   MTtable124 %>% dplyr::filter(source=="KEGG") %>% 
    dplyr::select(p_value,term_name,intersection_size) %>%
    filter(term_name %in% kegg_term) %>% 
    slice_min(., n=3 ,order_by=p_value) %>%
    mutate(log_val = -log10(p_value)) %>%
    mutate(order_val=rev(as.numeric(rownames(.))))
library(ComplexHeatmap)
  
keggDEGlong <-list("DNR"=DAdegk,
                 "DOX"=DXdegk,
                 "EPI"=EPdegk,
                 "MTX"=MTdegk,
                 "Athracyclines"=DDEdegk,
                 "TOP2Bi"=DDEMdegk)

keggtable <- data.table::rbindlist(keggDEGlong, idcol = "deg")


kegg_mat<- keggtable%>% 
  dplyr::select(deg,log_val,term_name) %>% 
  pivot_wider(id_cols = everything(),names_from="term_name",values_from="log_val") %>% 
 column_to_rownames('deg') %>% 
  as.matrix()# 

Heatmap(kegg_mat, name = "KEGG Pathway -log p values", cluster_rows = FALSE, cluster_columns = FALSE)

  
```