---
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"))
```

## Differential analysis with DESeq

```{r}
count_df_cds = insert_profiles %>%
  filter(type == "cds") %>%
  select(gene, n_insertions, sample) %>%
  pivot_wider(names_from = "sample", values_from = "n_insertions")
```

```{r}
# first all cds
diff_analysis_cds = NULL
gsea_cds=NULL

for(strain in c("1728", "1734", "1736", "1740", "1744", "1747", "1751")) {
  counts = as.data.frame(count_df_cds) %>% dplyr::select(`1657_1`, `1657_2`, strain)
  rownames(counts) = count_df_cds$gene
  coldata = data.frame(type = c("wt", "wt", "cnv"), sample = c('1657_1','1657_2',strain), row.names = colnames(counts))
  coldata$type = factor(coldata$type, levels = c('wt', 'cnv'))
  dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = coldata,
                              design = ~ type)
  dds <- DESeq(dds)
  res <- results(dds, alpha=0.05)
  diff_analysis_cds = bind_rows(diff_analysis_cds, 
                                    as_tibble(res) %>% 
                                      mutate(gene = rownames(res), strain = strain))
  geneList = res$log2FoldChange#[-1] #remove GRESHAMGFP
  gene.df <- bitr(rownames(res), fromType = "ORF", #[-1]
        toType = "ENTREZID",
        OrgDb = org.Sc.sgd.db::org.Sc.sgd.db)
  names(geneList) = gene.df$ENTREZID
  geneList = sort(geneList, decreasing = TRUE)
  ego <- gseGO(geneList     = geneList,
             OrgDb        = org.Sc.sgd.db::org.Sc.sgd.db,
             keyType = "ENTREZID",
              ont          = "BP",
              minGSSize    = 10,
              maxGSSize    = 500,
              pvalueCutoff = 0.05,
              pAdjustMethod="fdr",
              verbose      = FALSE,
              by="fgsea",
             seed = 1)
  ego = clusterProfiler::simplify(ego)
  gsea_cds = bind_rows(gsea_cds, as_tibble(ego@result) %>% mutate(sample = strain))
}
#write out so this doesn't have to run again
# these will probably also be supplementary files
write_csv(gsea_cds, paste0(data_dir,"/gsea_cds.csv"))
write_csv(diff_analysis_cds, paste0(data_dir,"/diff_analysis_cds.csv"))
```


```{r fig 3c}
diff_analysis_cds = read_csv(paste0(data_dir,"/diff_analysis_cds.csv"))

f3c = diff_analysis_cds %>%
  filter(padj < 0.05) %>%
  mutate(strain = as.character(strain)) %>%
  left_join(strain_names) %>%
  mutate(common_name = labtools::yeast_systematic_to_common(gene)) %>%
  filter(common_name != "GRESHAMGFP") %>%
  left_join(cnv_genes, by = c("strain" = "sample", "gene" = "gene", "strain_names" = "strain_names")) %>%
  mutate(copy_number = if_else(is.na(copy_number),"",as.character(copy_number))) %>%
  ggplot(aes(strain_names, common_name, fill = log2FoldChange)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
  theme_grey() +
  theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_line(size = 0.5, linetype = "solid",
                                   colour = "black")) + 
  xlab("") +
  ylab("") + 
  geom_text(aes(strain_names, common_name, label=copy_number)) +
  labs(fill = "log2FoldChange\nfor p.adj < 0.05") 

# check size ####
ggsave(paste0(fig_dir,"/Fig3C.pdf"), plot = f3c, width = 9, height = 5, units = "in")