ABL_SM_CRISPR_Cut_Analyses.Rmd

---
title: "ABL_SM_CRISPR_Cut_Analyses"
author: "Haider Inam"
date: '2023-03-22'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = normalizePath(".."))

library(stringr)
library(dplyr)
library(ggplot2)
library(plotly)
library(tictoc)
library(doParallel)
library(foreach)
library(RColorBrewer)


# 

```
```{r}
#Cleanup code for plotting
source("code/plotting/cleanup.R")
```

```{r,eval=F}
source("code/variantcaller/add_l298l.R")
for(i in c(1:18)){
  sample=paste("sample",i,sep = "")
  # sample="sample1"
  
  input_df_nol298l=read.csv(paste("data/Consensus_Data/novogene_lane18/",sample,"/nol298l/duplex/variant_caller_outputs/variants_unique_ann.csv",sep=""))
  input_df_l298l=read.csv(paste("data/Consensus_Data/novogene_lane18/",sample,"/l298l/duplex/variant_caller_outputs/variants_unique_ann.csv",sep=""))
  
  output_df=add_l298l(input_df_nol298l,input_df_l298l)
  write.csv(output_df,
            paste("data/Consensus_Data/novogene_lane18/",sample,"/duplex/variant_caller_outputs/variants_unique_ann.csv",sep = ""))  
}

input_df_nol298l=read.csv("data/Consensus_Data/novogene_lane18/sample9/nol298l/duplex/variant_caller_outputs/variants_unique_ann.csv")
input_df_l298l=read.csv("data/Consensus_Data/novogene_lane18/sample9/l298l/duplex/variant_caller_outputs/variants_unique_ann.csv")
  
output_df=add_l298l(input_df_nol298l,input_df_l298l)
write.csv(output_df,"data/Consensus_Data/novogene_lane18/sample9/duplex/variant_caller_outputs/variants_unique_ann.csv")  

```

```{r, eval=F}
source("code/compare_screens.R")
source("code/plotting/cleanup.R")
source("code/plotting/heatmap_plotting_function.R")
# rm(list=ls())

comparisons=read.csv("data/Consensus_Data/novogene_lane18/TwistRegion1Screen_Comparisons_Todo.csv")
comparisons=comparisons%>%filter(Completed%in%"FALSE")

for(i in 1:nrow(comparisons)){
  dirname=comparisons$dirname[i]
  pathname=paste("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/",dirname,sep = "")
  # Create directory if it doesn't already exist
  if (!file.exists(pathname)){
    dir.create(pathname)
} 
before_screen1_identifier=unlist(strsplit(comparisons$before_screen1_identifier[i],","))
after_screen1_identifier=unlist(strsplit(comparisons$after_screen1_identifier[i],","))
before_screen2_identifier=unlist(strsplit(comparisons$before_screen2_identifier[i],","))
after_screen2_identifier=unlist(strsplit(comparisons$after_screen2_identifier[i],","))
  # length(after_screen1_identifier)
# screen_compare_means=compare_screens(comparisons$before_screen1_identifier[i],
#                                      comparisons$after_screen1_identifier[i],
#                                      comparisons$before_screen2_identifier[i],
#                                      comparisons$after_screen2_identifier[i])
screen_compare_means=compare_screens(before_screen1_identifier,
                                     after_screen1_identifier,
                                     before_screen2_identifier,
                                     after_screen2_identifier)

screen_compare_means_forexport=apply(screen_compare_means,2,as.character)
# write.csv(screen_compare_means_forexport,file = paste("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/",dirname,"/screen_comparison_",dirname,".csv",sep=""))

# Plot 1. What does the heatmap look like from the average of the net growth rate?
heatmap_plotting_function(screen_compare_means,242,321,fill_variable = "netgr_obs_mean",fill_name = "Net growth rate")
# ggsave(paste("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/",dirname,"/plot1_heatmap.pdf",sep=""),width=10,height=6,units="in",useDingbats=F)
# screen_compare_means2=screen_compare_means%>%filter(alt_codon%in%twist$Codon)

# Plot 2a: What do the correlations look like for net growth rate (show mutants in text)?
ggplot(screen_compare_means,aes(x=netgr_obs_screen1,y=netgr_obs_screen2,color=resmuts,label=species))+geom_text(size=2.5)+geom_abline()+cleanup+stat_cor(method = "pearson")+labs(color="Known\nResistant\nMutant")+scale_x_continuous("Net growth rate screen 1")+scale_y_continuous("Net growth rate screen 2")
# ggsave(paste("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/",dirname,"/plot2a_Netgrowthrate_correlations_text.pdf",sep=""),width=6,height=4,units="in",useDingbats=F)

# Plot 2b: What do the correlations look like for enrichment scores (show mutants in points)?
ggplot(screen_compare_means,aes(x=netgr_obs_screen1,y=netgr_obs_screen2,label=species))+geom_point(color="black",shape=21,size=2,aes(fill=resmuts))+geom_abline()+cleanup+stat_cor(method = "pearson")+labs(fill="Known\nResistant\nMutant")+scale_x_continuous("Net growth rate screen 1")+scale_y_continuous("Net growth rate screen 2")
# ggsave(paste("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/",dirname,"/plot2b_Netgrowthrate_correlations_points.pdf",sep=""),width=6,height=4,units="in",useDingbats=F)
# Plot 2c: What do the correlations look like for enrichment scores (show mutants in text)?
ggplot(screen_compare_means,aes(x=score_screen1,y=score_screen2,color=resmuts,label=species))+geom_text(size=2.5)+geom_abline()+cleanup+ stat_cor(method = "pearson")+labs(color="Known\nResistant\nMutant")+scale_x_continuous("Enrichment score screen 1")+scale_y_continuous("Enrichment score screen 2")
# ggsave(paste("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/",dirname,"/plot2c_Enrichmentscores_correlations_text.pdf",sep=""),width=6,height=4,units="in",useDingbats=F)

# Plot 3a: Plots: what are the overall net growth rate distributions?
ggplot(screen_compare_means,aes(x=netgr_obs_mean,fill=resmuts))+geom_density(alpha=0.7)+cleanup+labs(fill="Known\nResistant\nMutant")+scale_x_continuous("Mean net growth rate of screens")
# ggsave(paste("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/",dirname,"/plot3a_Netgrowthrate_distributions_resmuts.pdf",sep=""),width=6,height=4,units="in",useDingbats=F)
# Plot 3b: Plots: what are the net growth rate distributions?

library(reshape2)
screen_compare_melt=melt(screen_compare_means%>%dplyr::select(species,netgr_obs_screen1,netgr_obs_screen2),id.vars = "species",measure.vars =c("netgr_obs_screen1","netgr_obs_screen2"),variable.name = "Condition",value.name = "netgr_obs")


ggplot(screen_compare_melt,aes(x=netgr_obs,fill=Condition))+
  geom_density(alpha=0.7)+
  cleanup+
  scale_x_continuous("Net growth rate observed")+
  scale_fill_discrete(labels=c("Screen 1","Screen 2"))
# ggsave(paste("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/",dirname,"/plot3b_Netgrowthrate_distributions.pdf",sep=""),width=6,height=4,units="in",useDingbats=F)

}

# dirname="K562_Medium_rep1vs2"

```


Doing some quick analysis on buried hydrophobic vs exposed residues

Quick analysis to do on the secondary structure:
The buried hydrophobic core is probably going to be more susceptible to polar/hydrophilic substitutions, and not as much to hydrophobic substitutions... do you see that in the il3 independence data? Answer: yes, I see a slight singal but it's not huge
```{r}
dssp=read.csv("data/DSSP_SolventAccessibility_ABL/2hyy_dspp.csv",header = T)
dssp=dssp%>%mutate(RESIDUE=as.numeric(RESIDUE),
                   ACC=as.numeric(ACC),AA=gsub("<ca>","",AA),
                   SS=case_when(STRUCTURE%in%c("E","B")~"b-sheet",
                                STRUCTURE%in%c("H","G","I")~"a-helix",
                                STRUCTURE%in%"T"~"turn",
                                T~"undefined"))%>%dplyr::select(-"STRUCTURE")

dssp=dssp%>%rename(protein_start=RESIDUE)

il3=read.csv("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/cross-replicate/baf3_IL3_rep1vsrep2_ft/screen_comparison_baf3_IL3_low_rep1vsrep2_ft.csv")
il3=il3%>%dplyr::select(species,protein_start,netgr_obs_mean)
il3=il3%>%group_by(protein_start)%>%summarize(netgr_obs_mean=mean(netgr_obs_mean))

il3_dssp=merge(il3,dssp,by="protein_start")
il3_dssp=il3_dssp%>%mutate(exposed=case_when(ACC>=40~"Exposed",
                                            T~"Buried"))
# il3_dssp=il3_dssp%>%filter(exposed%in%"Exposed")
ggplot(il3_dssp,aes(x=ACC))+geom_histogram()



ggplot(il3_dssp,aes(x=protein_start,y=netgr_obs_mean,color=exposed))+
  geom_point()+
  facet_wrap(~exposed,nrow=2)+
  scale_y_continuous("Mean net growth rate at residue")+
  scale_x_continuous("Residue on the ABL Kinase")+
  labs(color="Solvent \nAccessibility")+
  cleanup+theme(legend.position = c(.9,.85))+
  theme(
  strip.background = element_blank(),
  strip.text.x = element_blank()
)+theme(legend.background = element_rect(
                                  size=0.5, linetype="solid", 
                                  colour ="black"))
# ggsave("data/DSSP_SolventAccessibility_ABL/il3solvent_accessibility_v1.pdf",width=6,height = 4,units="in",useDingbats=F)


ggplot(il3_dssp,aes(x=protein_start,y=netgr_obs_mean,color=exposed))+
  geom_point()+
  # facet_wrap(~exposed,nrow=2)+
  scale_y_continuous("Mean net growth rate at residue")+
  scale_x_continuous("Residue on the ABL Kinase")+
  labs(color="Solvent \nAccessibility")+cleanup
# ggsave("data/DSSP_SolventAccessibility_ABL/il3solvent_accessibility_v2.pdf",width=6,height = 4,units="in",useDingbats=F)


ggplot(il3_dssp,aes(x=ACC,y=netgr_obs_mean,color=exposed))+
  geom_point()+
  scale_y_continuous("Mean net growth rate at residue")+
  scale_x_continuous("DSSP Solvent Accessibility at Residue")+
  cleanup+
  labs(color="Solvent \nAccessibility")+
  theme(legend.position = c(.9,.25))+
  theme(legend.background = element_rect(
                                  size=0.5, linetype="solid", 
                                  colour ="black"))
# ggsave("data/DSSP_SolventAccessibility_ABL/il3solvent_accessibility_v3.pdf",width=6,height = 4,units="in",useDingbats=F)



il3=read.csv("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/cross-replicate/baf3_IL3_rep1vsrep2_ft/screen_comparison_baf3_IL3_low_rep1vsrep2_ft.csv")
il3=il3%>%dplyr::select(species,ref_aa,protein_start,alt_aa,netgr_obs_mean)


il3_dssp=merge(il3,dssp,by.x=c("protein_start","ref_aa"),by.y=c("protein_start","AA"))
il3_dssp=il3_dssp%>%mutate(exposed=case_when(ACC>=40~"Exposed",
                                            T~"Buried"))

# il3_dssp=il3_dssp%>%mutate(hydrophobic=case_when(alt_aa%in%c("A","V","I","M","L","F","Y","W")~"hydrophobic",
#                                                  T~"other"))
ggplot(il3_dssp,aes(x=alt_aa,y=netgr_obs_mean))+geom_boxplot()+facet_wrap(~exposed)
ggplot(il3_dssp,aes(x=exposed,y=netgr_obs_mean))+geom_boxplot()+facet_wrap(~exposed)


ggplot(il3_dssp,aes(x=factor(protein_start),y=netgr_obs_mean,fill=exposed))+
  geom_boxplot()+
  # facet_wrap(~exposed,nrow=2)+
  scale_y_continuous("Mean net growth rate at residue")+
  scale_x_discrete("Residue on the ABL Kinase")+
  labs(fill="Solvent \nAccessibility")+cleanup+
  theme(axis.text.x=element_text(angle=90, hjust=1))

```



Doing some quick analysis on buried hydrophobic vs exposed residues
```{r,eval=F,include=F}
dssp=read.csv("data/DSSP_SolventAccessibility_ABL/2gqg_dspp.csv",header = T)
dssp=dssp%>%mutate(RESIDUE=as.numeric(RESIDUE),ACC=as.numeric(ACC),AA=gsub("<ca>","",AA))
dssp=dssp%>%rename(protein_start=RESIDUE)

il3=read.csv("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/cross-replicate/baf3_IL3_rep1vsrep2_ft/screen_comparison_baf3_IL3_low_rep1vsrep2_ft.csv")
il3=il3%>%select(species,protein_start,netgr_obs_mean)
il3=il3%>%group_by(protein_start)%>%summarize(netgr_obs_mean=mean(netgr_obs_mean))

il3_dssp=merge(il3,dssp,by="protein_start")
il3_dssp=il3_dssp%>%mutate(exposed=case_when(ACC>=40~"Exposed",
                                            T~"Buried"))
ggplot(il3_dssp,aes(x=ACC))+geom_histogram()



ggplot(il3_dssp,aes(x=protein_start,y=netgr_obs_mean,color=exposed))+
  geom_point()+
  facet_wrap(~exposed,nrow=2)+
  scale_y_continuous("Mean net growth rate at residue")+
  scale_x_continuous("Residue on the ABL Kinase")+
  labs(color="Solvent \nAccessibility")+
  cleanup+theme(legend.position = c(.9,.85))+
  theme(
  strip.background = element_blank(),
  strip.text.x = element_blank()
)+theme(legend.background = element_rect(
                                  size=0.5, linetype="solid", 
                                  colour ="black"))
# ggsave("data/DSSP_SolventAccessibility_ABL/il3solvent_accessibility_v1.pdf",width=6,height = 4,units="in",useDingbats=F)


ggplot(il3_dssp,aes(x=protein_start,y=netgr_obs_mean,color=exposed))+
  geom_point()+
  # facet_wrap(~exposed,nrow=2)+
  scale_y_continuous("Mean net growth rate at residue")+
  scale_x_continuous("Residue on the ABL Kinase")+
  labs(color="Solvent \nAccessibility")+cleanup
# ggsave("data/DSSP_SolventAccessibility_ABL/il3solvent_accessibility_v2.pdf",width=6,height = 4,units="in",useDingbats=F)


ggplot(il3_dssp,aes(x=ACC,y=netgr_obs_mean,color=exposed))+
  geom_point()+
  scale_y_continuous("Mean net growth rate at residue")+
  scale_x_continuous("DSSP Solvent Accessibility at Residue")+
  cleanup+
  labs(color="Solvent \nAccessibility")+
  theme(legend.position = c(.9,.25))+
  theme(legend.background = element_rect(
                                  size=0.5, linetype="solid", 
                                  colour ="black"))
# ggsave("data/DSSP_SolventAccessibility_ABL/il3solvent_accessibility_v3.pdf",width=6,height = 4,units="in",useDingbats=F)
```
Doing some quick DDG analysis
```{r}
ddg=read.table("data/DDG_ABL/2hyy_ddg.tsv")
colnames(ddg)=c("species","ddg","ddg_sd")
ddg=ddg%>%mutate(protein_start=234+as.numeric(gsub("[^0-9]","",species)),
                 ref_aa=substr(species,1,1),
                 alt_aa=sub(".+(.)$", "\\1", species),
                 species=paste(ref_aa,protein_start,alt_aa,sep = ""))
# class(ddg$protein_start)
ddg=ddg%>%filter(protein_start%in%c(242:322))

il3=read.csv("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/cross-replicate/baf3_IL3_rep1vsrep2_ft/screen_comparison_baf3_IL3_low_rep1vsrep2_ft.csv")
il3=il3%>%select(species,protein_start,netgr_obs_mean)
ddgil3=merge(ddg,il3,by=c("species","protein_start"))
ggplot(ddgil3,aes(x=ddg,y=netgr_obs_mean))+geom_point()
ggplot(ddgil3,aes(x=ddg,y=netgr_obs_mean,label=species))+geom_text()
```

```{r}
source("code/resmuts_adder.R",include=F,eval=F)
ic50data_all_sum=read.csv("output/ic50data_all_confidence_intervals_raw_data.csv",row.names = 1)
imat=read.csv("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/cross-replicate/baf3_Imat_low_rep1vsrep2_ft/screen_comparison_baf3_Imat_low_rep1vsrep2_ft.csv")
# imat=read.csv("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/cross-replicate/k562_Imat_medium_rep1vsrep2/screen_comparison_k562_imat_medium_rep1vsrep2.csv")
imat=imat%>%select(species,ref_aa,protein_start,alt_aa,netgr_obs_mean)
# x=imat%>%filter(species%in%"V299L")
imat=resmuts_adder(imat)
resmuts_merged=merge(imat%>%filter(resmuts%in%"TRUE"),ic50data_all_sum%>%dplyr::select(species,netgr_pred_model,netgr_pred_mean),by="species")


ggplot(resmuts_merged,aes(x=netgr_pred_mean,y=netgr_obs_mean,label=species))+geom_text()+theme_bw()+geom_abline()
```


```{r}
data=read.csv("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/cross-replicate/baf3_Imat_high_rep1vsrep2_ft/screen_comparison_baf3_Imat_high_rep1vsrep2_ft.csv",header = T,stringsAsFactors = F)
data=read.csv("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/cross-replicate/k562_imat_medium_rep1vsrep2/screen_comparison_k562_imat_medium_rep1vsrep2.csv",header = T,stringsAsFactors = F)
source("code/res_residues_adder.R")
source("code/cosmic_data_adder.R")
# data=read.csv("output/ABLEnrichmentScreens/Imat_Enrichment_bgmerged_2.22.23.csv",header = T,stringsAsFactors = F)
# class(data)
data=res_residues_adder(data)
data=cosmic_data_adder(data)
# data=data%>%rowwise()%>%mutate(netgr_mean=mean(netgr_obs.x,netgr_obs.y))
data=data%>%rowwise()%>%mutate(netgr_mean=mean(netgr_obs_screen1,netgr_obs_screen2))
data$resmut_cosmic="neither"
data[data$cosmic_present%in%TRUE,"resmut_cosmic"]="cosmic"
data[data$resmuts%in%TRUE,"resmut_cosmic"]="resmut"
data$resmut_cosmic=factor(data$resmut_cosmic,levels=c("neither","cosmic","resmut"))


ggplot(data%>%filter(protein_start>=242,protein_start<=494,n_nuc_min%in%1),aes(x=reorder(species,-netgr_obs_mean),y=netgr_obs_mean,fill=resmut_cosmic))+geom_col()+theme(axis.text.x=element_text(angle=90, hjust=1))+scale_y_continuous(name="Enrichcment Score")+scale_x_discrete(name="Mutant")+scale_fill_manual(name="Resistance Status",labels=c("Never Seen","Sanger Mutant","Known Resistant Mutant"),values=c("gray","orange","red"))+cleanup+theme(legend.position = "none")


ggplot(data%>%group_by(species)%>%mutate(n_species=n())%>%filter(protein_start>=242,protein_start<=494,n_species%in%1),aes(x=reorder(species,-netgr_mean),y=netgr_mean,fill=resmut_cosmic))+geom_col()+theme(axis.text.x=element_text(angle=90, hjust=1))+scale_y_continuous(name=bquote('Net Growth Rate '(Hours^-1)))+scale_x_discrete(name="Mutant")+scale_fill_manual(name="Resistance Status",labels=c("Never Seen","Sanger Mutant","Known Resistant Mutant"),values=c("gray","orange","red"))+cleanup+theme(legend.position = "none")
ggsave("ImatinibEnrichment_Netgr_Distribution.pdf",width=8,height=4,units="in",useDingbats=F)

ggplot(data%>%filter(protein_start>=242,protein_start<=494,n_nuc_min%in%1,score_mean>=0),aes(x=reorder(species,-netgr_mean),y=netgr_mean,fill=resmut_cosmic))+geom_col()+
  scale_y_continuous(name=bquote('Net Growth Rate '(Hours^-1)))+
  scale_x_discrete(name="Mutant")+
  scale_fill_manual(name="Resistance Status",labels=c("Never Seen","Sanger Mutant","Known Resistant Mutant"),values=c("gray","orange","red"))+
  theme(axis.text.x=element_text(angle=90, hjust=1),
        axis.text=element_text(size=6),
        panel.grid.major = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.background = element_blank())+
  theme(legend.position = "none")
ggsave("ImatinibEnrichment_Netgr_Distribution_zoom.pdf",width=6,height=4,units="in",useDingbats=F)

ggplot(data%>%filter(!ct_screen1_after%in%.5,!ct_screen2_after%in%.5,protein_start>=242,protein_start<=494),aes(x=netgr_obs_screen1,y=netgr_obs_screen2))+
  geom_point(color="black",shape=21,aes(fill=resmut_cosmic))+
  geom_abline()+
  scale_x_continuous(name="Net Growth Rate D0 D2")+
  scale_y_continuous(name="Net Growth Rate D0 D4")+
  labs(fill="Sanger\n Mutant")+
  scale_fill_manual(values=c("gray90","orange","red"))+
  cleanup+theme(legend.position = "none")
ggsave("ImatinibEnrichment_Plot_netgr_4.13.23.pdf",width=4,height=4,units="in",useDingbats=F)
x=data%>%filter(!ct_screen1_after%in%.5,!ct_screen2_after%in%.5,protein_start>=242,protein_start<=322)
cor(x$netgr_obs_screen1,x$netgr_obs_screen2)


ggplot(data%>%filter(protein_start>=242,protein_start<=494,n_nuc_min%in%1),aes(x=resmut_cosmic,y=netgr_mean,fill=resmut_cosmic))+
  geom_violin(color="black")+
  # geom_boxplot(color="black")+
  geom_jitter(color="black", size=1,width=.1, alpha=0.9)+
  scale_fill_manual(values=c("gray90","orange","red"))+scale_y_continuous(name=bquote('Net Growth Rate '(Hours^-1)))+scale_x_discrete("Resistance Status",labels=c("Never\n Seen","Sanger","Known\nResistant"))+cleanup+theme(legend.position = "none")
# ggsave("ImatinibEnrichment_netgr_BoxPlot_2.21.23.pdf",width=4,height=4,units="in",useDingbats=F)

plotly=ggplot(data%>%filter(protein_start>=242,protein_start<=494,n_nuc_min%in%1),aes(x=reorder(species,-netgr_mean),y=netgr_mean,fill=resmut_cosmic))+geom_col()+theme(axis.text.x=element_text(angle=90, hjust=1))+scale_y_continuous(name="Enrichcment Score")+scale_x_discrete(name="Mutant")+guides(fill = guide_legend(title = "Clinically\n Observed \n Resmut"))+scale_fill_manual(values=c("gray","orange","red"))
ggplotly(plotly)

ggplot(data%>%filter(!ct_screen1_after%in%.5,!ct_screen2_after%in%.5,protein_start>=242,protein_start<=494),aes(x=netgr_obs_screen1,y=netgr_obs_screen2))+
  geom_point(color="black",shape=21,aes(fill=resmut_cosmic))+
  geom_abline()+
  scale_x_continuous(name="Net Growth Rate D0 D2")+
  scale_y_continuous(name="Net Growth Rate D0 D4")+
  labs(fill="Sanger\n Mutant")+
  scale_fill_manual(values=c("gray90","orange","red"))+
  theme(legend.position = "none")

```
Plotting the short et al data for Twinstrand
```{r}
data_old=read.csv("output/ABLEnrichmentScreens/Imat_Enrichment_bgmerged_2.22.23.csv",header = T,stringsAsFactors = F)
data_old=data_old%>%rowwise()%>%mutate(netgr_obs_mean=mean(netgr_obs_screen1,netgr_obs_screen2))
data_old=data_old%>%filter(protein_start%in%c(323:494))
data_new=read.csv("output/ABLEnrichmentScreens/ABL_Region1_Lane18_Comparisons/cross-replicate/baf3_Imat_low_rep1vsrep2_ft/screen_comparison_baf3_Imat_low_rep1vsrep2_ft.csv",header = T,stringsAsFactors = F)
data_new=data_new%>%filter(protein_start%in%c(242:322))
data_merged=rbind(data_new,data_old)
short_mutants=read.csv("data/Short_et_al_fig1/short_et_al_3.12.23.csv",header = T)
short_mutants$preexisting_all=T
short_mutants=short_mutants%>%dplyr::select(-Index)%>%rename(species=Species)
data_merged=merge(data_merged,short_mutants,by="species",all.x = T)
data_merged[!data_merged$preexisting_all%in%T,"preexisting_all"]=F
nrow(data_merged%>%filter(preexisting_all%in%T))
source("code/resmuts_adder.R")
data_merged=resmuts_adder(data_merged)
data_merged$resmut_cosmic="neither"
data_merged[data_merged$preexisting_all%in%TRUE,"resmut_cosmic"]="preexisting"
data_merged[data_merged$resmuts%in%TRUE,"resmut_cosmic"]="resmut"
data_merged$resmut_cosmic=factor(data_merged$resmut_cosmic,levels=c("neither","preexisting","resmut"))


ggplot(data_merged%>%filter(protein_start>=242,protein_start<=322,n_nuc_min%in%1),aes(x=resmut_cosmic,y=netgr_obs_mean,fill=resmut_cosmic))+
  geom_violin(color="black")+
  # geom_boxplot(color="black")+
  geom_jitter(color="black", size=1,width=.1, alpha=0.9)+
  scale_fill_manual(values=c("gray90","#9F416D","#F39826"))+scale_y_continuous(name=bquote('Net Growth Rate '(Hours^-1)))+scale_x_discrete("Resistance Status",labels=c("Undetected","Pre-existing ALL","Known\nResistant"))+cleanup+theme(legend.position = "none")
ggsave("data/Short_et_al_fig1/shortetal_netgr_BoxPlot_4.14.23.pdf",width=4,height=4,units="in",useDingbats=F)

# x=data_merged%>%filter(protein_start>=242,protein_start<=322,n_nuc_min%in%1,preexisting_all%in%T)
```