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


## Compare to other studies

Compare to differentially expressed genes in Torres et al 2007

"All expression data were therefore renormalized to 
account for the extra chromosome by averaging the log2 ratios of all genes not contained
on the disomic chromosome. This number was then subtracted from all log2 ratios in
order to make the average log2 ratios zero for genes not on the disomic chromosome.

With this correction, genes contained on all disomic chromosomes over several replicate
experiments increased in expression by an average of 1.8-fold. Replicate flasks of the
wild-type strain were also compared in order to find the extent of noise in gene
expression pattern under these conditions. A cutoff of two standard deviations from the
mean (changes of 1.8-fold) was used to remove genes for which the experimental noise
could explain the variation in gene expression. Genes were included in the further
analysis if their expression exceeded this cutoff in at least one experiment in which gene
dosage was not increased compared to wild type (3124 genes). "


```{r torresab}
torres2007 = read_csv(paste0(data_dir, "/Torres2007TableS5.csv")) %>%
  mutate(gene = str_extract(NAME, "(Y[A-Z]{2}[0-9-A-Z]+)")) %>%
  select(-NAME)

torres2007ab = torres2007 %>%
  filter(!is.na(gene)) %>%
  select(gene, contains("Fig 1AB")) %>%
   mutate(mean_wt = (`Fig 1AB: wt (A11311) batch RNA 1` + 
                       `Fig 1AB: wt (A11311) batch RNA 2` + 
                       `Fig 1AB: wt (A11311) batch RNA 3` + 
                       `Fig 1AB: wt (A11311) batch RNA 4`)/4) %>%
   mutate(across(starts_with("Fig 1AB: Dis"), ~ .x - mean_wt))

torres2007ab %>%
  left_join(gene.y.order.gasch) %>%
  pivot_longer(starts_with("Fig 1AB"), names_to = "strain", values_to = "log2FoldChange") %>%
  ggplot(aes(strain, order, fill=log2FoldChange)) +
  geom_tile() +
  scale_fill_gradient2(low="blue",mid="white", high="red", breaks = c(-3,0,3)) +
  ylab("") +
  xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y=element_blank(),
        axis.text.x = element_blank(), axis.ticks.x=element_blank())

# cor between my strains and torres AB for ESR genes
torres2007abESR = torres2007ab %>%
  mutate(mean_torresab = rowMeans(select(torres2007ab, where(is.numeric)), na.rm = TRUE)) %>%
  select(gene, mean_torresab) %>%
  left_join(diffexpress_by_strain %>%
              mutate(strain = as.character(strain))) %>%
  right_join(gene.y.order.gasch) %>%
  filter(!is.na(strain)) %>%
  filter(!is.na(mean_torresab))
```

```{r fig 5a}
# this is a figure I think####
# put some significance
# add in actual growth rates
growth_rates = read_csv(paste0(data_dir, "/rel_growth_YPGal.csv")) 

fig5a = torres2007abESR %>%
  group_by(strain) %>%
  rstatix::cor_test(mean_torresab, log2FoldChange) %>%
  left_join(strain_names) %>%
  left_join(growth_rates) %>%
  ggplot(aes(mean_r_norm, cor, color = strain_names)) +
  geom_point(size=5) + 
  scale_color_manual(values = strain_cols) +
  ylab("ESR genes correlation with\nTorres 2007 aneuploids in batch") +
  xlab("Relative growth rate\n(batch YPGal)") +
  theme(legend.position="none") 

fig5a
```

```{r sfig 7}
cors = torres2007abESR %>%
  group_by(strain) %>%
  rstatix::cor_test(mean_torresab, log2FoldChange)

sfig7 = torres2007abESR %>%
  left_join(strain_names) %>%
  ggplot(aes(mean_torresab, log2FoldChange)) +
  geom_point(alpha = 0.5) +
  facet_wrap(~strain_names) +
  xlab("log2FoldChange ESR genes Torres 2007 batch aneuploids") +
  ylab("log2FoldChange ESR genes CNV strain") +
  geom_text(data = cors %>% left_join(strain_names), 
            aes(-0.5, -3, label = paste("r = ", round(cor,2), "\n",
                                     "95% CI = [", round(conf.low,2),", ", round(conf.high,2), "]\n",
                                   "p =", p)),
            size=4) +
  geom_smooth(method ="lm", se = F)

ggsave(paste0(fig_dir,"/SFig7.png"), plot = sfig7, width = 17, height = 10, units = "in")
ggsave(paste0(fig_dir,"/SFig7.pdf"), plot = sfig7, width = 17, height = 10, units = "in")
```


Torres fig 1c

. All expression data were therefore renormalized to 
account for the extra chromosome by averaging the log2 ratios of all genes not contained
on the disomic chromosome. This number was then subtracted from all log2 ratios in
order to make the average log2 ratios zero for genes not on the disomic chromosome."

the remaining gene expression changes included a transcription pattern shared by most of the aneuploid strains and not detectable or not present in exponentially growing cultures, nor in wild-type cells or cdc28-4 and cdc23-1 mutants grown in the chemostat under phosphate-limiting conditions.
Of the 4963 genes whose expression change was greater than the control threshold (factor of 1.3) in at least one strain, 397 genes showed changed expression in 10 or more of the 14 aneuploid strains."
```{r torres c}
torres2007c = torres2007 %>%
  select(gene, contains("Fig 1C")) %>%
  filter(!is.na(gene)) %>%
  select(-contains("cdc")) %>%
  mutate(across(starts_with("Fig 1C: Dis"), ~ .x - `Fig 1C: wt (A11311) chem RNA`))

# cor between my strains and torres C for ESR genes
torres2007cESR = torres2007c %>%
  select(-contains("cdc")) %>%
  select(-"Fig 1C: wt (A11311) chem RNA") %>%
  mutate(mean_torresc = rowMeans(select(torres2007c, where(is.numeric)), na.rm = TRUE)) %>%
  select(gene, mean_torresc) %>%
  left_join(diffexpress_by_strain %>%
              mutate(strain = as.character(strain))) %>%
  right_join(gene.y.order.gasch) %>%
  filter(!is.na(strain)) %>%
  filter(!is.na(mean_torresc))
```

```{r fig 5b}
fig5b = torres2007cESR %>%
  group_by(strain) %>%
  rstatix::cor_test(mean_torresc, log2FoldChange) %>%
  left_join(strain_names) %>%
  left_join(growth_rates) %>%
  ggplot(aes(mean_r_norm, cor, color = strain_names)) +
  geom_point(size=5) + 
  scale_color_manual(values = strain_cols) +
  ylab("ESR genes correlation with\nTorres 2007 aneuploids in chemostat") +
  xlab("Relative growth rate\n(batch YPGal)") +
  theme(legend.position="none") 

fig5b
```

```{r sfig8}
cors = torres2007cESR %>%
  group_by(strain) %>%
  rstatix::cor_test(mean_torresc, log2FoldChange)

sfig8 = torres2007cESR %>%
  left_join(strain_names) %>%
  ggplot(aes(mean_torresc, log2FoldChange)) +
  geom_point(alpha = 0.5) +
  facet_wrap(~strain) +
  xlab("log2FoldChange ESR genes Torres 2007 chemostat aneuploids") +
  ylab("log2FoldChange ESR genes CNV strain") +
  geom_text(data = cors %>% left_join(strain_names), 
            aes(-0.5, -3, label = paste("r = ", round(cor,2), "\n",
                                     "95% CI = [", round(conf.low,2),", ", round(conf.high,2), "]\n",
                                   "p =", p)),
            size=4) +
  geom_smooth(method ="lm", se = F)

ggsave(paste0(fig_dir,"/SFig8.png"), plot = sfig8, width = 17, height = 10, units = "in")
ggsave(paste0(fig_dir,"/SFig8.pdf"), plot = sfig8, width = 17, height = 10, units = "in")
```


Compare to differentially expressed genes in Tsai et al 2019

this is the whole transcriptome, find "CAGE": from paper: "22 genes, termed common aneuploidy gene expression (CAGE), displaying significantly differential expression relative to haploid, were identified across all five aneuploid populations" 
```{r tsai}
tsai_all = read_csv(paste0(data_dir,"/Tsai2019NIHMS1526814-supplement-Sup_Table_1.csv"), skip=1)
tsai = tsai_all %>% filter(`padj*` < 0.05)

tsai=tsai %>%
  select(Genes, log2FoldChange) %>%
  mutate(strain = "Tsai2019") %>%
  arrange(log2FoldChange) %>%
  mutate(order = 1:nrow(tsai)) %>%
  rename(gene = "Genes") %>%
  mutate(group = rep(1:37, each=6))

tsai %>%
  bind_rows(tsai %>% select(gene, order, group) %>% 
              left_join(diffexpress_by_strain, by = "gene") %>%
              mutate(strain = as.character(strain)) %>%
              select(gene, log2FoldChange, strain, order, group)) %>%
  filter(!is.na(strain)) %>%
  group_by(group, strain) %>%
  summarise(gmean = mean(log2FoldChange)) %>% 
  ggplot(aes(group, gmean, color = strain)) +
  #geom_line(stat = "identity") +
  geom_bar(stat = "identity") +
  facet_wrap(~strain, ncol=1)
  

tsai %>%
  bind_rows(tsai %>% select(gene, order, group) %>% 
              left_join(diffexpress_by_strain, by = "gene") %>%
              mutate(strain = as.character(strain)) %>%
              select(gene, log2FoldChange, strain, order, group)) %>%
  filter(!is.na(strain)) %>%
  ggplot(aes(strain, order, fill= log2FoldChange)) +
  geom_tile() +
  scale_fill_gradient2(low="blue",mid="white", high="red", breaks = c(-3,0,3)) +
  ylab("") +
  xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y=element_blank()) #axis.text.x = element_blank(), axis.ticks.x=element_blank()
```
 
```{r fig 5c}
fig5c = tsai %>%
  select(-strain) %>%
  rename(log2FC_tsai = log2FoldChange) %>%
  left_join(diffexpress_by_strain %>%
              filter(!is.na(strain))%>%
              rename(log2FC_this = log2FoldChange)) %>%
  mutate(strain = as.character(strain)) %>%
  group_by(strain) %>%
  filter(!is.na(log2FC_this)) %>%
  rstatix::cor_test(log2FC_tsai, log2FC_this) %>%
  left_join(strain_names) %>%
  left_join(growth_rates) %>%
  ggplot(aes(mean_r_norm, cor, color = strain_names)) +
  geom_point(size=5) + 
  scale_color_manual(values = strain_cols[4:10]) +
  ylab("CAGE genes correlation\nwith Tsai 2019") +
  xlab("Relative growth rate\n(batch YPGal)") +
  theme(legend.title=element_blank())

fig5c
```

```{r sfig9}
cors = tsai %>%
  select(-strain) %>%
  rename(log2FC_Tsai = log2FoldChange) %>%
  left_join(diffexpress_by_strain) %>%
  filter(!is.na(strain)) %>%
  group_by(strain) %>%
  rstatix::cor_test(log2FC_Tsai, log2FoldChange)

sfig9 = tsai %>%
  select(-strain) %>%
  rename(log2FC_Tsai = log2FoldChange) %>%
  left_join(diffexpress_by_strain) %>%
  filter(!is.na(strain)) %>%
  left_join(strain_names) %>%
  ggplot(aes(log2FC_Tsai, log2FoldChange)) +
  geom_point(alpha = 0.5) +
  facet_wrap(~strain_names) +
  xlab("log2FoldChange CAGE genes Tsai 2019 aneuploids") +
  ylab("log2FoldChange CAGE genes CNV strain") +
  geom_text(data = cors %>% left_join(strain_names), 
            aes(-0.5, -3, label = paste("r = ", round(cor,2), "\n",
                                     "95% CI = [", round(conf.low,2),", ", round(conf.high,2), "]\n",
                                   "p =", p)),
            size=4) +
  geom_smooth(method ="lm", se = F)

ggsave(paste0(fig_dir,"/SFig9.png"), plot = sfig9, width = 17, height = 10, units = "in")
ggsave(paste0(fig_dir,"/SFig9.pdf"), plot = sfig9, width = 17, height = 10, units = "in")
```

```{r sfig10}
#check against torres AB

# 10A
cor = tsai %>%
     select(-strain) %>%
     rename(log2FC_tsai = log2FoldChange) %>%
     left_join(torres2007ab %>%
                   select(-contains("wt")) %>%
                   mutate(mean_torresab = rowMeans(select(torres2007ab, where(is.numeric)), na.rm = TRUE)) %>%
                   select(gene, mean_torresab) ) %>%
     filter(!is.na(mean_torresab)) %>%
     rstatix::cor_test(log2FC_tsai, mean_torresab) 

sfig10a = tsai %>%
     select(-strain) %>%
     rename(log2FC_tsai = log2FoldChange) %>%
     left_join(torres2007ab %>%
                   select(-contains("wt")) %>%
                   mutate(mean_torresab = rowMeans(select(torres2007ab, where(is.numeric)), na.rm = TRUE)) %>%
                   select(gene, mean_torresab) ) %>%
     filter(!is.na(mean_torresab)) %>%
  ggplot(aes(log2FC_tsai, mean_torresab)) +
  geom_point() + 
  xlab("log2FoldChange CAGE genes\nTsai 2019 aneuploids") +
  ylab("log2FoldChange CAGE genes\nTorres 2007 aneuploids in batch") +
  geom_text(data = cor, 
            aes(-0.5, -0.75, label = paste("r = ", round(cor,2), "\n",
                                     "95% CI = [", round(conf.low,2),", ", round(conf.high,2), "]\n",
                                   "p =", p)),
            size=5) +
  geom_smooth(method ="lm", se = F)

# 10B
cor = tsai %>%
     select(-strain) %>%
     rename(log2FC_tsai = log2FoldChange) %>%
     left_join(torres2007ab %>%
                   select(-contains("wt")) %>%
                   mutate(mean_torresab = rowMeans(select(torres2007ab, where(is.numeric)), na.rm = TRUE)) %>%
                   select(gene, mean_torresab) ) %>%
     filter(!is.na(mean_torresab)) %>%
  filter(!(gene %in% gasch_data$UID)) %>%
     rstatix::cor_test(log2FC_tsai, mean_torresab) 

sfig10b = tsai %>%
     select(-strain) %>%
     rename(log2FC_tsai = log2FoldChange) %>%
     left_join(torres2007ab %>%
                   select(-contains("wt")) %>%
                   mutate(mean_torresab = rowMeans(select(torres2007ab, where(is.numeric)), na.rm = TRUE)) %>%
                   select(gene, mean_torresab) ) %>%
     filter(!is.na(mean_torresab)) %>%
    filter(!(gene %in% gasch_data$UID)) %>%
  ggplot(aes(log2FC_tsai, mean_torresab)) +
  geom_point() + 
  xlab("log2FoldChange CAGE genes (excluding ESR)\nTsai 2019 aneuploids") +
  ylab("log2FoldChange CAGE genes\n(excluding ESR) Torres 2007 aneuploids\nin batch") +
  geom_text(data = cor, 
            aes(-0.5, -0.75, label = paste("r = ", round(cor,2), "\n",
                                     "95% CI = [", round(conf.low,2),", ", round(conf.high,2), "]\n",
                                   "p =", p)),
            size=5) +
  geom_smooth(method ="lm", se = F)

#check against torres C
# sfig 10C
cor = tsai %>%
  select(-strain) %>%
  rename(log2FC_tsai = log2FoldChange) %>%
  left_join(torres2007c %>%
              select(-contains("cdc")) %>%
              select(-"Fig 1C: wt (A11311) chem RNA") %>%
              mutate(mean_torresc = rowMeans(select(torres2007c, where(is.numeric)), na.rm = TRUE)) %>%
              select(gene, mean_torresc) ) %>%
  filter(!is.na(mean_torresc)) %>%
  rstatix::cor_test(log2FC_tsai, mean_torresc) 

sfig10c = tsai %>%
  select(-strain) %>%
  rename(log2FC_tsai = log2FoldChange) %>%
  left_join(torres2007c %>%
              select(-contains("cdc")) %>%
              select(-"Fig 1C: wt (A11311) chem RNA") %>%
              mutate(mean_torresc = rowMeans(select(torres2007c, where(is.numeric)), na.rm = TRUE)) %>%
              select(gene, mean_torresc) ) %>%
  filter(!is.na(mean_torresc)) %>%
  #rstatix::cor_test(log2FC_tsai, mean_torresc) %>%
  ggplot(aes(log2FC_tsai, mean_torresc)) +
  geom_point()  + 
  xlab("log2FoldChange CAGE genes\nTsai 2019 aneuploids") +
  ylab("log2FoldChange CAGE genes\nTorres 2007 aneuploids in chemostats") +
  geom_text(data = cor, 
            aes(2, -1.5, label = paste("r = ", round(cor,2), "\n",
                                     "95% CI = [", round(conf.low,2),", ", round(conf.high,2), "]\n",
                                   "p =", p)),
            size=5) +
  geom_smooth(method ="lm", se = F)

# sfig 10D
cor = tsai %>%
  select(-strain) %>%
  rename(log2FC_tsai = log2FoldChange) %>%
  left_join(torres2007c %>%
              select(-contains("cdc")) %>%
              select(-"Fig 1C: wt (A11311) chem RNA") %>%
              mutate(mean_torresc = rowMeans(select(torres2007c, where(is.numeric)), na.rm = TRUE)) %>%
              select(gene, mean_torresc) ) %>%
  filter(!is.na(mean_torresc)) %>%
  filter(!(gene %in% gasch_data$UID)) %>%
  rstatix::cor_test(log2FC_tsai, mean_torresc) 

sfig10d = tsai %>%
  select(-strain) %>%
  rename(log2FC_tsai = log2FoldChange) %>%
  left_join(torres2007c %>%
              select(-contains("cdc")) %>%
              select(-"Fig 1C: wt (A11311) chem RNA") %>%
              mutate(mean_torresc = rowMeans(select(torres2007c, where(is.numeric)), na.rm = TRUE)) %>%
              select(gene, mean_torresc) ) %>%
  filter(!is.na(mean_torresc)) %>%
  filter(!(gene %in% gasch_data$UID)) %>%
  ggplot(aes(log2FC_tsai, mean_torresc)) +
  geom_point() + 
  #scale_color_manual(values = strain_cols) +
  xlab("log2FoldChange CAGE genes (excluding ESR)\nTsai 2019 aneuploids") +
  ylab("log2FoldChange CAGE genes\n(excluding ESR) Torres 2007 aneuploids\nin chemostats") +
  geom_text(data = cor, 
            aes(1.5, -2, label = paste("r = ", round(cor,2), "\n",
                                     "95% CI = [", round(conf.low,2),", ", round(conf.high,2), "]\n",
                                   "p =", p)),
            size=5) +
  geom_smooth(method ="lm", se = F)

sfig10 = sfig10a + sfig10b + sfig10c + sfig10d + plot_annotation(tag_levels = 'A')
ggsave(paste0(fig_dir,"/SFig10.pdf"), plot = sfig10, width = 25, height = 13, units = "in")
ggsave(paste0(fig_dir,"/SFig10.png"), plot = sfig10, width = 25, height = 13, units = "in")
```



Look at Hose et al 2020
May need to cluster
```{r}
hose2020 = read_delim(paste0(data_dir, "/Hose2020elife-52063-fig2-data1-v2.txt")) %>% 
  rename(gene = "UID...1") %>%
  select(gene, contains("YPS"), contains("WA"), contains("W303"), contains("Torres 2007")) %>%
  mutate(across(-matches('gene'), ~  as.numeric(.x)))

#### CLUSTER ####
#" 861 genes with altered expression in both YPS1009_Chr12 ssd1Δ and NCYC110_Chr8 ssd1Δ mutants compared to their isogenic wild-type aneuploids" 
# so cluster based on those 2 strains
mat = as.matrix(hose2020[,c(4,7)])
rownames(mat) = hose2020$gene

# remove rows with too many NA distances --> 11 total
# from https://stackoverflow.com/questions/61469201/pheatmap-won-t-cluster-rows-na-nan-inf-in-foreign-function-call-arg-10
giveNAs = which(is.na(as.matrix(dist(mat))),arr.ind=TRUE)
tab = sort(table(c(giveNAs)),decreasing=TRUE)
checkNA = sapply(1:length(tab),function(i){
sum(is.na(as.matrix(dist(mat[-as.numeric(names(tab[1:i])),]))))
})
rmv = names(tab)[1:min(which(checkNA==0))]
mat = mat[-as.numeric(rmv),]

pbc.dist <- dist(mat)

pbc.cut.dist <- as.matrix(pbc.dist)
pbc.cut.dist[upper.tri(pbc.cut.dist, diag = T)] <- NA
pbc.cut.dist.cut <- apply(pbc.cut.dist, 1, cut.knn, nn=400)

pbc.cut.dist.cut <- 1 / pbc.cut.dist.cut
pbc.cut.dist.cut[is.na(pbc.cut.dist.cut)] <- 0
pbc.graph <- igraph::graph.adjacency(pbc.cut.dist.cut, mode = "undirected", weighted = TRUE, diag = TRUE)

# Leiden clustering Traag et al. 2019 From Louvain to Leiden: guaranteeing well-connected communities
pbc.cluster <- cluster_leiden(pbc.graph, resolution_parameter = 0.1)
pbc.gene.cluster.hose <- data.frame(Gene=rownames(mat), Cluster=0)
for (i in 1:length(pbc.cluster)) {pbc.gene.cluster.hose[pbc.gene.cluster.hose$Gene %in% pbc.cluster[[i]], "Cluster"] <- i}

# Heirarchial clustering
gene.y.hclust <- hclust(pbc.dist, 'centroid') #single # centroid
gene.y.order.hose <- as.character(gene.y.hclust$labels[gene.y.hclust$order])
gene.y.order.hose <- tibble(gene = gene.y.order.hose, order = 1:length(gene.y.order.hose))

rownames(pbc.gene.cluster.hose) = pbc.gene.cluster.hose$Gene
pbc.gene.cluster.hose= subset(pbc.gene.cluster.hose, select = c("Cluster") )
pbc.gene.cluster.hose$Cluster = as.character(pbc.gene.cluster.hose$Cluster)


hose2020 %>%
  select(gene, contains("YPS")) %>% #, contains("WA")) %>%
  replace(is.na(.), 0) %>%
  left_join(gene.y.order.hose) %>%
  pivot_longer(-c(gene, order), names_to = "strain", values_to = "log2FoldChange") %>%
  mutate(exp = "hose2020") %>%
  bind_rows(diffexpress_by_strain %>% left_join(gene.y.order.hose) %>%
              select(gene, order, strain, log2FoldChange) %>%
              mutate(strain = as.character(strain), exp = "this study")) %>%
  ggplot(aes(strain, order, fill= log2FoldChange)) +
  geom_tile() +
  scale_fill_gradient2(low="blue",mid="white", high="red", breaks = c(-3,0,3)) +
  ylab("") +
  xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y=element_blank()) + #,
  facet_wrap(~exp, scales="free_x")


# correlations
hose2020 %>%
  select(gene, `YPS1009 aneuploid wt  vs YPS1009 euploid wt  AVG n=4`) %>%
  left_join(diffexpress_by_strain %>%
              mutate(strain = as.character(strain)) %>%
              rename(log2FC_this = log2FoldChange)) %>%
  filter(!is.na(strain))%>%
  ggplot(aes(`YPS1009 aneuploid wt  vs YPS1009 euploid wt  AVG n=4`, log2FC_this)) +
  geom_point() +
  facet_wrap(~strain) +
  xlab("log2FoldChange Hose aneuploid vs euploid 2020") +
  ylab("log2FoldChange this study") +
  geom_smooth(method = "lm")

hose2020 %>%
  select(gene, `YPS1009 aneuploid ssd1  vs YPS1009 aneuploid wt  AVG n=4`) %>%
  left_join(diffexpress_by_strain %>%
              mutate(strain = as.character(strain)) %>%
              rename(log2FC_this = log2FoldChange)) %>%
  filter(!is.na(strain))%>%
  ggplot(aes(`YPS1009 aneuploid ssd1  vs YPS1009 aneuploid wt  AVG n=4`, log2FC_this)) +
  geom_point() +
  facet_wrap(~strain) +
  xlab("log2FoldChange Hose aneuploid ssd1 vs aneuploid wt 2020") +
  ylab("log2FoldChange this study") +
  geom_smooth(method = "lm")
 

hose2020 %>%
  select(gene, `Torres 2007 - Fig 1AB: Dis XI (A13771) batch RNA 1`) %>%
  left_join(diffexpress_by_strain %>%
              mutate(strain = as.character(strain)) %>%
              rename(log2FC_this = log2FoldChange)) %>%
  filter(!is.na(strain))%>%
  ggplot(aes(`Torres 2007 - Fig 1AB: Dis XI (A13771) batch RNA 1`, log2FC_this)) +
  geom_point() +
  facet_wrap(~strain) +
  xlab("log2FoldChange Torres 2007 - Fig 1AB: Dis XI (A13771) batch RNA 1") +
  ylab("log2FoldChange this study") +
  geom_smooth(method = "lm")
  
hose2020 %>%
  select(gene, `Torres 2007 - Fig 1AB: Dis XI (A13771) batch RNA 2`) %>%
  left_join(diffexpress_by_strain %>%
              mutate(strain = as.character(strain)) %>%
              rename(log2FC_this = log2FoldChange)) %>%
  filter(!is.na(strain))%>%
  ggplot(aes(`Torres 2007 - Fig 1AB: Dis XI (A13771) batch RNA 2`, log2FC_this)) +
  geom_point() +
  facet_wrap(~strain) +
  xlab("log2FoldChange Torres 2007 - Fig 1AB: Dis XI (A13771) batch RNA 2") +
  ylab("log2FoldChange this study") +
  geom_smooth(method = "lm")
```
Hose 2020 proteins that significantly changed
```{r}
hose2020sig_prot = read_delim(paste0(data_dir, "/Hose2020elife-52063-fig3-data1-v2.txt")) %>%
  rename(gene = "UID") %>%
  select(gene, contains("RNA")) %>%
  pivot_longer(-gene, names_to = "strain", values_to = "log2FoldChange") %>%
  mutate(log2FoldChange = as.numeric(log2FoldChange))
 
hose2020sig_prot %>%
  filter(strain == "YPS1009 aneuploid wt AVG vs YPS1009 euploid wt AVG RNA") %>%
  mutate(strain = "Hose wt aneuploid") %>%
  bind_rows(diffexpress_by_strain %>%
              filter(gene %in% hose2020sig$gene) %>%
              filter(padj < 0.05) %>%
              select(gene, strain, log2FoldChange) %>%
              mutate(strain = as.character(strain))) %>%
  ggplot(aes(strain, gene, fill= log2FoldChange)) +
  geom_tile() +
  scale_fill_gradient2(low="blue",mid="white", high="red", breaks = c(-3,0,3)) +
  ylab("") +
  xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y=element_blank(),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))


hose2020sig %>%
  filter(strain == "YPS1009 aneuploid wt AVG vs YPS1009 euploid wt AVG RNA") %>%
  select(-strain) %>%
  mutate(Hose = "Hose wt aneuploid") %>%
  left_join(diffexpress_by_strain %>%
              rename(log2FC_this = "log2FoldChange")) %>%
  ggplot(aes(log2FoldChange, log2FC_this)) +
  geom_point() +
  facet_wrap(~strain)
```



Look at ESR Gasch 2000 

```{r}
library(gaschYHS)
data(fig3)
gasch_data = fig3[2:nrow(fig3),colnames(fig3) != ""]
for(i in 91:ncol(gasch_data)) {
  gasch_data[,i] = as.numeric(as.character(gasch_data[,i]))
}
mat = as.matrix(gasch_data[,3:ncol(gasch_data)])

#I think I'm going to have to cluster annoyingly

pbc.dist <- dist(mat)

pbc.cut.dist <- as.matrix(pbc.dist)
pbc.cut.dist[upper.tri(pbc.cut.dist, diag = T)] <- NA
pbc.cut.dist.cut <- apply(pbc.cut.dist, 1, cut.knn, nn=100)

pbc.cut.dist.cut <- 1 / pbc.cut.dist.cut
pbc.cut.dist.cut[is.na(pbc.cut.dist.cut)] <- 0
pbc.graph <- igraph::graph.adjacency(pbc.cut.dist.cut, mode = "undirected", weighted = TRUE, diag = TRUE)

# Leiden clustering Traag et al. 2019 From Louvain to Leiden: guaranteeing well-connected communities
pbc.cluster <- cluster_leiden(pbc.graph, resolution_parameter = 0.1)
pbc.gene.cluster.gasch <- data.frame(Gene=rownames(mat), Cluster=0)
for (i in 1:length(pbc.cluster)) {pbc.gene.cluster.gasch[pbc.gene.cluster.gasch$Gene %in% pbc.cluster[[i]], "Cluster"] <- i}

# Heirarchial clustering
gene.y.hclust <- hclust(pbc.dist, 'complete')
gene.y.order.gasch <- as.character(gene.y.hclust$labels[gene.y.hclust$order])
gene.y.order.gasch <- tibble(gene = gene.y.order.gasch, order = 1:length(gene.y.order.gasch))

rownames(pbc.gene.cluster.gasch) = pbc.gene.cluster.gasch$Gene
pbc.gene.cluster.gasch = subset(pbc.gene.cluster.gasch, select = c("Cluster") )
pbc.gene.cluster.gasch$Cluster = as.character(pbc.gene.cluster.gasch$Cluster)


gasch_data %>%
  rename(gene = UID) %>%
  left_join(gene.y.order.gasch) %>%
  pivot_longer(-c(gene, NAME, order), names_to = "strain", values_to = "log2FoldChange") %>%
  ggplot(aes(strain, order, fill= log2FoldChange)) +
  geom_tile() +
  scale_fill_gradient2(low="blue",mid="white", high="red", breaks = c(-3,0,3)) +
  ylab("") +
  xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y=element_blank(),
        axis.text.x = element_blank(), axis.ticks.x=element_blank())
```

```{r}
diffexpress_by_strain %>%
  right_join(gene.y.order.gasch) %>%
  ggplot(aes(as.character(strain), order, fill= log2FoldChange)) +
  geom_tile() +
  scale_fill_gradient2(low="blue",mid="white", high="red")+#, breaks = c(-3,0,3)) +
  ylab("") +
  xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y=element_blank())
```

```{r}
t = gasch_data %>%
  rename(gene = UID) %>%
  mutate(mean_across = rowMeans(select(gasch_data, where(is.numeric)), na.rm = TRUE)) %>%
  left_join(gene.y.order.gasch) %>%
  arrange(order) %>%
  mutate(group = rep(1:124, each=7)) %>%
  group_by(group) %>%
  summarise(gmean = mean(mean_across)) %>% 
  mutate(strain = "Gasch2000\n(ESR)")

diffexpress_by_strain %>%
  right_join(gasch_data %>%
  rename(gene = UID) %>%
  mutate(mean_across = rowMeans(select(gasch_data, where(is.numeric)), na.rm = TRUE)) %>%
  left_join(gene.y.order.gasch) %>%
  arrange(order) %>%
  mutate(group = rep(1:124, each=7)) %>%
    select(gene, order, group)) %>%
  group_by(strain, group) %>%
  summarise(gmean = mean(log2FoldChange)) %>%
  filter(!is.na(strain)) %>%
  select(group, gmean, strain) %>%
  mutate(strain = as.character(strain)) %>%
  bind_rows(t) %>%
  ggplot(aes(group, gmean)) +
  #geom_line(stat = "identity") +
  geom_bar(stat = "identity") +
  facet_wrap(~strain, ncol=1)

diffexpress_by_strain %>%
  right_join(gasch_data %>%
  rename(gene = UID) %>%
  mutate(mean_across = rowMeans(select(gasch_data, where(is.numeric)), na.rm = TRUE)) %>%
  left_join(gene.y.order.gasch) %>%
  arrange(order) %>%
  mutate(group = rep(1:124, each=7)) %>%
    select(gene, order, group)) %>%
  group_by(strain, group) %>%
  summarise(gmean = mean(log2FoldChange)) %>%
  filter(!is.na(strain)) %>%
  select(group, gmean, strain) %>%
  mutate(strain = as.character(strain)) %>%
  bind_rows(t) %>%
  ggplot(aes(group, gmean)) +
  #geom_line(stat = "identity") +
  geom_bar(stat = "identity") +
  facet_wrap(~strain, ncol=1)

gasch_data %>%
  rename(gene = UID) %>%
  mutate(mean_gasch = rowMeans(select(gasch_data, where(is.numeric)), na.rm = TRUE)) %>%
  select(gene, mean_gasch) %>%
  left_join(diffexpress_by_strain %>%
              mutate(strain = as.character(strain))) %>%
  filter(!is.na(strain))%>%
  ggplot(aes(mean_gasch, log2FoldChange)) +
  geom_point() +
  facet_wrap(~strain) +
  xlab("log2FoldChange Gasch 2000 ESR") +
  ylab("log2FoldChange this study") +
  geom_smooth(method = "lm")
```

```{r}
# genes both in the ESR and significant in one or more of my strains (73 genes)
diff_express_sig %>%
  filter(gene %in% gasch_data$UID) %>%
  pivot_longer(-gene, names_to = "strain", values_to = "log2FoldChange") %>%
  mutate(exp = "This study") %>%
  bind_rows(gasch_data %>%
              rename(gene = UID) %>%
              filter(gene %in% s_genes) %>%
              select(-NAME) %>%
              pivot_longer(-c(gene), names_to = "strain", values_to = "log2FoldChange") %>%
              mutate(exp = "Gasch2000")) %>%
  bind_rows(torres2007ab %>% 
              pivot_longer(-gene, names_to = "strain", values_to = "log2FoldChange") %>%
              filter(gene %in% s_genes & gene %in% gasch_data$UID) %>%
              mutate(exp = "Torres2007AB")) %>%
  bind_rows(torres2007c %>% 
              pivot_longer(-gene, names_to = "strain", values_to = "log2FoldChange") %>%
              filter(gene %in% s_genes & gene %in% gasch_data$UID) %>%
              mutate(exp = "Torres2007Chemo")) %>%
  bind_rows(tsai_all %>%
              rename(gene = "Genes") %>%
              mutate(strain = "tsai") %>%
              select(gene, strain, log2FoldChange) %>%
              filter(gene %in% s_genes & gene %in% gasch_data$UID) %>%
              mutate(exp = "Tsai2019")) %>%
  bind_rows(hose2020 %>%
              select(gene, "YPS1009 aneuploid wt  vs YPS1009 euploid wt  AVG n=4", 
                     "WA1 aneuploid wt vs WA1 euploid wt  AVG n=3") %>%
              filter(gene %in% s_genes & gene %in% gasch_data$UID) %>%
              pivot_longer(-gene, names_to = "strain", values_to = "log2FoldChange") %>%
              mutate(log2FoldChange = as.numeric(log2FoldChange), 
                     exp = "Hose2020")) %>%
  mutate(common_name = labtools::yeast_systematic_to_common(gene)) %>%
  ggplot(aes(strain, gene, fill= log2FoldChange)) +
  geom_tile() +
  scale_fill_gradient2(low="blue",mid="white", high="red")+#, breaks = c(-3,0,3)) +
  ylab("") +
  xlab("") +
  theme_classic() +
  #theme(axis.text.y = element_blank(), axis.ticks.y=element_blank()) +
  facet_wrap(~exp, scales = "free_x", nrow = 1) 
  
```

```{r fig 5de}
# significant in one or more of my strains (436 genes)
data5de = diff_express_sig %>%
  pivot_longer(-gene, names_to = "strain", values_to = "log2FoldChange") %>%
  mutate(exp = "This study") %>%
  bind_rows(torres2007ab %>% 
              pivot_longer(-gene, names_to = "strain", values_to = "log2FoldChange") %>%
              filter(gene %in% s_genes) %>%
              mutate(exp = "Torres 2007\nbatch")) %>%
  bind_rows(torres2007c %>% 
              pivot_longer(-gene, names_to = "strain", values_to = "log2FoldChange") %>%
              filter(gene %in% s_genes) %>%
              mutate(exp = "Torres 2007\nchemo")) %>%
  bind_rows(tsai_all %>%
              rename(gene = "Genes") %>%
              mutate(strain = "tsai") %>%
              select(gene, strain, log2FoldChange) %>%
              filter(gene %in% s_genes) %>%
              mutate(exp = "Tsai 2019")) %>%
  bind_rows(hose2020 %>%
              select(gene, "YPS1009 aneuploid wt  vs YPS1009 euploid wt  AVG n=4", 
                     "WA1 aneuploid wt vs WA1 euploid wt  AVG n=3") %>%
              filter(gene %in% s_genes) %>%
              pivot_longer(-gene, names_to = "strain", values_to = "log2FoldChange") %>%
              mutate(log2FoldChange = as.numeric(log2FoldChange), 
                     exp = "Hose 2020")) %>%
  bind_rows(hose2020 %>%
              select(gene, "YPS1009 aneuploid ssd1  vs YPS1009 aneuploid wt  AVG n=4", 
                     "WA1 aneuploid ssd1 vs WA1 aneuploid wt AVG n=3") %>%
              filter(gene %in% s_genes) %>%
              pivot_longer(-gene, names_to = "strain", values_to = "log2FoldChange") %>%
              mutate(log2FoldChange = as.numeric(log2FoldChange), 
                     exp = "Hose 2020\nssd1 Δ")) %>%
  mutate(common_name = labtools::yeast_systematic_to_common(gene)) %>%
  left_join(gene.y.order)

fig5d = data5de %>%
  filter(exp == "This study") %>%
  left_join(strain_names) %>%
  mutate(exp = " \n ") %>%
  ggplot(aes(strain_names, order, fill= log2FoldChange)) +
  geom_tile() +
  scale_fill_gradient2(low="blue",mid="white", high="red", na.value="transparent")+#, breaks = c(-3,0,3)) +
  ylab("") +
  xlab("") +
  theme_void(base_size=24) +
  theme(axis.text.y = element_blank(), axis.ticks.y=element_blank(),
        axis.text.x = element_text(size = 21, margin = margin(t = 0), vjust=7)) +
  theme(legend.position="none",
        panel.border = element_rect(colour = "black", fill=NA, size=1)) 
  
fig5e = data5de %>%
  filter(exp != "This study") %>%
  ggplot(aes(strain, order, fill= log2FoldChange)) +
  geom_tile() +
  scale_fill_gradient2(low="blue",mid="white", high="red", na.value="transparent")+#, breaks = c(-3,0,3)) +
  ylab("") +
  xlab("") +
  theme_void(base_size=24) +
  facet_wrap(~factor(exp, levels = c("Torres 2007\nbatch", "Torres 2007\nchemo" ,"Tsai 2019", "Hose 2020", "Hose 2020\nssd1 Δ")),
             scales = "free_x", nrow = 1, strip.position="bottom") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1))

```

```{r fig 5}
layout <- "
ABC
DEE
DEE
"
fig5=fig5a + fig5b + fig5c + fig5d +fig5e + plot_layout(design = layout) + plot_annotation(tag_levels = 'A')

ggsave(paste0(fig_dir,"/Fig5.pdf"), plot = fig5, width = 25, height = 13, units = "in")
ggsave(paste0(fig_dir,"/Fig5.png"), plot = fig5, width = 25, height = 13, units = "in")
```



Is the ESR induced in euploid?
```{r}
counts_per_mill %>%
  right_join(gene.y.order.gasch) %>%
  mutate(strain = str_sub(sample, 4,7)) %>%
  group_by(strain, gene, order) %>%
  mutate(mean_umi_mill = mean(umi_million)) %>%
  filter(!is.na(strain)) %>%
  ggplot(aes(strain, order, fill = mean_umi_mill)) +
  geom_tile() +
  scale_fill_gradient2(low="blue",mid="white", high="red")+#, breaks = c(0,5,10)) +
  ylab("") +
  xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y=element_blank())

counts_per_mill %>%
  right_join(gene.y.order.gasch) %>%
  mutate(strain = str_sub(sample, 4,7)) %>%
  group_by(strain, gene, order) %>%
  mutate(mean_umi_mill = mean(umi_million)) %>%
  filter(!is.na(strain)) %>%
  ggplot(aes(strain, order, fill = log(mean_umi_mill))) +
  geom_tile() +
  scale_fill_gradientn(colours = terrain.colors(10)) +
#  scale_fill_gradient2(low="blue",mid="white", high="red", breaks = c(0,5,10)) +
  ylab("") +
  xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y=element_blank())
```

