Knowles2019.Rmd

---
title: "Knowles 2019 data"
author: "ERM"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
	echo = TRUE,
	message = FALSE,
	warning = FALSE
)
```
## Knowles data

```{r ensembl conversion, eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
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')
knowfig4 <- read.csv("data/knowfig4.csv",row.names = 1)
knowfig5 <- read.csv("data/knowfig5.csv",row.names = 1)

# covert ensg to entrezid -------------------------------------------------

knowles5 <- getBM(attributes=my_attributes,filters ='ensembl_gene_id',
                  values =knowfig5, mart = ensembl)
#377

knowles4 <- getBM(attributes=my_attributes,filters ='ensembl_gene_id',
                  values =knowfig4, mart = ensembl)
#524

knowles4 <-knowles4 %>% 
  distinct(entrezgene_id,.keep_all = TRUE)
  #521
knowles5 <-knowles5 %>% 
  distinct(entrezgene_id,.keep_all = TRUE) 
#377




saveRDS(knowles4,"output/knowles4.RDS")
saveRDS(knowles5,"output/knowles5.RDS")
  #Creating the Dox me/ re QTL data set:

DOXreQTLs <- toplistall %>% 
  dplyr::filter(time =="24_hours") %>% 
  mutate(sigcount = if_else(adj.P.Val < 0.05,'sig','notsig'))%>%
  mutate(meSNP=if_else(ENTREZID %in%knowles4$entrezgene_id,"y","no")) %>% 
  mutate(reQTL=if_else(ENTREZID %in%knowles5$entrezgene_id,"a","b")) %>% 
   dplyr::filter(id == "Doxorubicin") %>% 
   dplyr::filter(sigcount=='sig') %>% 
   dplyr::filter(reQTL=="a") %>% 
   dplyr::select(time, id, ENTREZID,SYMBOL)

saveRDS(DOXreQTLs,"output/DOXreQTLs.RDS")
DOXmeSNPs <- toplistall %>% 
  dplyr::filter(time =="24_hours") %>% 
  mutate(sigcount = if_else(adj.P.Val < 0.05,'sig','notsig'))%>%
  mutate(meSNP=if_else(ENTREZID %in%knowles4$entrezgene_id,"y","no")) %>% 
  mutate(reQTL=if_else(ENTREZID %in%knowles5$entrezgene_id,"a","b")) %>% 
   dplyr::filter(id == "Doxorubicin") %>% 
   dplyr::filter(sigcount=='sig') %>% 
   dplyr::filter(meSNP=="y") %>% 
   dplyr::select(time, id, ENTREZID,SYMBOL)
saveRDS(DOXmeSNPs,"output/DOXmeSNPs.RDS")
#Creating the DNR me/re QTL data sets:
DNRreQTLs <- toplistall %>% 
dplyr::filter(time =="24_hours") %>% 
  mutate(sigcount = if_else(adj.P.Val < 0.05,'sig','notsig'))%>%
  mutate(meSNP=if_else(ENTREZID %in%knowles4$entrezgene_id,"y","no")) %>% 
  mutate(reQTL=if_else(ENTREZID %in%knowles5$entrezgene_id,"a","b")) %>%
  dplyr::filter(id =="Daunorubicin") %>% 
  dplyr::filter(sigcount=='sig') %>% 
   dplyr::filter(reQTL=="a") %>% 
   dplyr::select(time, id, ENTREZID,SYMBOL)
saveRDS(DNRreQTLs,"output/DNRreQTLs.RDS")
DNRmeSNPs <- toplistall %>% 
dplyr::filter(time =="24_hours") %>% 
  mutate(sigcount = if_else(adj.P.Val < 0.05,'sig','notsig'))%>%
  mutate(meSNP=if_else(ENTREZID %in%knowles4$entrezgene_id,"y","no")) %>% 
  mutate(reQTL=if_else(ENTREZID %in%knowles5$entrezgene_id,"a","b")) %>%
  dplyr::filter(id =="Daunorubicin") %>% 
  dplyr::filter(sigcount=='sig') %>% 
   dplyr::filter(meSNP=="y") %>% 
   dplyr::select(time, id, ENTREZID,SYMBOL)
saveRDS(DNRmeSNPs,"output/DNRmeSNPs.RDS")

#Creating the Epi me/re QTL data sets:
EPIreQTLs <- toplistall %>% 
  dplyr::filter(time =="24_hours") %>% 
  mutate(sigcount = if_else(adj.P.Val < 0.05,'sig','notsig'))%>%
  mutate(meSNP=if_else(ENTREZID %in%knowles4$entrezgene_id,"y","no")) %>% 
  mutate(reQTL=if_else(ENTREZID %in%knowles5$entrezgene_id,"a","b")) %>%
   dplyr::filter(id =="Epirubicin") %>% 
   dplyr::filter(sigcount=='sig') %>% 
   dplyr::filter(reQTL=="a") %>% 
   dplyr::select(time, id, ENTREZID,SYMBOL)
saveRDS(EPIreQTLs,"output/EPIreQTLs.RDS")

EPImeSNPs <- toplistall %>% 
  dplyr::filter(time =="24_hours") %>% 
  mutate(sigcount = if_else(adj.P.Val < 0.05,'sig','notsig'))%>%
  mutate(meSNP=if_else(ENTREZID %in%knowles4$entrezgene_id,"y","no")) %>% 
  mutate(reQTL=if_else(ENTREZID %in%knowles5$entrezgene_id,"a","b")) %>%
   dplyr::filter(id =="Epirubicin") %>% 
   dplyr::filter(sigcount=='sig') %>% 
    dplyr::filter(meSNP=="y") %>% 
   dplyr::select(time, id, ENTREZID,SYMBOL)
saveRDS(EPImeSNPs,"output/EPImeSNPs.RDS")


#Creating the MTX me/re QTL data sets:

MTXreQTLs <- toplistall %>% 
  dplyr::filter(time =="24_hours") %>% 
  mutate(sigcount = if_else(adj.P.Val < 0.05,'sig','notsig'))%>%
  mutate(meSNP=if_else(ENTREZID %in%knowles4$entrezgene_id,"y","no")) %>% 
  mutate(reQTL=if_else(ENTREZID %in%knowles5$entrezgene_id,"a","b")) %>%
   dplyr::filter(id =="Mitoxantrone") %>% 
   dplyr::filter(sigcount=='sig') %>% 
   dplyr::filter(reQTL=="a") %>% 
   dplyr::select(time, id, ENTREZID,SYMBOL)
saveRDS(MTXreQTLs,"output/MTXreQTLs.RDS")

MTXmeSNPs <- toplistall %>% 
  dplyr::filter(time =="24_hours") %>% 
  mutate(sigcount = if_else(adj.P.Val < 0.05,'sig','notsig'))%>%
  mutate(meSNP=if_else(ENTREZID %in%knowles4$entrezgene_id,"y","no")) %>% 
  mutate(reQTL=if_else(ENTREZID %in%knowles5$entrezgene_id,"a","b")) %>%
   dplyr::filter(id =="Mitoxantrone") %>% 
   dplyr::filter(sigcount=='sig') %>% 
    dplyr::filter(meSNP=="y") %>% 
   dplyr::select(time, id, ENTREZID,SYMBOL)
saveRDS(MTXmeSNPs,"output/MTXmeSNPs.RDS")

##saved for later positioning

reQTLcombine=list(DNRreQTLs$ENTREZID,DOXreQTLs$ENTREZID,EPIreQTLs$ENTREZID,MTXreQTLs$ENTREZID)

length(unique(c(DNRreQTLs$ENTREZID,DOXreQTLs$ENTREZID,EPIreQTLs$ENTREZID,MTXreQTLs$ENTREZID)))

print("number of unique genes expressed in pairwise DE set")
QTLoverlaps <- VennDiagram::get.venn.partitions(reQTLcombine)
DoxonlyQTL <-  QTLoverlaps$..values..[[14]]
     ggVennDiagram::ggVennDiagram(reQTLcombine, category.names = c("DNR-187\nsigDEexpressed","DOX-180 \nsigDEexpressed","EPI-176\nsigDEexpressed","MTX-40\nsigDEexpressed"))  

```


``` {r loading in data sets}

backGL <- read.csv("data/backGL.txt",     row.names =1)
drug_palNoVeh <- c("#8B006D","#DF707E","#F1B72B", "#3386DD","#707031")
drug_palc <- 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)
mat_drug <- c("Daunorubicin","Doxorubicin","Epirubicin","Mitoxantrone")
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))
level_order2 <- c('75','87','77','79','78','71')


toplistall <- readRDS("data/toplistall.RDS")
knowles4 <-readRDS("output/knowles4.RDS")
knowles5 <-readRDS("output/knowles5.RDS")
DOXmeSNPs <- readRDS("output/DOXmeSNPs.RDS")
DOXreQTLs <- readRDS("output/DOXreQTLs.RDS")
DNRmeSNPs <- readRDS("output/DNRmeSNPs.RDS")
DNRreQTLs <- readRDS("output/DNRreQTLs.RDS")
EPImeSNPs <- readRDS("output/EPImeSNPs.RDS")
EPIreQTLs <- readRDS("output/EPIreQTLs.RDS")
MTXmeSNPs <- readRDS("output/MTXmeSNPs.RDS")
MTXreQTLs <- readRDS("output/MTXreQTLs.RDS")

toplist24hr <-   toplistall %>%
  mutate(id = as.factor(id)) %>%
  mutate(time=factor(time, levels=c("3_hours","24_hours"))) %>%
  filter(time=="24_hours")

toplist3hr <-   toplistall %>%
  mutate(id = as.factor(id)) %>%
  mutate(time=factor(time, levels=c("3_hours","24_hours"))) %>%
  filter(time=="3_hours")

```
Determining the genetic basis of anthracycline-cardiotoxicity by molecular response QTL mapping in induced cardiomyocytes  
David A Knowles, Courtney K Burrows†, John D Blischak,
Kristen M Patterson, Daniel J Serie, Nadine Norton, Carole Ober,
Jonathan K Pritchard, Yoav Gilad

Knowles $~~et~ al.~$ eLife 2018;7:e33480. DOI: https://doi.org/10.7554/eLife.33480  

My first question was about transcriptional response at the 24 hour mark with my treatments.   3 hour RNA-seq had low levels of DEGs,so my focus is at 24 hours.  This also happens to be when the Knowles paper collected their RNA-seq data


Supplementary 4 contains a list of 518 SNPs within 1 Mb of TSS, which had a detectable marginal effect on expression (5% FDR). When converted from ensembl gene id to entrez gene id,  my list of genes for this supplement = `r length(knowles4$entrezgene_id)`.  I will call these meSNPS for marginal effect QTLs.  

Supplementary 5 contains a list of 376 response eQTLs (reQTLs).  These are variants that were associated with modulation of transcriptomic response to doxorubicin treatment. After database name conversion, I have `r length(knowles5$entrezgene_id)` genes.
``` {r QTL enrichment numbers}

meSNPdf <- toplistall %>% 
  mutate(id = as.factor(id)) %>%
  mutate(time=factor(time, levels=c("3_hours","24_hours"))) %>%
  mutate(sigcount = if_else(adj.P.Val < 0.05,'sig','notsig'))%>%
  mutate(meSNP = if_else( ENTREZID %in% knowles4$entrezgene_id, "y" , "no")) %>% 
  dplyr::select(ENTREZID,id,time,meSNP,sigcount) %>%
  group_by(id,time,meSNP,sigcount) %>% 
  tally() %>% 
  pivot_wider(id_cols = c(id,time,sigcount), names_from = meSNP,names_glue = "meSNP_{meSNP}", values_from = n)

meSNPdf %>% 
  kable(., caption= "Number of genes from minimally expressed SNPs found in DEGs (sig and nonsig) broken down by drug,time, and significance") %>% 
  kable_paper("striped", full_width = FALSE) %>% 
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped"),font_size = 18) %>% 
  scroll_box(width = "80%", height = "400px")



reQTLdf <- toplistall %>% 
  mutate(id = as.factor(id)) %>%
  mutate(time=factor(time, levels=c("3_hours","24_hours"))) %>%
  mutate(sigcount = if_else(adj.P.Val < 0.05,'sig','notsig'))%>%
  mutate(reQTL = if_else( ENTREZID %in% knowles5$entrezgene_id, "y" , "no")) %>% 
  dplyr::select(ENTREZID,id,time,reQTL,sigcount) %>%
  group_by(id,time,reQTL,sigcount) %>% 
  tally() %>% 
  pivot_wider(id_cols = c(id,time,sigcount), names_from = reQTL,names_glue = "reQTL_{reQTL}", values_from = n)

reQTLdf %>% 
  kable(., caption= "Number of genes from response eQTLS found in DEGs (sig and nonsig) broken down by drug,time, and significance") %>% 
  kable_paper("striped", full_width = FALSE) %>% 
  kable_styling(full_width = FALSE, position = "left",bootstrap_options = c("striped"),font_size = 18) %>% 
  scroll_box(width = "80%", height = "400px")

```
###  chi square
```{r chi square test}

chi_squarek4 <- toplistall %>% 
  mutate(id = as.factor(id)) %>%
  dplyr::filter(id!="Trastuzumab") %>% 
  mutate(sigcount = if_else(adj.P.Val <0.05,'sig','notsig'))%>%
  mutate(meSNP=if_else(ENTREZID %in%knowles4$entrezgene_id,"y","no")) %>% 
  group_by(id,time) %>% 
  summarise(pvalue= chisq.test(meSNP, sigcount)$p.value) 

chi_squarek5 <- 
  toplistall %>% 
  mutate(id = as.factor(id)) %>%
  dplyr::filter(id!="Trastuzumab") %>% 
  mutate(sigcount = if_else(adj.P.Val <0.05,'sig','notsig'))%>%
  mutate(reQTL=if_else(ENTREZID %in%knowles5$entrezgene_id,"a","b")) %>% 
  group_by(id,time) %>% 
  summarise(pvalue= chisq.test(reQTL, sigcount)$p.value) 

chi_squarek4 %>% 
  kable(., caption= "meQTLs (mininmally expressed QTLS) chi-squared pvalues") %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "right",bootstrap_options = c("striped"),font_size = 18) %>% 
  scroll_box(width = "50%", height = "400px")

chi_squarek5 %>% 
  kable(., caption= "reQTLS in knowles data (reQTLS) chi-squared pvalues") %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "right",bootstrap_options = c("striped"),font_size = 18) %>% 
  scroll_box(width = "50%", height = "400px")

```










```{r abs FC }
toplistall %>% 
  filter(time=="24_hours") %>% 
  group_by(time, id) %>% 
  mutate(sigcount = if_else(adj.P.Val < 0.05,'sig','notsig'))%>%
  mutate(meSNP=if_else(ENTREZID %in%knowles4$entrezgene_id,"y","no")) %>% 
  mutate(reQTL=if_else(ENTREZID %in%knowles5$entrezgene_id,"a","b")) %>% 
  ggplot(., aes(x=id, y=abs(logFC)))+
  geom_boxplot(aes(fill=id))+
  ggpubr::fill_palette(palette =drug_palNoVeh)+
  guides(fill=guide_legend(title = "Treatment"))+
  facet_wrap(~time)+
  theme_bw()+
  xlab("")+
  ylab("abs |Log Fold Change|")+
  theme_bw()+
  facet_wrap(~time)+
  #ggtitle("All QTLs from all expressed genes (n=507)")+
  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 = element_text(size = 8, color = "black", angle = 0),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))
  
```





```{r Knowles comparison between sets}

Knowles_count <- 
  toplistall %>%
  filter(id!='Trastuzumab') %>% 
  mutate(id = as.factor(id)) %>%
  mutate(time=factor(time, levels=c("3_hours","24_hours"))) %>%
  group_by(time, id) %>%
  mutate(K4 = if_else(ENTREZID %in% knowles4$entrezgene_id,1,0))%>%
  mutate(K5 = if_else(ENTREZID %in% knowles5$entrezgene_id,1,0))%>%
  filter(adj.P.Val<0.05) %>%
  dplyr::summarize(n=n(), K4=sum(K4), K5=sum(K5)) %>% 
  as.tibble() %>% 
  dplyr::select(time,id,K4,K5) %>% rename("K4_y"='K4',"K5_y"='K5') %>% 
  mutate(time = case_match(time, '3_hours'~'3 hrs',
                           '24_hours'~'24 hrs',.default = time)) %>% 
  mutate(K4_n= 503-K4_y, K5_n=371-K5_y) %>% 
  pivot_longer(!c(time,id), names_to='QTL',values_to="gene_count") %>%
  separate(QTL,into=c("QTL_type",'group'),sep = '_') %>%  
  mutate(QTL_type =case_match(QTL_type, 'K4'~'meQTLs',
                         'K5'~'reQTLs',.default = QTL_type)) %>% 
  mutate(time=factor(time, levels=c("3 hrs","24 hrs"))) %>% 
  group_by(id,time,QTL_type) %>% 
  mutate(percent=gene_count/sum(gene_count)*100)
  
  ggplot(Knowles_count, aes(x=QTL_type,y=gene_count, group=group,  fill=group))+
    geom_col(position='fill')+
    guides(fill="none")+
    facet_wrap(time~id,nrow=2,ncol=4)+
    theme_classic()+
   scale_color_manual(values=drug_palNoVeh)+
    scale_fill_manual(values=drug_palNoVeh)+
       scale_y_continuous(label = scales::percent)

  ###This works
  chi_frame <- Knowles_count %>%
  pivot_wider(id_cols = c(time, id,QTL_type), names_from = group, values_from = gene_count) 

testDNR3chi <- matrix(c(6,13,497, 358), nrow=2, ncol=2, byrow=FALSE,
dimnames=list(c("k4", "k5"),c( "y", "n"))) 
    DNR_3chi <- chisq.test(testDNR3chi,correct=FALSE)$p.value
testDNR24chi <- matrix(c(199,187,304, 184), nrow=2, ncol=2, byrow=FALSE,
dimnames=list(c("k4", "k5"),c( "y", "n")))  
    DNR_24chi <- chisq.test(testDNR24chi,correct=FALSE)$p.value

testDOX3chi <- matrix(c(0,0,503, 371), nrow=2, ncol=2, byrow=FALSE,
dimnames=list(c("k4", "k5"),c( "y", "n")))  
    DOX_3chi <- chisq.test(testDOX3chi,correct=FALSE)$p.value
testDOX24chi <- matrix(c(184,180,318, 191), nrow=2, ncol=2, byrow=FALSE,
dimnames=list(c("k4", "k5"),c( "y", "n"))) 
    DOX_24chi <- chisq.test(testDOX24chi,correct=FALSE)$p.value

testEPI3chi <- matrix(c(4,4,499, 367), nrow=2, ncol=2, byrow=FALSE,
dimnames=list(c("k4", "k5"),c( "y", "n"))) 
    EPI_3chi <- chisq.test(testEPI3chi,correct=FALSE)$p.value
testEPI24chi <- matrix(c(172,176,331, 195), nrow=2, ncol=2, byrow=FALSE,
dimnames=list(c("k4", "k5"),c( "y", "n"))) 
    EPI_24chi <- chisq.test(testEPI24chi,correct=FALSE)$p.value
  
testMTX3chi <- matrix(c(1,0,502, 370), nrow=2, ncol=2, byrow=FALSE,
dimnames=list(c("k4", "k5"),c( "y", "n"))) 
    MTX_3chi <- chisq.test(testMTX3chi,correct=FALSE)$p.value
testMTX24chi <- matrix(c(30,40,473, 331), nrow=2, ncol=2, byrow=FALSE,
dimnames=list(c("k4", "k5"),c( "y", "n"))) 
    MTX_24chi <- chisq.test(testMTX24chi,correct=FALSE)$p.value
K4nK5chi <- data.frame(treatment=c('DNR_3','DNR_24','DOX_3','DOX_24','EPI_3','EPI_24','MTX_3','MTX_24'), chi_p.value=c(DNR_3chi,DNR_24chi,DOX_3chi,DOX_24chi,EPI_3chi,EPI_24chi,MTX_3chi,MTX_24chi))

K4nK5chi %>% kable(., caption= "Chi Square p. values between proportions of meQTLs and reQTLS by time and treatment") %>% 
  kable_paper("striped", full_width = TRUE) %>%  
  kable_styling(full_width = FALSE, position = "right",bootstrap_options = c("striped"),font_size = 18) %>% 
  scroll_box( height = "500px")
 
```

























## GWAS
```{r GOI analysis uploading, echo= FALSE,eval=FALSE}
GWAS_goi <- c('RARG', 'ITGB7', 'TNS2','ZNF740','SLC28A3','RMI1',
'FEDORA' ,'GDF5','FRS2','HDDC2','EEF1B2')

library(biomaRt)
ensembl <- useMart("ensembl", dataset="hsapiens_gene_ensembl")
my_chr <- c(1:22, 'M', 'X', 'Y')  
my_attributes <- c('entrezgene_id', 'ensembl_gene_id', 'hgnc_symbol')


GWAS_goi<- getBM(attributes=my_attributes,filters ='hgnc_symbol',
         values = GWAS_goi, mart = ensembl)
GWAS_goi<-GWAS_goi %>% distinct(entrezgene_id,.keep_all = TRUE) %>% add_row(entrezgene_id='124903732',ensembl_gene_id='ENSG00000260788', hgnc_symbol="RP11-298D21.1
")

  
write.csv(GWAS_goi,"output/GWAS_goi.csv")


``` 
```{r data sets loading, echo=FALSE, message=FALSE}
library(ComplexHeatmap)
library(tidyverse)
library(ggsignif)
library(biomaRt)
library(RColorBrewer)
library(cowplot)
# library(ggpubr)
library(scales)
library(sjmisc)
library(kableExtra)
library(broom)
library(ggstats)
toplistall.csv <- read.csv("output/toplistall.csv", row.names = 1)

drug_palNoVeh <- c("#8B006D","#DF707E","#F1B72B", "#3386DD","#707031","#41B333")
toplistall <- readRDS("data/toplistall.RDS")
knowles4all <-readRDS("output/knowles4.RDS")
knowles5all <-readRDS("output/knowles5.RDS")
DOXreQTLs <- readRDS("output/DOXreQTLs.RDS")
GWAS_goi <- read.csv("output/GWAS_goi.csv", row.names = 1)


toplist24hr <-   toplistall %>%
  mutate(id = as.factor(id)) %>%
  mutate(time=factor(time, levels=c("3_hours","24_hours"))) %>%
  filter(time=="24_hours")

toplist3hr <-   toplistall %>%
  mutate(id = as.factor(id)) %>%
  mutate(time=factor(time, levels=c("3_hours","24_hours"))) %>%
  filter(time=="3_hours")



  
  # guides(fill=guide_legend(title = "Treatment"))+
  # facet_wrap(~time)+
  # theme_bw()+
  # xlab("")+
  # ylab("|Log Fold Change|")+
  # theme_bw()+
  # ggtitle("|Log Fold| for all genes in NR set")+
  # theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
  #       axis.title = element_text(size = 15, color = "black"),
  #       axis.line = element_line(linewidth = 1.5),
  #       strip.background = element_rect(fill = "#C77CFF"),
  #       axis.text = element_text(size = 8, color = "black", angle = 15),         strip.text.x = element_text(size = 12, color = "black", face = "bold"))


  
```


```{r GWAS genes of interest}

##get the abs FC of all GOI
GWASabsFCsig <- 
  toplistall %>% 
  mutate(absFC=abs(logFC)) %>% 
  mutate(id = as.factor(id)) %>%
  filter(id !="Trastuzumab") %>% 
  mutate(time=factor(time, levels=c("3_hours","24_hours"))) %>%
  filter(ENTREZID %in% GWAS_goi$entrezgene_id) %>% 
  dplyr::select(ENTREZID ,time, id, absFC, adj.P.Val, SYMBOL) %>% 
  filter(adj.P.Val <0.05) %>% 
  mutate(id =case_match(id,'Daunorubicin'~'DNR','Doxorubicin'~'DOX','Epirubicin'~'EPI','Mitoxantrone'~'MTX', .default = id)) %>% 
  mutate(time =case_match(time,"3_hours"~'3_hr',"24_hours"~'24_hr',.default = time)) %>% 
 unite("trt_time",id:time,sep = '_') %>% 
  pivot_wider(id_cols=trt_time, names_from = SYMBOL, values_from =adj.P.Val)# %>% 
  
 
GWASabsFC <- toplistall %>% 
  mutate(absFC=abs(logFC)) %>% 
  mutate(id = as.factor(id)) %>%
  filter(id !="Trastuzumab") %>% 
  # mutate(time=factor(time, levels=c("3_hours","24_hours"))) %>%
  filter(ENTREZID %in% GWAS_goi$entrezgene_id) %>% 
  dplyr::select(SYMBOL ,time, id, absFC) %>% 
  mutate(id =case_match(id,'Daunorubicin'~'DNR','Doxorubicin'~'DOX','Epirubicin'~'EPI','Mitoxantrone'~'MTX', .default = id)) %>% 
  mutate(time =case_match(time,"3_hours"~'3_hr',"24_hours"~'24_hr',.default = time)) %>% 
 unite("trt_time",id:time,sep = '_') %>% 
  pivot_wider(id_cols=trt_time, names_from = SYMBOL, values_from = absFC) %>% 
  column_to_rownames(var="trt_time")# %>%

Heatmap( t(as.matrix(GWASabsFC)), name = "Absolute\nFold Change\nvalues", 
         cluster_rows = FALSE, cluster_columns = FALSE,
         # row_order = c(1,5,2,6,3,7,4,8),
         row_order = c('RARG',  
                          'TNS2', 
                          'ZNF740',
                          'SLC28A3',
                          'RMI1',
                          'FRS2', 
                          'HDDC2','EEF1B2'))


GWASabsFCsig %>% 
  mutate_at(.vars = 2:8, .funs= scientific_format()) %>% 
  kable(., caption= "Significance of GWAS sets") %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",font_size = 18) %>% 
  scroll_box(width = "80%", height = "400px")



```



```{r   stacked barplot}
library(data.table)
DEG_cormotif <- readRDS("data/DEG_cormotif.RDS")
list2env(DEG_cormotif,envir=.GlobalEnv)
backGL <- read.csv("data/backGL.txt")
siglist <- readRDS("data/siglist.RDS")
list2env(siglist,envir=.GlobalEnv)

k4test <- knowles4all  %>% 
  distinct(entrezgene_id,.keep_all = TRUE) %>% 
  mutate(DNR_24=if_else(entrezgene_id %in% sigVDA24$ENTREZID,'y','no')) %>% 
  mutate(DOX_24=if_else(entrezgene_id %in% sigVDX24$ENTREZID,'y','no')) %>% 
  mutate(EPI_24=if_else(entrezgene_id %in% sigVEP24$ENTREZID,'y','no')) %>% 
  mutate(MTX_24=if_else(entrezgene_id %in% sigVMT24$ENTREZID,'y','no')) %>% 
  mutate(DNR_3=if_else(entrezgene_id %in% sigVDA3$ENTREZID,'y','no')) %>% 
  mutate(DOX_3=if_else(entrezgene_id %in% sigVDX3$ENTREZID,'y','no')) %>% 
  mutate(EPI_3=if_else(entrezgene_id %in% sigVEP3$ENTREZID,'y','no')) %>% 
  mutate(MTX_3=if_else(entrezgene_id %in% sigVMT3$ENTREZID,'y','no')) %>%
  dplyr::select(entrezgene_id,DNR_24,DNR_3,DOX_24,DOX_3,EPI_24,EPI_3,MTX_24,MTX_3) 
  
k5test <- knowles5all  %>% 
  distinct(entrezgene_id,.keep_all = TRUE) %>% 
  mutate(DNR_24=if_else(entrezgene_id %in% sigVDA24$ENTREZID,'y','no')) %>% 
  mutate(DOX_24=if_else(entrezgene_id %in% sigVDX24$ENTREZID,'y','no')) %>% 
  mutate(EPI_24=if_else(entrezgene_id %in% sigVEP24$ENTREZID,'y','no')) %>% 
  mutate(MTX_24=if_else(entrezgene_id %in% sigVMT24$ENTREZID,'y','no')) %>% 
  mutate(DNR_3=if_else(entrezgene_id %in% sigVDA3$ENTREZID,'y','no')) %>% 
  mutate(DOX_3=if_else(entrezgene_id %in% sigVDX3$ENTREZID,'y','no')) %>% 
  mutate(EPI_3=if_else(entrezgene_id %in% sigVEP3$ENTREZID,'y','no')) %>% 
  mutate(MTX_3=if_else(entrezgene_id %in% sigVMT3$ENTREZID,'y','no')) %>%
  dplyr::select(entrezgene_id,DNR_24,DNR_3,DOX_24,DOX_3,EPI_24,EPI_3,MTX_24,MTX_3) 
  


length(toplist24hr$ENTREZID[toplist24hr$id=='Daunorubicin'])


testdf <- data.frame(colnames(k4test)) %>% 
  slice(2:9) %>% 
  rename('trx'=colnames.k4test.) %>% 
  mutate(K4_y=c(length(k4test$entrezgene_id[k4test$DNR_24=='y']),
               length(k4test$entrezgene_id[k4test$DNR_3=='y']),
               length(k4test$entrezgene_id[k4test$DOX_24=='y']),
               length(k4test$entrezgene_id[k4test$DOX_3=='y']),
               length(k4test$entrezgene_id[k4test$EPI_24=='y']),
               length(k4test$entrezgene_id[k4test$EPI_3=='y']),
               length(k4test$entrezgene_id[k4test$MTX_24=='y']),
               length(k4test$entrezgene_id[k4test$MTX_3=='y'])),
         K4_n=c(length(k4test$entrezgene_id[k4test$DNR_24=='no']),
               length(k4test$entrezgene_id[k4test$DNR_3=='no']),
               length(k4test$entrezgene_id[k4test$DOX_24=='no']),
               length(k4test$entrezgene_id[k4test$DOX_3=='no']),
               length(k4test$entrezgene_id[k4test$EPI_24=='no']),
               length(k4test$entrezgene_id[k4test$EPI_3=='no']),
               length(k4test$entrezgene_id[k4test$MTX_24=='no']),
               length(k4test$entrezgene_id[k4test$MTX_3=='no'])),
         K5_y=c(length(k5test$entrezgene_id[k5test$DNR_24=='y']),
               length(k5test$entrezgene_id[k5test$DNR_3=='y']),
               length(k5test$entrezgene_id[k5test$DOX_24=='y']),
               length(k5test$entrezgene_id[k5test$DOX_3=='y']),
               length(k5test$entrezgene_id[k5test$EPI_24=='y']),
               length(k5test$entrezgene_id[k5test$EPI_3=='y']),
               length(k5test$entrezgene_id[k5test$MTX_24=='y']),
               length(k5test$entrezgene_id[k5test$MTX_3=='y'])),
         K5_n=c(length(k5test$entrezgene_id[k5test$DNR_24=='no']),
               length(k5test$entrezgene_id[k5test$DNR_3=='no']),
               length(k5test$entrezgene_id[k5test$DOX_24=='no']),
               length(k5test$entrezgene_id[k5test$DOX_3=='no']),
               length(k5test$entrezgene_id[k5test$EPI_24=='no']),
               length(k5test$entrezgene_id[k5test$EPI_3=='no']),
               length(k5test$entrezgene_id[k5test$MTX_24=='no']),
               length(k5test$entrezgene_id[k5test$MTX_3=='no'])))
 

Knowles_count <- testdf %>% 
  pivot_longer(!trx, names_to='QTL',values_to="gene_count") %>% 
  separate(trx,into=c("treatment",'time'),sep = '_') %>%
  separate(QTL,into=c("QTL_type",'group'),sep = '_') %>% 
  mutate(time= case_match(time, '3'~'3 hrs','24'~'24 hrs',.default = time)) %>% 
  mutate(time=factor(time, levels=c("3 hrs","24 hrs"))) %>% 
  mutate(QTL_type =case_match(QTL_type, 'K4'~'meSNPs','K5'~'reQTLs',.default = QTL_type)) %>% 
  group_by(treatment,time,QTL_type) %>% 
  mutate(percent=gene_count/sum(gene_count)*100)
  
  ggplot(Knowles_count, aes(x=QTL_type,y=gene_count, group=group, col=group, fill=group))+
    geom_col(position='fill')+
    guides(fill="none")+
    facet_wrap(treatment~time,nrow=2,ncol=4)+
    theme_classic()+
   scale_color_manual(values=drug_palNoVeh)+
    scale_fill_manual(values=drug_palNoVeh)+
       scale_y_continuous(label = scales::percent)
    
  
 

 
 
 
```



```{r doxeqtls  }
summaryDOXQTL %>% kable(., caption= "Number of DOXeqtl that are significant in each set") %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",font_size = 18) %>% 
  scroll_box(width = "80%", height = "400px")
nosigsumDOXQTL %>% kable(., caption= "Number of DOXeqtl that are significant in each set") %>% 
  kable_paper("striped", full_width = FALSE) %>%  
  kable_styling(full_width = FALSE, position = "left",font_size = 18) %>% 
  scroll_box(width = "80%", height = "400px")


```