---
title: "Pancancerdnds"
author: "Lucie_Gourmet"
date: "09/11/2021"
output: html_document
---

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

```{r}
library(tidyverse)
library(ggplot2)
library(dndscv)
library(ggpubr)
library(ggrepel)
rm(list = ls())
```

```{r function2runandplot}
#Luis
#load("~/Dropbox (Personal)/PI_position/Projects/Aim2_DriverLandscape_Copy20Dec2021/Final_list_escape_genes/annot_table_with_immune_category.RData")

#Lucie
load("~/Dropbox/Projects/DriverLandscape/Final_list_escape_genes/annot_table_with_immune_category.RData")

annot_table<-annot_table_with_immune_category

#Luis
refdb_run ="~/Dropbox (Personal)/PI_position/Projects/DriverLandscape/Data/RefCDS_human_GRCh38.p12.rda"



#Lucie
refdb_run ="~/Dropbox/Projects/DriverLandscape/Data/RefCDS_human_GRCh38.p12.rda"
all<-c("YES","NO")
run_sel_plot<-function(x=all){
  #Separating escape from non escaped
  df_x<-annot_table %>% filter(escapedall==x)
  
  df_esc_x<-df_x %>% filter(escapedall=="YES") %>% dplyr::select(sampleID,chr,pos,ref,mut) 
  df_non_x<-df_x %>% filter(escapedall=="NO")%>% dplyr::select(sampleID,chr,pos,ref,mut)
  
  df_esc_x_dndscv<-dndscv(df_esc_x, refdb = refdb_run ,cv=NULL)
  df_non_x_dndscv<-dndscv(df_non_x, refdb = refdb_run ,cv=NULL)
  
  df_x_dndscv<-list(df_esc_x_dndscv,df_non_x_dndscv)
  
  return(df_x_dndscv)
}
###List 1 -> ESCAPED
###List 2 -> NONESCAPED

run_plot_cor2phen<-function(x){
  blacklist<-c("B2M","CALR","CANX", "CD4", "CD74","CD8A","CD8B","CIITA","CREB1","CTSB","CTSL", "CTSS", "ERAP1","ERAP2", "FAS", "HLA-A","HLA-B", "HLA-C", "HLA-DMA", "HLA-DMB", "HLA-DOA", "HLA-DOB","HLA-DPA1", "HLA-DPB1","HLA-DQA1","HLA-DQA2","HLA-DQB1","HLA-DRA","HLA-DRB1","HLA-DRB3","HLA-DRB4","HLA-DRB5","HLA-E","HLA-F","HLA-G","HSP90AA1","HSP90AB1","HSPA1A","HSPA1B","HSPA1L","HSPA2","HSPA4","HSPA5","HSPA6","HSPA8","HSPBP1","IFI30","IFNG","IRF1","KIR2DL1","KIR2DL2","KIR2DL3","KIR2DL4","KIR2DS1","KIR2DS2","KIR2DS4","KIR2DS5","KIR3DL1","KIR3DL2","KIR3DL3","KLRC1","KLRC2","KLRC3","KLRC4","KLRD1","LGMN","MEX3B","NFYA","NFYB","NFYC","PDIA3","PSMA7","PSMB10","PSMB11","PSMB6","PSMB8","PSMB9","PSME1","PSME2","PSME3","PSMF1","RFX5","RFXANK","RFXAP","TAP1","TAP2","TAPBP","TNF", "TTN", "TTN")
df_1<-x[[1]]$sel_cv
df_2<-x[[2]]$sel_cv
df_1 <- df_1 %>% filter(!str_detect(gene_name, "^OR|^CS"))
df_2 <- df_2 %>% filter(!str_detect(gene_name, "^OR|^CS"))
dfboth<-left_join(df_1, df_2,by = c("gene_name"),suffix=c(".escapedall",".nonescapedall"))
dfboth <- dfboth %>% mutate(logdndsescapedall=log2(wmis_cv.escapedall),logdndsnonescapedall=log2(wmis_cv.nonescapedall),logMuts=log10(n_syn.escapedall+n_syn.nonescapedall+1))
  dfboth$Significant<-""
dfboth$Significant <- ifelse((dfboth$qmis_cv.nonescapedall < 0.1 & dfboth$qmis_cv.escapedall < 0.1), "SigBoth", ifelse(dfboth$qmis_cv.nonescapedall < 0.1,"Signonescapedall", ifelse( dfboth$qmis_cv.escapedall < 0.1,"Sigescapedall","NotSig")))
dfboth.filt<-dfboth %>% filter(dfboth$qmis_cv.nonescapedall <0.1|dfboth$qmis_cv.escapedall <0.1)
dfboth.filt<-subset(dfboth.filt, !(gene_name %in% blacklist))

p<-ggscatter(dfboth.filt,
             y = "logdndsnonescapedall", x = "logdndsescapedall",
             color="Significant",fill= "Significant",
             #add = "reg.line", 
             shape = 21, size = 3,  # Add regressin line
             add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
             conf.int = FALSE, # Add confidence interval
             cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
             cor.coeff.args = list(method = "pearson", label.sep = "\n"),
             
             label = "gene_name", repel = TRUE)+
  geom_hline(yintercept=0, linetype="dashed",
             color = "red", size=0.3)+
  geom_vline(xintercept=0, linetype="dashed",
             color = "red", size=0.3) +
  geom_abline(intercept=0,slope=1, linetype="twodash", color="grey5", size=0.3) 

return(p)
}

run_volcano_all<-function(x){
  
  blacklist<-c("B2M","CALR","CANX", "CD4", "CD74","CD8A","CD8B","CIITA","CREB1","CTSB","CTSL", "CTSS", "ERAP1","ERAP2", "FAS", "HLA-A","HLA-B", "HLA-C", "HLA-DMA", "HLA-DMB", "HLA-DOA", "HLA-DOB","HLA-DPA1", "HLA-DPB1","HLA-DQA1","HLA-DQA2","HLA-DQB1","HLA-DRA","HLA-DRB1","HLA-DRB3","HLA-DRB4","HLA-DRB5","HLA-E","HLA-F","HLA-G","HSP90AA1","HSP90AB1","HSPA1A","HSPA1B","HSPA1L","HSPA2","HSPA4","HSPA5","HSPA6","HSPA8","HSPBP1","IFI30","IFNG","IRF1","KIR2DL1","KIR2DL2","KIR2DL3","KIR2DL4","KIR2DS1","KIR2DS2","KIR2DS4","KIR2DS5","KIR3DL1","KIR3DL2","KIR3DL3","KLRC1","KLRC2","KLRC3","KLRC4","KLRD1","LGMN","MEX3B","NFYA","NFYB","NFYC","PDIA3","PSMA7","PSMB10","PSMB11","PSMB6","PSMB8","PSMB9","PSME1","PSME2","PSME3","PSMF1","RFX5","RFXANK","RFXAP","TAP1","TAP2","TAPBP","TNF", "TTN")
  
  df_1<-x[[1]]$sel_cv
  df_2<-x[[2]]$sel_cv
  
  df_1 <- df_1 %>% filter(!str_detect(gene_name, "^OR|^CS"))
  df_2 <- df_2 %>% filter(!str_detect(gene_name, "^OR|^CS"))
  
  
  #####trick to plot genes with pvalue equals 0
  minp1<-min(df_1 %>% filter(pallsubs_cv>0) %>% pull(pallsubs_cv))
  df_1<-df_1 %>% mutate(ifelse(pallsubs_cv==0,minp1-runif(1)*minp1,pallsubs_cv))
  
  minp2<-min(df_2 %>% filter(pallsubs_cv>0) %>% pull(pallsubs_cv))
  df_2<-df_2 %>% mutate(ifelse(pallsubs_cv==0,minp2-runif(1)*minp2,pallsubs_cv))
  ######End of trick######
  
  
  
  df_1$log2FolddNdS<-log2(df_1$wallsubs_cv)
  df_2$log2FolddNdS<-log2(df_2$wallsubs_cv)
  
  df_1$Significant<-""
  df_1$Significant <- ifelse(df_1$qallsubs_cv < 0.1, paste("FDR <","0.1"), "Not Sig")
  df_1<-subset(df_1, !(gene_name %in% blacklist))
  
  
  p1 <- ggplot(df_1, aes(x = log2FolddNdS, y = -log10(pallsubs_cv))) +
    geom_point(aes(color = Significant)) +
    scale_color_manual(values = c("red", "grey")) +
    theme_bw(base_size = 12) + theme(legend.position = "bottom") +
    geom_text_repel(
      data = subset(df_1, qallsubs_cv < 0.1),
      aes(label = gene_name),
      size = 3, max.overlaps = Inf,
      box.padding = unit(0.35, "lines"),
      point.padding = unit(0.3, "lines")
    )   
  #p1
  
  df_2$Significant<-""
  df_2$Significant <- ifelse(df_2$qallsubs_cv < 0.1, paste("FDR <","0.1"), "Not Sig")
  df_2<-subset(df_2, !(gene_name %in% blacklist))
  
  p2 <- ggplot(df_2, aes(x = log2FolddNdS, y = -log10(pallsubs_cv))) +
    geom_point(aes(color = Significant)) +
    scale_color_manual(values = c("red", "grey")) +
    theme_bw(base_size = 12) + theme(legend.position = "bottom") +
    geom_text_repel(
      data = subset(df_2, qallsubs_cv < 0.1),
      aes(label = gene_name),
      size = 3,max.overlaps = Inf,
      box.padding = unit(0.35, "lines"),
      point.padding = unit(0.3, "lines")
    )   
  #p2
  
  p3<-ggarrange(p1,p2,labels="AUTO",common.legend = T)
  return(p3)
}


run_volcano<-function(x){
  
  blacklist<-c("B2M","CALR","CANX", "CD4", "CD74","CD8A","CD8B","CIITA","CREB1","CTSB","CTSL", "CTSS", "ERAP1","ERAP2", "FAS", "HLA-A","HLA-B", "HLA-C", "HLA-DMA", "HLA-DMB", "HLA-DOA", "HLA-DOB","HLA-DPA1", "HLA-DPB1","HLA-DQA1","HLA-DQA2","HLA-DQB1","HLA-DRA","HLA-DRB1","HLA-DRB3","HLA-DRB4","HLA-DRB5","HLA-E","HLA-F","HLA-G","HSP90AA1","HSP90AB1","HSPA1A","HSPA1B","HSPA1L","HSPA2","HSPA4","HSPA5","HSPA6","HSPA8","HSPBP1","IFI30","IFNG","IRF1","KIR2DL1","KIR2DL2","KIR2DL3","KIR2DL4","KIR2DS1","KIR2DS2","KIR2DS4","KIR2DS5","KIR3DL1","KIR3DL2","KIR3DL3","KLRC1","KLRC2","KLRC3","KLRC4","KLRD1","LGMN","MEX3B","NFYA","NFYB","NFYC","PDIA3","PSMA7","PSMB10","PSMB11","PSMB6","PSMB8","PSMB9","PSME1","PSME2","PSME3","PSMF1","RFX5","RFXANK","RFXAP","TAP1","TAP2","TAPBP","TNF", "TTN")
  
  df_1<-x[[1]]$sel_cv
  df_2<-x[[2]]$sel_cv
  
df_1 <- df_1 %>% filter(!str_detect(gene_name, "^OR|^CS"))
df_2 <- df_2 %>% filter(!str_detect(gene_name, "^OR|^CS"))


#####trick to plot genes with pvalue equals 0
  minp1<-min(df_1 %>% filter(pmis_cv>0) %>% pull(pmis_cv))
  df_1<-df_1 %>% mutate(ifelse(pmis_cv==0,minp1-runif(1)*minp1,pmis_cv))
  
  minp2<-min(df_2 %>% filter(pmis_cv>0) %>% pull(pmis_cv))
  df_2<-df_2 %>% mutate(ifelse(pmis_cv==0,minp2-runif(1)*minp2,pmis_cv))
  ######End of trick######
  
  
  
  df_1$log2FolddNdS<-log2(df_1$wmis_cv)
  df_2$log2FolddNdS<-log2(df_2$wmis_cv)
  
  df_1$Significant<-""
  df_1$Significant <- ifelse(df_1$qmis_cv < 0.1, paste("FDR <","0.1"), "Not Sig")
  df_1<-subset(df_1, !(gene_name %in% blacklist))

  
  p1 <- ggplot(df_1, aes(x = log2FolddNdS, y = -log10(pmis_cv))) +
  geom_point(aes(color = Significant)) +
  scale_color_manual(values = c("red", "grey")) +
  theme_bw(base_size = 12) + theme(legend.position = "bottom") +
  geom_text_repel(
    data = subset(df_1, qmis_cv < 0.1),
    aes(label = gene_name),
    size = 3, max.overlaps = Inf,
    box.padding = unit(0.35, "lines"),
    point.padding = unit(0.3, "lines")
  )   
  #p1
  
  df_2$Significant<-""
  df_2$Significant <- ifelse(df_2$qmis_cv < 0.1, paste("FDR <","0.1"), "Not Sig")
  df_2<-subset(df_2, !(gene_name %in% blacklist))
  
  p2 <- ggplot(df_2, aes(x = log2FolddNdS, y = -log10(pmis_cv))) +
  geom_point(aes(color = Significant)) +
  scale_color_manual(values = c("red", "grey")) +
  theme_bw(base_size = 12) + theme(legend.position = "bottom") +
  geom_text_repel(
    data = subset(df_2, qmis_cv < 0.1),
    aes(label = gene_name),
    size = 3,max.overlaps = Inf,
    box.padding = unit(0.35, "lines"),
    point.padding = unit(0.3, "lines")
  )   
  #p2
  
  p3<-ggarrange(p1,p2,labels="AUTO",common.legend = T)
  return(p3)
}

run_volcano_trunc<-function(x){
   blacklist<-c("B2M","CALR","CANX", "CD4", "CD74","CD8A","CD8B","CIITA","CREB1","CTSB","CTSL", "CTSS", "ERAP1","ERAP2", "FAS", "HLA-A","HLA-B", "HLA-C", "HLA-DMA", "HLA-DMB", "HLA-DOA", "HLA-DOB","HLA-DPA1", "HLA-DPB1","HLA-DQA1","HLA-DQA2","HLA-DQB1","HLA-DRA","HLA-DRB1","HLA-DRB3","HLA-DRB4","HLA-DRB5","HLA-E","HLA-F","HLA-G","HSP90AA1","HSP90AB1","HSPA1A","HSPA1B","HSPA1L","HSPA2","HSPA4","HSPA5","HSPA6","HSPA8","HSPBP1","IFI30","IFNG","IRF1","KIR2DL1","KIR2DL2","KIR2DL3","KIR2DL4","KIR2DS1","KIR2DS2","KIR2DS4","KIR2DS5","KIR3DL1","KIR3DL2","KIR3DL3","KLRC1","KLRC2","KLRC3","KLRC4","KLRD1","LGMN","MEX3B","NFYA","NFYB","NFYC","PDIA3","PSMA7","PSMB10","PSMB11","PSMB6","PSMB8","PSMB9","PSME1","PSME2","PSME3","PSMF1","RFX5","RFXANK","RFXAP","TAP1","TAP2","TAPBP","TNF", "TTN")
  df_1<-x[[1]]$sel_cv
  df_2<-x[[2]]$sel_cv
  df_1 <- df_1 %>% filter(!str_detect(gene_name, "^OR|^CS"))
df_2 <- df_2 %>% filter(!str_detect(gene_name, "^OR|^CS"))
  #####trick to plot genes with pvalue equals 0
  minp1<-min(df_1 %>% filter(ptrunc_cv>0) %>% pull(ptrunc_cv))
  df_1<-df_1 %>% mutate(ifelse(ptrunc_cv==0,minp1-runif(1)*minp1,ptrunc_cv))
  
  minp2<-min(df_2 %>% filter(ptrunc_cv>0) %>% pull(pmis_cv))
  df_2<-df_2 %>% mutate(ifelse(ptrunc_cv==0,minp2-runif(1)*minp2,ptrunc_cv))
  ######End of trick######
  
  df_1$log2FolddNdS<-log2(df_1$wnon_cv)
  df_2$log2FolddNdS<-log2(df_2$wnon_cv)
  
  df_1$Significant<-""
  df_1$Significant <- ifelse(df_1$qtrunc_cv < 0.1, paste("FDR <","0.1"), "Not Sig")
  df_1<-subset(df_1, !(gene_name %in% blacklist))
  p1 <- ggplot(df_1, aes(x = log2FolddNdS, y = -log10(ptrunc_cv))) +
  geom_point(aes(color = Significant)) +
  scale_color_manual(values = c("red", "grey")) +
  theme_bw(base_size = 12) + theme(legend.position = "bottom") +
  geom_text_repel(
    data = subset(df_1, qtrunc_cv < 0.1),
    aes(label = gene_name),
    size = 3,max.overlaps = Inf,
    box.padding = unit(0.35, "lines"),
    point.padding = unit(0.3, "lines")
  )   
  #p1
  
  df_2$Significant<-""
  df_2$Significant <- ifelse(df_2$qtrunc_cv < 0.1, paste("FDR <","0.1"), "Not Sig")
  df_2<-subset(df_2, !(gene_name %in% blacklist))
  p2 <- ggplot(df_2, aes(x = log2FolddNdS, y = -log10(ptrunc_cv))) +
  geom_point(aes(color = Significant)) +
  scale_color_manual(values = c("red", "grey")) +
  theme_bw(base_size = 12) + theme(legend.position = "bottom") +
  geom_text_repel(
    data = subset(df_2, qtrunc_cv < 0.1),
    aes(label = gene_name),
    size = 3,max.overlaps = Inf,
    box.padding = unit(0.35, "lines"),
    point.padding = unit(0.3, "lines")
  )   
  #p2
  
  p3<-ggarrange(p1,p2,labels="AUTO",common.legend = T)
  return(p3)
}

run_plot_cor2phen_trunc<-function(x){
  blacklist<-c("B2M","CALR","CANX", "CD4", "CD74","CD8A","CD8B","CIITA","CREB1","CTSB","CTSL", "CTSS", "ERAP1","ERAP2", "FAS", "HLA-A","HLA-B", "HLA-C", "HLA-DMA", "HLA-DMB", "HLA-DOA", "HLA-DOB","HLA-DPA1", "HLA-DPB1","HLA-DQA1","HLA-DQA2","HLA-DQB1","HLA-DRA","HLA-DRB1","HLA-DRB3","HLA-DRB4","HLA-DRB5","HLA-E","HLA-F","HLA-G","HSP90AA1","HSP90AB1","HSPA1A","HSPA1B","HSPA1L","HSPA2","HSPA4","HSPA5","HSPA6","HSPA8","HSPBP1","IFI30","IFNG","IRF1","KIR2DL1","KIR2DL2","KIR2DL3","KIR2DL4","KIR2DS1","KIR2DS2","KIR2DS4","KIR2DS5","KIR3DL1","KIR3DL2","KIR3DL3","KLRC1","KLRC2","KLRC3","KLRC4","KLRD1","LGMN","MEX3B","NFYA","NFYB","NFYC","PDIA3","PSMA7","PSMB10","PSMB11","PSMB6","PSMB8","PSMB9","PSME1","PSME2","PSME3","PSMF1","RFX5","RFXANK","RFXAP","TAP1","TAP2","TAPBP","TNF", "TTN", "TTN")
  df_1<-x[[1]]$sel_cv
  df_2<-x[[2]]$sel_cv
  df_1 <- df_1 %>% filter(!str_detect(gene_name, "^OR|^CS"))
  df_2 <- df_2 %>% filter(!str_detect(gene_name, "^OR|^CS"))
  df_1 <- df_1 %>% filter(!gene_name %in% blacklist)
  df_2 <- df_2 %>% filter(!gene_name %in% blacklist)
  dfboth<-left_join(df_1, df_2,by = c("gene_name"),suffix=c(".escapedall",".nonescapedall"))
  dfboth <- dfboth %>% mutate(logdndsescapedall=log2(wnon_cv.escapedall),logdndsnonescapedall=log2(wnon_cv.nonescapedall),logMuts=log10(n_syn.escapedall+n_syn.nonescapedall+1))
  dfboth$Significant<-""
  dfboth$Significant <- ifelse((dfboth$qtrunc_cv.nonescapedall < 0.1 & dfboth$qtrunc_cv.escapedall < 0.1), "SigBoth", ifelse(dfboth$qtrunc_cv.nonescapedall < 0.1,"Signonescapedall", ifelse( dfboth$qtrunc_cv.escapedall < 0.1,"Sigescapedall","NotSig")))
  dfboth.filt<-dfboth %>% filter(dfboth$qtrunc_cv.nonescapedall <0.1|dfboth$qtrunc_cv.escapedall <0.1)
  dfboth.filt<-subset(dfboth.filt, !(gene_name %in% blacklist))
  
  p<-ggscatter(dfboth.filt,
               y = "logdndsnonescapedall", x = "logdndsescapedall",
               color="Significant",fill= "Significant",
               #add = "reg.line", 
               shape = 21, size = 3,  # Add regressin line
               add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
               conf.int = FALSE, # Add confidence interval
               cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
               cor.coeff.args = list(method = "pearson", label.sep = "\n"),
               
               label = "gene_name", repel = TRUE)+
    geom_hline(yintercept=0, linetype="dashed",
               color = "red", size=0.3)+
    geom_vline(xintercept=0, linetype="dashed",
               color = "red", size=0.3) +
    geom_abline(intercept=0,slope=1, linetype="twodash", color="grey5", size=0.3) 
  
  return(p)
}


```

```{r run4each}

#get pancancer correlation results without stratification

dndscv_result_pan<-dndscv(annot_table, refdb = refdb_run ,cv=NULL)
save(dndscv_result_pan, file="dndscv_result_pan.RData")


df_1<-dndscv_result_pan$sel_cv
df_1 <- df_1 %>% filter(!str_detect(gene_name, "^OR|^CS"))

dfboth <- df_1%>% mutate(Significant_mis=ifelse(qmis_cv < 0.1 , "Sig", "Notsig"), Significant_trunc=ifelse(qtrunc_cv < 0.1 , "Sig", "Notsig"))
dfboth.filt<-dfboth %>% filter(Significant_mis=="Sig" | Significant_trunc=="Sig")
dfboth.filt<-subset(dfboth.filt, !(gene_name %in% blacklist))

pancancer_without_strat_mis<-dfboth.filt %>% filter(Significant_mis== "Sig")
pancancer_without_strat_trunc<-dfboth.filt %>% filter(Significant_trunc== "Sig")

save(pancancer_without_strat_mis, file="pancancer_without_strat_mis.RData")
save(pancancer_without_strat_trunc, file="pancancer_without_strat_trunc.RData")



#get pancancer correlation results with stratification
run_sel_plot<-function(){
  #Separating escape from non escaped
  
  df_esc_x<-df_x %>% filter(escapedall=="YES") %>% dplyr::select(sampleID,chr,pos,ref,mut) 
  df_non_x<-df_x %>% filter(escapedall=="NO")%>% dplyr::select(sampleID,chr,pos,ref,mut)
  
  df_esc_x_dndscv<-dndscv(df_esc_x, refdb = refdb_run ,cv=NULL)
  df_non_x_dndscv<-dndscv(df_non_x, refdb = refdb_run ,cv=NULL)
  
  df_x_dndscv<-list(df_esc_x_dndscv,df_non_x_dndscv)
  
  return(df_x_dndscv)
}



pancancer_results<-run_sel_plot()
save(pancancer_results, file="pancancer_results.RData")
#for missense mut
p1<-run_plot_cor2phen(pancancer_results)
 blacklist<-c("B2M","CALR","CANX", "CD4", "CD74","CD8A","CD8B","CIITA","CREB1","CTSB","CTSL", "CTSS", "ERAP1","ERAP2", "FAS", "HLA-A","HLA-B", "HLA-C", "HLA-DMA", "HLA-DMB", "HLA-DOA", "HLA-DOB","HLA-DPA1", "HLA-DPB1","HLA-DQA1","HLA-DQA2","HLA-DQB1","HLA-DRA","HLA-DRB1","HLA-DRB3","HLA-DRB4","HLA-DRB5","HLA-E","HLA-F","HLA-G","HSP90AA1","HSP90AB1","HSPA1A","HSPA1B","HSPA1L","HSPA2","HSPA4","HSPA5","HSPA6","HSPA8","HSPBP1","IFI30","IFNG","IRF1","KIR2DL1","KIR2DL2","KIR2DL3","KIR2DL4","KIR2DS1","KIR2DS2","KIR2DS4","KIR2DS5","KIR3DL1","KIR3DL2","KIR3DL3","KLRC1","KLRC2","KLRC3","KLRC4","KLRD1","LGMN","MEX3B","NFYA","NFYB","NFYC","PDIA3","PSMA7","PSMB10","PSMB11","PSMB6","PSMB8","PSMB9","PSME1","PSME2","PSME3","PSMF1","RFX5","RFXANK","RFXAP","TAP1","TAP2","TAPBP","TNF", "TTN", "TTN")
df_1<-pancancer_results[[1]]$sel_cv
df_2<-pancancer_results[[2]]$sel_cv
df_1 <- df_1 %>% filter(!str_detect(gene_name, "^OR|^CS"))
df_2 <- df_2 %>% filter(!str_detect(gene_name, "^OR|^CS"))
dfboth<-left_join(df_1, df_2,by = c("gene_name"),suffix=c(".escapedall",".nonescapedall"))
dfboth <- dfboth %>% mutate(logdndsescapedall=log2(wmis_cv.escapedall),logdndsnonescapedall=log2(wmis_cv.nonescapedall),logMuts=log10(n_syn.escapedall+n_syn.nonescapedall+1))
  dfboth$Significant<-""
dfboth$Significant <- ifelse((dfboth$qmis_cv.nonescapedall < 0.1 & dfboth$qmis_cv.escapedall < 0.1), "SigBoth", ifelse(dfboth$qmis_cv.nonescapedall < 0.1,"Signonescapedall", ifelse( dfboth$qmis_cv.escapedall < 0.1,"Sigescapedall","NotSig")))
dfboth.filt<-dfboth %>% filter(dfboth$qmis_cv.nonescapedall <0.1|dfboth$qmis_cv.escapedall <0.1)
dfboth.filt<-subset(dfboth.filt, !(gene_name %in% blacklist))

p<-ggscatter(dfboth.filt,
             y = "logdndsnonescapedall", x = "logdndsescapedall",
             color="Significant",fill= "Significant",
             #add = "reg.line", 
             shape = 21, size = 3,  # Add regressin line
             add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
             conf.int = FALSE, # Add confidence interval
             cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
             cor.coeff.args = list(method = "pearson", label.sep = "\n"),
             
             label = "gene_name", repel = TRUE)+
  geom_hline(yintercept=0, linetype="dashed",
             color = "red", size=0.3)+
  geom_vline(xintercept=0, linetype="dashed",
             color = "red", size=0.3) +
  geom_abline(intercept=0,slope=1, linetype="twodash", color="grey5", size=0.3) 
ggsave("pancancer_results_final.pdf",plot = p1, units = "cm",height = 15, width = 25)



library(openxlsx)
pancancer_results_stratification_mis<-dfboth.filt
write.xlsx(pancancer_results_stratification_mis, 'pancancer_results_stratification_mis.xlsx')
##################################count genes in escaped + and -
esc_all<-c("SigBoth", "Sigescapedall")
esc_specific<-c("Sigescapedall")

nonesc_all<-c("SigBoth", "Signonescapedall")
nonesc_specific<-c("Signonescapedall")

df_esc_all<-pancancer_results_stratification_mis%>% filter(Significant %in% esc_all)
df_esc_specific<-pancancer_results_stratification_mis%>% filter(Significant %in% esc_specific)


df_nonesc_all<-pancancer_results_stratification_mis%>% filter(Significant %in% nonesc_all)
df_nonesc_specific<-pancancer_results_stratification_mis%>% filter(Significant %in% nonesc_specific)
##################################c

p_cor_all<-run_plot_cor2phen(pancancer_results)
ggsave("Fig2E_Manuscript_pancancer_correlation.pdf",plot = p_cor_all, units = "cm",height = 15, width = 15)
p_volcano_mis<-run_volcano(pancancer_results)
ggsave("Fig2CandD_Pancancer_volcano_mis_cutoff.pdf",plot = p_volcano_mis, units = "cm",height = 15, width = 25)


#for truncating mutations
trunc<-run_volcano_trunc(pancancer_results)
ggsave("3_SF_Pancancer_volcano_trunc_cutoff.pdf",plot = trunc, units = "cm",height = 15, width = 25)

corr_trunc<-run_plot_cor2phen_trunc(pancancer_results)
ggsave("SF_Pancancer_pancancer_correlation__trunc.pdf",plot = corr_trunc, units = "cm",height = 15, width = 25)





############# Venn diagram for trunc

all_results<-dndscv_result_pan$sel_cv
escape_results<-pancancer_results[[1]]$sel_cv
nonescape_results<-pancancer_results[[2]]$sel_cv


df_non_esc<-full_join(nonescape_results, escape_results, by=c("gene_name"), suffix=c("_non", "_esc"))

df_non_esc_all<-full_join(df_non_esc, all_results, by=c("gene_name"), suffix=c("","_all"))
df_non_esc_all_filt<-df_non_esc_all%>%filter(!gene_name %in% blacklist)%>%filter(!str_detect(gene_name, "^OR|^CS"))
sig_non<-df_non_esc_all_filt%>% filter(qtrunc_cv_non<0.1)%>%pull(gene_name)
sig_esc<-df_non_esc_all_filt%>% filter(qtrunc_cv_esc<0.1)%>%pull(gene_name)
sig_all<-df_non_esc_all_filt%>% filter(qtrunc_cv<0.1)%>%pull(gene_name)

sig_esc_and_all<-df_non_esc_all_filt%>% filter(qtrunc_cv_esc<0.1 & qtrunc_cv<0.1)%>%pull(gene_name)
sig_non_and_all<-df_non_esc_all_filt%>% filter(qtrunc_cv_non<0.1 & qtrunc_cv<0.1)%>%pull(gene_name)
sig_non_and_esc<-df_non_esc_all_filt%>% filter(qtrunc_cv_non<0.1 & qtrunc_cv_esc<0.1)%>%pull(gene_name)
sig_non_and_esc_all<-df_non_esc_all_filt%>% filter(qtrunc_cv_non<0.1 & qtrunc_cv_esc<0.1 & qtrunc_cv<0.1)%>%pull(gene_name)


library("VennDiagram")

# move to new plotting page
grid.newpage()

# create Venn diagram with three sets
p<-draw.triple.venn(area1=41, area2=99, area3=64, 
                    n12=37, n23=59, n13=31, n123=31, 
                    category=c("Escape","All","Nonescape"),
                    col="Black",fill=c("Red","Yellow","Turquoise"))
grid.draw(p)
ggsave("Venn_diagram_drivers_trunc.pdf",plot = p_venn, units = "cm",height = 20, width = 20)



















##################################count genes in escaped + and -
blacklist<-c("B2M","CALR","CANX", "CD4", "CD74","CD8A","CD8B","CIITA","CREB1","CTSB","CTSL", "CTSS", "ERAP1","ERAP2", "FAS", "HLA-A","HLA-B", "HLA-C", "HLA-DMA", "HLA-DMB", "HLA-DOA", "HLA-DOB","HLA-DPA1", "HLA-DPB1","HLA-DQA1","HLA-DQA2","HLA-DQB1","HLA-DRA","HLA-DRB1","HLA-DRB3","HLA-DRB4","HLA-DRB5","HLA-E","HLA-F","HLA-G","HSP90AA1","HSP90AB1","HSPA1A","HSPA1B","HSPA1L","HSPA2","HSPA4","HSPA5","HSPA6","HSPA8","HSPBP1","IFI30","IFNG","IRF1","KIR2DL1","KIR2DL2","KIR2DL3","KIR2DL4","KIR2DS1","KIR2DS2","KIR2DS4","KIR2DS5","KIR3DL1","KIR3DL2","KIR3DL3","KLRC1","KLRC2","KLRC3","KLRC4","KLRD1","LGMN","MEX3B","NFYA","NFYB","NFYC","PDIA3","PSMA7","PSMB10","PSMB11","PSMB6","PSMB8","PSMB9","PSME1","PSME2","PSME3","PSMF1","RFX5","RFXANK","RFXAP","TAP1","TAP2","TAPBP","TNF", "TTN")
df_1<-pancancer_results[[1]]$sel_cv
df_2<-pancancer_results[[2]]$sel_cv
df_1 <- df_1 %>% filter(!str_detect(gene_name, "^OR|^CS"))
df_2 <- df_2 %>% filter(!str_detect(gene_name, "^OR|^CS"))
dfboth<-left_join(df_1, df_2,by = c("gene_name"),suffix=c(".escapedall",".nonescapedall"))
dfboth <- dfboth %>% mutate(logdndsescapedall=log2(wnon_cv.escapedall),logdndsnonescapedall=log2(wnon_cv.nonescapedall),logMuts=log10(n_syn.escapedall+n_syn.nonescapedall+1))
dfboth$Significant<-""
dfboth$Significant <- ifelse((dfboth$qtrunc_cv.nonescapedall < 0.1 & dfboth$qtrunc_cv.escapedall < 0.1), "SigBoth", ifelse(dfboth$qtrunc_cv.nonescapedall < 0.1,"Signonescapedall", ifelse( dfboth$qtrunc_cv.escapedall < 0.1,"Sigescapedall","NotSig")))
dfboth.filt.trunc<-dfboth %>% filter(dfboth$qtrunc_cv.nonescapedall <0.1|dfboth$qtrunc_cv.escapedall <0.1)
dfboth.filt.trunc<-subset(dfboth.filt.trunc, !(gene_name %in% blacklist))

esc_all<-c("SigBoth", "Sigescapedall")
esc_specific<-c("Sigescapedall")

nonesc_all<-c("SigBoth", "Signonescapedall")
nonesc_specific<-c("Signonescapedall")

df_esc_all<-dfboth.filt.trunc%>% filter(Significant %in% esc_all)
df_esc_specific<-dfboth.filt.trunc%>% filter(Significant %in% esc_specific)


df_nonesc_all<-dfboth.filt.trunc%>% filter(Significant %in% nonesc_all)
df_nonesc_specific<-dfboth.filt.trunc%>% filter(Significant %in% nonesc_specific)


##################################count new driver genes
new_driver<-pancancer_results_stratification_mis%>% filter(Known_driver=="No")
new_driver_non<-pancancer_results_stratification_mis%>% filter(Known_driver=="No",Significant %in% nonesc_all)
new_driver_esc<-pancancer_results_stratification_mis%>% filter(Known_driver=="No",Significant %in% esc_all)




#df_all<-run_sel_plot(x = all)
#save(df_all, file="df_all.RData")

#df_1<-df_all[[1]]$sel_cv
#df_2<-df_all[[2]]$sel_cv
#df_1 <- df_1 %>% filter(!str_detect(gene_name, "^OR|^CS"))
#df_2 <- df_2 %>% filter(!str_detect(gene_name, "^OR|^CS"))
#dfboth<-left_join(df_1, df_2,by = c("gene_name"),suffix=c(".escapedall",".nonescapedall"))
#dfboth <- dfboth %>% mutate(logdndsescapedall=log2(wmis_cv.escapedall),logdndsnonescapedall=log2(wmis_cv.nonescapedall),logMuts=log10(n_syn.escapedall+n_syn.nonescapedall+1))
#  dfboth$Significant<-""
#dfboth$Significant <- ifelse((dfboth$qmis_cv.nonescapedall < 0.1 & dfboth$qmis_cv.escapedall < 0.1), "SigBoth", ifelse(dfboth$qmis_cv.nonescapedall < 0.1,"Signonescapedall", ifelse( dfboth$qmis_cv.escapedall < 0.1,"Sigescapedall","NotSig")))
#dfboth.filt<-dfboth %>% filter(dfboth$qmis_cv.nonescapedall <0.1|dfboth$qmis_cv.escapedall <0.1)
#blacklist<-c("B2M","CALR","CANX", "CD4", "CD74","CD8A","CD8B","CIITA","CREB1","CTSB","CTSL", "CTSS", "ERAP1","ERAP2", "FAS", "HLA-A","HLA-B", "HLA-C", "HLA-DMA", "HLA-DMB", "HLA-DOA", "HLA-DOB","HLA-DPA1", "HLA-DPB1","HLA-DQA1","HLA-DQA2","HLA-DQB1","HLA-DRA","HLA-DRB1","HLA-DRB3","HLA-DRB4","HLA-DRB5","HLA-E","HLA-F","HLA-G","HSP90AA1","HSP90AB1","HSPA1A","HSPA1B","HSPA1L","HSPA2","HSPA4","HSPA5","HSPA6","HSPA8","HSPBP1","IFI30","IFNG","IRF1","KIR2DL1","KIR2DL2","KIR2DL3","KIR2DL4","KIR2DS1","KIR2DS2","KIR2DS4","KIR2DS5","KIR3DL1","KIR3DL2","KIR3DL3","KLRC1","KLRC2","KLRC3","KLRC4","KLRD1","LGMN","MEX3B","NFYA","NFYB","NFYC","PDIA3","PSMA7","PSMB10","PSMB11","PSMB6","PSMB8","PSMB9","PSME1","PSME2","PSME3","PSMF1","RFX5","RFXANK","RFXAP","TAP1","TAP2","TAPBP","TNF", "TTN")
#remove genes from stratification
#df.filter<-subset(dfboth.filt, !(gene_name %in% blacklist))
#library(openxlsx)
#write.xlsx(df.filter, 'df.filter.xlsx')

#return(p)
#n<-unique(df_all[[1]]$sel_cv)
#df_all<-df_all %>% subset(df_all[[1]]$sel_cv, !(gene_name %in% blacklist))
#p_cor_all<-run_plot_cor2phen(df_all)
#ggsave("Fig2E_Manuscript_pancancer_correlation.pdf",plot = p_cor_all, units = "cm",height = 15, width = 15)
#p_volcano_mis<-run_volcano(df_all)
#trunc<-run_volcano_trunc(df_all)

#ggsave("Fig2CandD_Pancancer_volcano_mis_cutoff.pdf",plot = p_volcano_mis, units = "cm",height = 15, width = 25)
#ggsave("3_SF_Pancancer_volcano_trunc_cutoff.pdf",plot = trunc, units = "cm",height = 15, width = 25)





```

## R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>.

When you click the **Knit** button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

```{r cars}
blacklist<-c("B2M","CALR","CANX", "CD4", "CD74","CD8A","CD8B","CIITA","CREB1","CTSB","CTSL", "CTSS", "ERAP1","ERAP2", "FAS", "HLA-A","HLA-B", "HLA-C", "HLA-DMA", "HLA-DMB", "HLA-DOA", "HLA-DOB","HLA-DPA1", "HLA-DPB1","HLA-DQA1","HLA-DQA2","HLA-DQB1","HLA-DRA","HLA-DRB1","HLA-DRB3","HLA-DRB4","HLA-DRB5","HLA-E","HLA-F","HLA-G","HSP90AA1","HSP90AB1","HSPA1A","HSPA1B","HSPA1L","HSPA2","HSPA4","HSPA5","HSPA6","HSPA8","HSPBP1","IFI30","IFNG","IRF1","KIR2DL1","KIR2DL2","KIR2DL3","KIR2DL4","KIR2DS1","KIR2DS2","KIR2DS4","KIR2DS5","KIR3DL1","KIR3DL2","KIR3DL3","KLRC1","KLRC2","KLRC3","KLRC4","KLRD1","LGMN","MEX3B","NFYA","NFYB","NFYC","PDIA3","PSMA7","PSMB10","PSMB11","PSMB6","PSMB8","PSMB9","PSME1","PSME2","PSME3","PSMF1","RFX5","RFXANK","RFXAP","TAP1","TAP2","TAPBP","TNF", "TTN")
df_1<-pancancer_results[[1]]$sel_cv
df_2<-pancancer_results[[2]]$sel_cv
df_1 <- df_1 %>% filter(!str_detect(gene_name, "^OR|^CS"))
df_2 <- df_2 %>% filter(!str_detect(gene_name, "^OR|^CS"))
#####trick to plot genes with pvalue equals 0
minp1<-min(df_1 %>% filter(ptrunc_cv>0) %>% pull(ptrunc_cv))
df_1<-df_1 %>% mutate(ifelse(ptrunc_cv==0,minp1-runif(1)*minp1,ptrunc_cv))

minp2<-min(df_2 %>% filter(ptrunc_cv>0) %>% pull(pmis_cv))
df_2<-df_2 %>% mutate(ifelse(ptrunc_cv==0,minp2-runif(1)*minp2,ptrunc_cv))
######End of trick######

df_1$log2FolddNdS<-log2(df_1$wnon_cv)
df_2$log2FolddNdS<-log2(df_2$wnon_cv)

df_1$Significant<-""
df_1$Significant <- ifelse(df_1$qtrunc_cv < 0.1, paste("FDR <","0.1"), "Not Sig")
df_1<-subset(df_1, !(gene_name %in% blacklist))

df_esc_trunc<-df_1 %>% filter(Significant=="FDR < 0.1")


df_2$Significant<-""
df_2$Significant <- ifelse(df_2$qtrunc_cv < 0.1, paste("FDR <","0.1"), "Not Sig")
df_2<-subset(df_2, !(gene_name %in% blacklist))
df_nonesc_trunc<-df_2 %>% filter(Significant=="FDR < 0.1")


library(openxlsx)
write.xlsx(df_esc_trunc, 'df_esc_trunc.xlsx')
write.xlsx(df_nonesc_trunc, 'df_nonesc_trunc.xlsx')

#write.xlsx(df_esc_trunc, 'df_esc_trunc.xlsx')
#write.xlsx(df_nonesc_trunc, 'df_nonesc_trunc.xlsx')


```

## Including Plots

You can also embed plots, for example:

```{r pressure, echo=FALSE}
plot(pressure)
```

Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.
