---
title: "RNAseq in *GAP1* CNV strains in YPGal"
output: html_notebook
---

Cells were grown overnight in YPGal, then backdiluted and allowed to grow for xx hours in fresh YPGal before harvesting (harvesting in mid-log)

```{r}
library(EnhancedVolcano)
library(DESeq2)
library(org.Sc.sgd.db)
library(clusterProfiler)
library(GGally)
library(tidyverse)
library(ComplexHeatmap)
require(scales)
library(dendextend)
library(patchwork)
theme_set(theme_bw(base_size = 20))

dir = "/Volumes/GoogleDrive/My Drive/Gresham Lab_Grace/france_satay/Hermes_mutagenesis_paper/analysis_new/rnaseq_r"
setwd(dir)

source("functions.R")

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

## Check replicates
```{r}
counts_df = read_tsv(paste0(data_dir, "/coverage_table_rnaseq.txt")) %>% mutate(gene = `#gene`) %>%
  filter(!is.na(gene))
# convert "genex" to systematic id
for(i in 1:nrow(counts_df)) {
  y = gff$X9[str_detect(gff$X9, paste0("ID=",counts_df$gene[i]))] %>% str_split(";")
  z = str_remove(y[[1]][str_detect(y[[1]],"locus_tag")],"locus_tag=")
  counts_df$gene[i] = z
}
counts_df = counts_df %>% select(-`#gene`)
```

```{r}
# plot umis per sample
counts_df %>%
  select(-gene) %>%
  pivot_longer(everything(), names_to = "sample", values_to = "counts") %>%
  group_by(sample) %>%
  summarise(sum(counts)) %>%
  ggplot(aes(sample, `sum(counts)`)) +
  geom_bar(stat = "identity") +
  ylab("total umis") +
  xlab("") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  scale_y_continuous(labels = comma)
```

```{r}
# convert umi counts per million (normalizes for library size)
counts_per_mill = counts_df %>% 
  pivot_longer(cols = starts_with("DGY"), names_to = "sample", values_to = "counts") %>%
  mutate(strain = str_sub(sample, 4, -3)) %>%
  left_join(strain_names) %>%
  group_by(sample) %>%
  mutate(per_million = sum(counts)/1e6) %>% 
  ungroup() %>%
  mutate(umi_million = counts/per_million)

write_csv(counts_per_mill %>% filter(sample != "DGY1740-2"), "rnaseq_counts_per_mill.csv")
```

```{r}
# check replicate correlations 
# pearson corr
# get correlations
t = counts_per_mill %>%
  mutate(strain = str_sub(sample, 4, -7),
         rep = str_sub(sample, 9, -1)) %>%
  unite(sample_rep, strain_names, rep) %>%
  dplyr::select(gene, sample_rep, umi_million) %>% 
  pivot_wider(names_from = sample_rep, values_from = umi_million) 


reps = t %>%
    dplyr::select(eu_1, eu_2, eu_3)  
cor_plots=GGally::ggpairs(reps, 
                upper = list(continuous = GGally::wrap("cor", method = "pearson"))) +
    ggtitle(LETTERS[1])
ggsave(cor_plots, file=paste0(fig_dir, "/SFigXX", LETTERS[1], ".pdf"), height = 12, width = 12)

i=2
for(strain in c("aneu", "trip1", "trip2", "trip3", "trip4", "iso", "quad")) {
  reps = t %>%
    dplyr::select(colnames(t)[str_detect(colnames(t), strain)])  
  cor_plots=GGally::ggpairs(reps, 
                upper = list(continuous = GGally::wrap("cor", method = "pearson"))) +
    ggtitle(LETTERS[i])
  ggsave(cor_plots, file=paste0(fig_dir, "/SFigXX", LETTERS[i], ".pdf"), height = 12, width = 12)
  i=i+1
}

```
