---
title: "Transposon mutagenesis in *GAP1* CNV strains"
author: "Grace Avecilla"
output: html_notebook
---


```{r}
library(EnhancedVolcano)
library(DESeq2)
library(org.Sc.sgd.db)
library(clusterProfiler)
library(tidyverse)
library(patchwork)
library(ggbeeswarm)
library(ggpubr)
theme_set(theme_bw(base_size = 20))

data_dir = "/Volumes/GoogleDrive/My Drive/Gresham Lab_Grace/france_satay/Hermes_mutagenesis_paper/data"

fig_dir = "/Volumes/GoogleDrive/My Drive/Gresham Lab_Grace/france_satay/Hermes_mutagenesis_paper/figures"

source("/Volumes/GoogleDrive/My Drive/Gresham Lab_Grace/france_satay/Hermes_mutagenesis_paper/analysis_w_PS_corrections/functions.R")

seq_runs = c("bgi1","bgi2", "nyc1", "nyc2")
```

```{r}
# get all genes excluding dubious orfs and their copy number
strain_cns = read_csv(paste0(data_dir, "/gene_median_relative_depth_DNA_corrected_v3/gene_median_relative_depth_DNA-Table 1.csv")) %>% 
  select(Gene, contains("cor")) %>%
  pivot_longer(-Gene, names_to = "strain", values_to = "copy_number") %>%
  mutate(strain = str_extract(strain,"[0-9]{4}")) %>%
  left_join(strain_names)

all_genes = strain_cns$Gene
```

# Correlation between different sequencing methods
Using insertions per gene
```{r}
get_ipg = function(x, ver){
  data = read_tsv(paste0(data_dir, "/hpc_output/",ver, "/", x), 
                  col_types = cols(CDS=col_character(),`#insertion` = col_double())) %>%
    mutate(sample = str_sub(x, 1, -22), gene=NA,
           run=ver)
}
ipg_forcor = NULL
for(i in seq_runs) {
  ipg_files = list.files(path = paste0(data_dir, "/hpc_output/", i, "/"), pattern = '*insertionPerGene.txt')
  ipg_t = map(ipg_files, get_ipg, ver=i)
  ipg_df_t = do.call(rbind, ipg_t)
  ipg_forcor = ipg_forcor %>% 
    bind_rows(ipg_df_t %>% unite(sample_rep, sample, run))
}
```

 
# Get insertion profiles
```{r get insertion profiles, message=FALSE}
#get insertion profiles for all genes
insert_profiles = map(all_genes_systematic, get_insert_profile, data = rpp_df)
insert_profiles = do.call(rbind, insert_profiles)

# add in zero categories
t=insert_profiles %>% right_join(tibble(all_genes_systematic), by= c("gene"="all_genes_systematic")) %>%
  mutate_at(vars(sample, gene, type), factor) 
t=t %>% 
  tidyr::expand(sample, gene, type)
insert_profiles = insert_profiles %>% #filter(!is.na(sample)) %>% 
  right_join(t) %>%
  mutate(n_insertions = replace_na(n_insertions, 0), normalized_insertions = replace_na(normalized_insertions, 0)) %>%
  distinct()

# write out so this doesn't need to be done again
write_csv(insert_profiles, paste0(data_dir,"/insertion_profiles.csv"))
```


```{r get insertion profile csv}
insert_profiles = read_csv(paste0(data_dir,"/insertion_profiles.csv"))
```

```{r}
inserts_per_mill_cds = insert_profiles %>% 
  left_join(strain_names, c("sample" = "strain")) %>%
  left_join(summary_rpp) %>%
  filter(type %in% c("cds")) %>%
  group_by(sample, gene) %>%
  mutate(inserts = sum(n_insertions)) %>% ungroup() %>%
  select(sample, strain_names, gene, inserts, total_sites) %>%
    distinct() %>%
  mutate(inserts_per_mill = inserts/(total_sites/1e6))
```


```{r fig 3a}
# check all that have inserts
f3a=inserts_per_mill_cds %>% 
  mutate(median_inserts_per_mill = median(inserts_per_mill)) %>%
  filter(gene %in% pull(x,gene)) %>%
  group_by(gene) %>%
  filter(!(sample %in% c("1657_1", "1657_2"))) %>%
  mutate(all_strains = sum(inserts)) %>%
  filter(all_strains > 0) %>%
  mutate(all_g0 = all(inserts > 0)) %>%
  filter(all_g0 == T) %>%
  ungroup() %>%
  left_join(cnv_genes) %>%
  mutate(common_name = labtools::yeast_systematic_to_common(gene)) %>%
  filter(!is.na(common_name)) %>%
  ggplot(aes(common_name, inserts_per_mill, fill = strain_names)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.8) +
  scale_fill_manual(values = strain_cols[4:10]) +
  geom_abline(aes(intercept = median_inserts_per_mill, slope=0),
                color="azure4", size=0.75) +
  geom_text(aes(label = copy_number, x = common_name, y = inserts_per_mill), 
            position = position_dodge(width = 0.8), vjust = -0.6, size=5.5) +
  ylab("Insertions per million") +
  theme(legend.title=element_blank()) +
  xlab("")
  

ggsave(paste0(fig_dir,"/Fig3a.pdf"), plot = f3a, width = 12, height = 6.5, units = "in")

```
