---
title: "Miscellaneous analysis for pangenome stuff"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  cache = TRUE, cache.path = '.cache_knitr/misc/',
  fig.width = 7, fig.height = 4, dpi = 350,
  fig.path='figure/0-misc/',
  dev = c('png', 'pdf')
)
```

## Setup
```{r, cache=FALSE}
library(GenomicRanges)
library(scales)
library(dplyr)
library(tidyr)
library(ggplot2)
library(tibble)
library(readr)
library(glue)
library(ggridges)
library(ggupset)

theme_set(theme_classic())
theme_update(
  panel.background = element_rect(fill = "transparent", colour = NA),
  plot.background = element_rect(fill = "transparent", colour = NA),
  legend.background = element_rect(fill = "transparent", colour = NA)
)

source("modules/frequently_used_items.R")
source("modules/get_chrom_sizes.R")
```

## Chromosome sizes
Prepare data frame of the chrom sizes. We also change the order of the species to follow ecological clades in the presentation.
```{r}
chrom_sizes_df <- bind_rows(chrom_sizes_list, .id = "species")
chrom_sizes_df$species <- forcats::fct_relevel(chrom_sizes_df$species,
  'mayZeb', 'troMau', 'astCal', 'otoArg', 'aulStu', 'copChr', 'rhaChi', 'rhaChi2')

head(chrom_sizes_df)
```

Ridgeplot of the contig size histogram.
```{r contigdist_ridges, fig.width = 4.7, fig.height = 6}
ggplot(chrom_sizes_df, aes(x = length, y = species, fill = species)) +
  geom_density_ridges(scale = 0.92, stat='binline') +
  theme_ridges() + 
  scale_x_log10("contig/scaffold lengths (bp)",
        breaks = trans_breaks("log10", function(x) 10^x),
        labels = trans_format("log10", math_format(10^.x))) +
  coord_cartesian(clip = "off") +
  theme(legend.position = "none") +
  ylab(NULL) +
  scale_y_discrete(limits=rev, breaks = color_schema_species$species_abrv) +
  scale_fill_manual(values = color_schema_species$color, breaks = color_schema_species$species_abrv)
```

Faceetted histogram version
```{r contigdist_facet, fig.width = 6, fig.height = 8}
ggplot(chrom_sizes_df, aes(x = length, fill = species)) +
  geom_histogram(bins=60) +
  facet_wrap(.~species, scales = 'free_y', ncol=1) +
  scale_x_log10("contig/scaffold lengths (bp)",
        breaks = trans_breaks("log10", function(x) 10^x),
        labels = trans_format("log10", math_format(10^.x))) +
  ylab('frequency') +
  theme(legend.position = "none") +
  scale_fill_manual(values = color_schema_species$color, breaks = color_schema_species$species_abrv)
```

## Some additional pairwise comparisons stuff
Most plots are in the Snakemake pipeline, but here are some one-off things.
```{r unnormalised_pairwise_stats, fig.width=10, fig.height=3.5}
df_bubble_stats_L50 <- readr::read_delim(
  paste0(Sys.getenv("HOME"),
    "/code/malawi_transposon/cloud/workflow/malawi_haplochromines_conserv/results/",
    "bubble_agg_stats-L50.txt") ) %>% 
  mutate(
    species_pair = paste(species1, species2, sep='-'),
    species1 = factor(species1, levels = color_schema_species$species_abrv),
    species2 = factor(species2, levels = color_schema_species$species_abrv)
  ) %>% 
  pivot_longer(c(n_bubble, mean_bubble_size, total_bubble_size))

df_bubble_stats_L50 %>% 
  ggplot() +
  geom_jitter(aes(x=species1, y=value, color=species2, shape=species2), width=0.12) +
  facet_wrap(.~name, nrow=1, scales = 'free_y') +
  xlab("backbone") + ylab(NULL) +
  labs(color = "aligned query", shape = "aligned query") +
  scale_color_manual(values = color_schema_species$color, breaks = color_schema_species$species_abrv) +
  scale_shape_manual(values = c(19,19,19,15,15,15,17,17), breaks = color_schema_species$species_abrv) +
  theme(legend.position="top", axis.text.x = element_text(angle = 45, hjust = 1)) +
  guides(colour = guide_legend(nrow = 1))
```


## Missing BUSCO genes
```{r}
busco_dir <- paste0(Sys.getenv("HOME"), 
  "/code/malawi_transposon/cloud/genome/longread-decompress/renamed/BUSCO_supp_tables/")

busco_lookup <- read_delim(paste0(busco_dir, "links_to_ODB10.txt"),
  comment = '#', delim = '\t', col_names = c('busco_id', 'gene_description', 'orthodb_link'))

read_busco_misses <- function(species_name) {
  x <- readLines(glue(busco_dir, "{species_name}.txt"))
  return(x[!startsWith(x, '#')])
}

busco_misses <- lapply(color_schema_species$species_abrv, read_busco_misses)
names(busco_misses) <- color_schema_species$species_abrv

busco_misses[c('astCal', 'troMau')]

find_overlapping_misses <- function(list_of_misses) {
  overlap_hits <- list_of_misses[[1]]
  for (i in seq(2, length(list_of_misses))) {
    overlap_hits <- intersect(overlap_hits, list_of_misses[[i]])
  }
  return(overlap_hits)
}

find_overlapping_misses(busco_misses)

busco_lookup %>% 
  filter(busco_id %in% find_overlapping_misses(busco_misses)) %>% 
  write_csv("./local/missed_BUSCO/busco_misses_all.csv")

busco_lookup %>% 
  filter(busco_id %in% find_overlapping_misses(busco_misses[c('astCal', 'mayZeb', 'troMau', 'aulStu', 'rhaChi')])) %>% 
  write_csv("./local/missed_BUSCO/busco_misses_pacbio.csv")

busco_lookup %>% 
  filter(busco_id %in% find_overlapping_misses(busco_misses[c('otoArg', 'copChr', 'rhaChi2')])) %>% 
  write_csv("./local/missed_BUSCO/busco_misses_ont.csv")
```

Upset plot of missed hits
```{r busco_misses, fig.width=6, fig.height=3.0}
df_misses <- mapply(
  function(x, y) tibble(busco_id=x, sample=y), 
  busco_misses, color_schema_species$species_abrv, 
  SIMPLIFY=FALSE) %>% 
  bind_rows() %>% 
  add_column(value=TRUE) %>% 
  pivot_wider(names_from = "sample", values_from = "value") %>% 
  replace(is.na(.), FALSE)

df_misses_agg <- df_misses %>% 
  pivot_longer(-busco_id, names_to='species') %>% 
  filter(value) %>% 
  select(-value) %>% 
  group_by(busco_id) %>% 
  summarise(species_w_missed_hits = list(species), n_species=n())

df_misses_agg %>% filter(n_species > 5)

df_misses_agg %>% 
  ggplot(aes(x = species_w_missed_hits, fill=n_species)) +
  geom_bar(fill='#88b821', width=.4) +
  geom_text(aes(label = ..count..), size=1.8, vjust=-0.8, 
            position = position_dodge(width = .9), stat = "count") +
  scale_x_upset(n_intersections=100) +
  xlab(NULL) +
  scale_y_continuous("Undetected BUSCOs", expand = expansion(mult = c(.05, .1)) ) +
  theme_combmatrix(combmatrix.panel.line.size = 0, combmatrix.panel.point.color.fill = "#397dbd")
```
