---
title: "Allele and node sharing between nonreference samples"
output: html_document
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  cache = TRUE, cache.path = '.cache_knitr/allele_sharing/',
  fig.width = 7, fig.height = 4, dpi = 350,
  fig.path='figure/2c-allele_sharing/',
  dev = c('png', 'pdf')
)
```

In this document, we are calculating the degree of allele sharing between the assemblies in the multiassembly graph.

## Setup
```{r, cache=FALSE}
library(scales)
library(factoextra)
library(reactable)
library(ComplexHeatmap)
library(tibble)
library(ggplot2)
library(readr)
library(tidyr)
library(dplyr)
library(glue)
library(ggcorrplot)
library(ggpubr)
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")
```

## Data import
We import the data frame that has been preprocessed in `1-preprocess.Rmd`.
```{r}
# note: n_paths_categorised, n_paths_bio_categorised, n_uncalled are categorical variables
df <- read_csv("./local/malawi_haplochromines-summary-preprocessed.csv",
  col_types = paste0('cciiiiicciiicilliidc', paste0(rep('ciicic', num_of_genomes), collapse=''), collapse=''),
  na = c("", "NA", ".")
)

# node lengths
df_node_lengths <- readr::read_delim(paste0(dir_sv, "malawi_haplochromines-graph.nodelengths"), col_names=c('node', 'len'))

glue("Read in preprocessed dataframe with {nrow(df)} bubbles.")
```

## Allele sharing between samples
Note that with the `replace_na` we put, it means we only include sites with allelic information for all samples.
```{r}
df_share <- df %>% 
  select(id, ends_with("_path"), ends_with("_len")) %>% 
  pivot_longer(-id, names_to = c("species", ".value"), names_pattern = "(.*)_(.*)") %>%
  drop_na() %>% 
  add_column(value = TRUE) %>% 
  pivot_wider(names_from=species) %>% 
  mutate_at(color_schema_species$species_abrv, ~replace_na(., FALSE))

df_share_agg <- df_share %>% 
  filter(astCal == FALSE) %>% 
  pivot_longer(-c(id, path, len), names_to='species') %>% 
  filter(value) %>% 
  select(- value) %>% 
  group_by(id, path, len) %>% 
  summarise(species_w_allele = list(species), n_species=n())
```

This is an UpSet plot showing the number of non-reference alleles / events shared across the non-backbone samples. It suggests that most alleles are singleton events.
```{r shared_allele_upset_count, fig.width=6.2, fig.height=4}
df_share_agg %>% 
  ggplot(aes(x = species_w_allele, fill=n_species)) +
  geom_bar(fill='#ad5090', 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=25) +
  xlab(NULL) +
  scale_y_continuous("Bubbles w shared nonreference SV", expand = expansion(mult = c(.05, .1)) ) +
  theme_combmatrix(combmatrix.panel.line.size = 0, combmatrix.panel.point.color.fill = "#d69d0d")
```

This shows the cumulative length of the non-reference alleles found in the other assemblies compared to the backbone. **I think this number is highly misleading, and should not be used.**
```{r shared_allele_cumul_length, fig.width=6.2, fig.height=3.5}
df_share_agg$species_w_allele_joined <- sapply(
  df_share_agg$species_w_allele, function(x) paste0(x, collapse=',')
)

df_share_agg2 <- df_share_agg %>% 
  group_by(species_w_allele_joined) %>% 
  summarise(n_alleles = n(), total_len = sum(len)) %>% 
  arrange(desc(total_len), species_w_allele_joined) %>% 
  mutate(species_w_allele_joined = factor(species_w_allele_joined, levels=rev(species_w_allele_joined)))

# show top few hits
df_share_agg2[1:20,] %>% 
  ggplot(aes(x = species_w_allele_joined, y = total_len/10**6)) +
  geom_bar(width=.75, color='black', fill='#6dc786', stat='identity') +
  geom_text(aes(label = format(round(total_len/10**6, 2), nsmall=2) ), 
            size=2.8, vjust=0.5, hjust=-0.1,
            position = position_dodge(width = .9), stat = "identity") +
  xlab(NULL) +
  scale_y_continuous(
    "Total length of non-reference alleles (Mbp)", 
    expand = expansion(mult = c(.05, .15))) +
  coord_flip()
```

Number of bubbles where there are private singleton nonreference nodes.
```{r}
glue("Bubbles with singleton non-references nodes: {count}",
  count = df_share_agg %>% filter(n_species==1) %>% .$id %>% unique() %>% length()
)
```

Number of bubbles where there are non-reference nodes that only appear in both Rhamphochromis.
```{r}
glue("Bubbles with non-references nodes shared by 2 Rhampho only: {count}",
  count = df_share_agg %>% filter(species_w_allele_joined == 'rhaChi,rhaChi2')  %>% .$id %>% unique() %>% length()
)
```


Calculate the length of the astCal alleles at the sites where all the nonreference alleles are the same. Note that at these bubbles, it does not necessarily mean that all the backbone nodes are unique to calliptera, because it is still possible for the nonreference samples to share those nodes at multi-allelic sites, based on our definition. However, the **count of sites** can be interpreted as bubbles where the backbone has private sequence for itself.
```{r}
glue("Bubbles with non-references nodes shared by all 7 but not in calliptera: {count}",
  count = df_share_agg %>% filter(n_species == 7) %>% .$id %>% unique() %>% length()
)
```

## Node sharing
The allele sharing is slightly misleading, because it may count some nodes multiple times if they show up on multiple alleles. I decide that it is worth doing a version for just the non-reference nodes.
```{r}
df_share_nodes <- df %>% 
  select(id, ends_with("_path"), ends_with("_len")) %>% 
  pivot_longer(-id, names_to = c("species", ".value"), names_pattern = "(.*)_(.*)") %>%
  drop_na() %>% 
  filter(path != "*") %>% 
  mutate(path_sub = gsub(">|<", ",", path)) %>% 
  mutate(path_sub = gsub("^,", "", path_sub)) %>% 
  select(species, path_sub) %>% 
  separate_longer_delim(path_sub, ',') %>% 
  group_by(path_sub, species) %>% 
  summarise(value = n()) %>% 
  mutate(value = TRUE) %>% 
  pivot_wider(names_from=species) %>% 
  replace(is.na(.), FALSE) %>% 
  rename(node = path_sub)

df_share_nodes_agg <- df_share_nodes %>% 
  filter(astCal == FALSE) %>%   # remove non-reference nodes
  pivot_longer(-node, names_to='species') %>% 
  #####
  # reorder the species (doesn't help the UpSet, but needed for cumulative)
  mutate(species = factor(species, levels=color_schema_species$species_abrv[2:8])) %>% 
  arrange(node, species) %>% 
  mutate(species = as.character(species)) %>% 
  #####
  filter(value) %>% select(-value) %>% 
  group_by(node) %>% summarise(species_w_node = list(species), n_species=n())
```

This is an UpSet plot showing the number of non-reference nodes shared across the non-backbone samples. It suggests that most nodes are private to one species.
```{r shared_node_upset_count, fig.width=6.2, fig.height=4}
df_share_nodes_agg %>% 
  ggplot(aes(x = species_w_node, fill=n_species)) +
  geom_bar(fill='#5f6de8', 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=25) +
  xlab(NULL) +
  scale_y_continuous("Shared nonreference nodes", expand = expansion(mult = c(.05, .1)) ) +
  theme_combmatrix(combmatrix.panel.line.size = 0, combmatrix.panel.point.color.fill = "#b35300")
```

This shows the cumulative length of the non-reference nodes found in the other assemblies compared to the backbone.
```{r shared_node_cumul_length, fig.width=6.2, fig.height=3.5}
df_share_nodes_agg$species_w_node_joined <- sapply(
  df_share_nodes_agg$species_w_node, function(x) paste0(x, collapse=',')
)

# aggregate and sum up node length by groups of species
df_share_nodes_agg2 <- df_share_nodes_agg %>%
  left_join(df_node_lengths) %>% 
  group_by(species_w_node_joined) %>% 
  summarise(n_nodes = n(), total_len = sum(len)) %>% 
  arrange(desc(total_len), species_w_node_joined) %>% 
  mutate(species_w_node_joined = factor(species_w_node_joined, levels=rev(species_w_node_joined)))

# show top few hits
df_share_nodes_agg2[1:20,] %>% 
  ggplot(aes(x = species_w_node_joined, y = total_len/10**6)) +
  geom_bar(width=.75, color='black', fill='#6dc786', stat='identity') +
  geom_text(aes(label = format(round(total_len/10**6, 2), nsmall=2) ), 
            size=2.8, vjust=0.5, hjust=-0.1,
            position = position_dodge(width = .9), stat = "identity") +
  xlab(NULL) +
  scale_y_continuous(
    "Total length of non-reference segments (Mbp)", 
    expand = expansion(mult = c(.05, .15))) +
  coord_flip()
```

Sanity checks:

* All the non-reference nodes should be no more than `291,473,128` bp spread over `254,301` segments. 
* All the variable nodes should be no more than `413,047,161` bp spread over `448,040` segments.

We might have filtered out some when removing sites with insufficient coverage across the assemblies, so numbers might be lower, but not higher.
```{r}
glue("All non-references nodes: ",
  "{nrow(tmp)} segments: {sum(tmp$len)} bp",
  tmp = df_node_lengths %>% filter(node %in% df_share_nodes_agg$node)
)

glue("All variable nodes used: ",
  "{nrow(tmp)} segments: {sum(tmp$len)} bp",
  tmp = df_node_lengths %>% filter(node %in% df_share_nodes$node)
)
```

Interesting stats worth quoting:

Non-reference nodes that are singletons, only appear in one assembly.
```{r}
glue("Singleton non-references nodes: ",
  "{sum(tmp$n_nodes)} segments: {sum(tmp$total_len)} bp",
  tmp = df_share_nodes_agg2 %>% filter(species_w_node_joined %in% color_schema_species$species_abrv)
)
```

Non-reference nodes that only appear in both Rhamphochromis.
```{r}
glue("Non-references nodes shared by both Rhamphochromis: ",
  "{sum(tmp$n_nodes)} segments: {sum(tmp$total_len)} bp",
  tmp = df_share_nodes_agg2 %>% filter(species_w_node_joined == "rhaChi,rhaChi2")
)
```


Non-reference nodes missing entirely in calliptera.
```{r}
glue("Non-references nodes shared by others, but not present in calliptera: ",
  "{sum(tmp$n_nodes)} segments: {sum(tmp$total_len)} bp",
  tmp = df_share_nodes_agg2 %>% filter(species_w_node_joined == "mayZeb,troMau,aulStu,otoArg,copChr,rhaChi,rhaChi2")
)
```

Backbone nodes that only appear in calliptera.
```{r}
df_share_nodes %>% 
  filter(astCal == TRUE) %>% 
  filter(aulStu==F, mayZeb==F, otoArg==F, copChr==F, rhaChi==F, rhaChi2==F, troMau==F) %>% 
  select(node) %>% 
  left_join(df_node_lengths) %>%
  ungroup() %>% 
  summarise(count = n(), total_len=sum(len))
```

