---
title: "Discovering SVs in the Malawi cichlids pan-genome graph"
output: html_document
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  cache = TRUE, cache.path = '.cache_knitr/analyse_sv_properties/',
  fig.width = 7, fig.height = 4, dpi = 350,
  fig.path='figure/2a-analyse_sv_properties/',
  dev = c('png', 'pdf')
)
```

In this document, we are utilising various metrics to understand the properties of the Lake Malawi cichlid pan-genome graph built using `minigraph` using default parameters. We aim to get an idea of the openness of the pangenome, how many sites of structural variation (SV) there are and how widespread/private alleles are.

## 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", ".")
)

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

summary(select(df, n_segments:lenbio_minmaxratio, -inv, -segment_list) )
```


## Principal component analysis
We want to run a PCA using the allele lengths at all the SVs. We extract the columns of the allele lengths, and calculate the `variance` row-wise as a new column.
```{r}
sv_lengths_var <- df %>% 
  select(ends_with('_len')) %>%
  rowwise() %>%
  mutate(variance = c_across(ends_with('_len')) %>% var()) %>%
  ungroup()

# rename the column names to be shorter
names(sv_lengths_var) <- gsub('_len', '', names(sv_lengths_var))
```

We exclude loci where all the alleles are of the same length (i.e. variance of zero). For this analysis to work, we have to filter out the bubbles where there are uncalled alleles (variance = NA) and those where alleles are all the same length (variance = 0). Both of these cause problems with PCA.
```{r}
# select those where the variance != 0 
# those could be substitutions / inversions
# or all species have the same allele
sv_lengths_nonzerovar <- sv_lengths_var %>%
  filter(variance != 0) %>%
  select(-variance)

glue('Number of length-variables SVs: {nrow(sv_lengths_nonzerovar)}')
```

We also like to exclude loci where the length difference between the alleles is not huge. This is because we are making a Z-score heatmap later on, so these variants may seem more different than they actually are.
```{r allele_len_diff_hist, fig.height=3}
max_allele_lengths <- apply(sv_lengths_nonzerovar, 1, max)
min_allele_lengths <- apply(sv_lengths_nonzerovar, 1, min)
allele_len_diff <- max_allele_lengths - min_allele_lengths

ggplot() + 
  geom_histogram(aes(x = allele_len_diff), fill='gray70', color='black') + 
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  xlab("Length difference between longest and shortest theoretical path") +
  ylab("Number of bubbles")

# print some stats of the length difference
summary(allele_len_diff)

# standard deviation
sd(allele_len_diff)

# standard error
sd(allele_len_diff)/sqrt(length(allele_len_diff))
```

Filter to keep those above the 25th percentile.
```{r}
sv_lengths_subset <- sv_lengths_nonzerovar[allele_len_diff > quantile(allele_len_diff, probs = .25), ]
glue('Number of length-variables SVs above 25th percentile: {nrow(sv_lengths_subset)}')
```

```{r pca_species, fig.width=4, fig.height=3.2}
res.pca <- prcomp(t(sv_lengths_subset), scale = TRUE)
fviz_eig(res.pca, choice='variance')
fviz_pca_ind(res.pca, 
  col.ind = "cos2", repel = TRUE,
  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")
)
```

We also plot the Pearson's correlation between the samples. We did a log transform of the data, so that there is more of a linear relationship and normal-ish distribution (assumption of Pearson correlation).
```{r corr_sv_sanity_checks}
ggplot() + 
  geom_point(aes(x = sv_lengths_var$astCal + 1, y = sv_lengths_var$rhaChi + 1)) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  xlab("allele length in astCal") +
  ylab("allele length in rhaChi")

ggplot() + 
  geom_histogram(aes(x = sv_lengths_var$astCal + 1), fill='gray70', color='black', bins = 15) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  xlab("allele length in astCal") 
```

We make some correlation plots using the `r nrow(df %>% filter(n_uncalled == 0))` SVs where every species has an allele called.
```{r corr_sv_pearson, fig.width=5, fig.height=5}
sv_log_transformed <- log10(drop_na(sv_lengths_var)[1:num_of_genomes] + 1)
corr <- cor(sv_log_transformed, method='pearson')
round(corr, 2)

ggcorrplot(corr, method = "square", lab = TRUE, 
  outline.color = 'gray10', digits = 2, hc.order = FALSE) +
  scale_fill_gradient2(midpoint=0.5, limits=c(0,1), low='#6D9EC1', mid='#E5DEDA', high='#E46726') +
  theme_minimal() +
  labs(fill='correlation') +
  xlab(NULL) + ylab(NULL) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        text = element_text(size=15),
        legend.position='bottom',
        legend.title=element_text(size=rel(0.6)),
        legend.text=element_text(size=rel(0.5)) )
```

Now here are some plots to show that the Rhamphochromis are closer to each other than expected. For every possible pair of species, we want to obtain the shortest and longest allele at each locus, and calculate their length difference. This analysis is obviously constrained to those sites where the species pair have different alleles.
```{r}
build_allele_len_diff_df <- function(x) {
  columns_of_species_pair <- sv_lengths_nonzerovar[x]
  
  # get abbrevations of the two species
  species_abrv_pair <- color_schema_species$species_initials[
    match(x, color_schema_species$species_abrv)
  ]
  
  output_df <- data.frame(
    species_pair = paste0(species_abrv_pair[1], '-', species_abrv_pair[2]),
    species1 = species_abrv_pair[1],
    species2 = species_abrv_pair[2],
    short_len = apply(columns_of_species_pair, MARGIN=1, min),
    long_len  = apply(columns_of_species_pair, MARGIN=1, max)
  ) 
  return(output_df)
}


df_lendiff_pairs <- apply(
  combn(color_schema_species$species_abrv, 2), MARGIN=2, build_allele_len_diff_df)  %>% 
  bind_rows() %>% 
  mutate(len_diff = long_len - short_len)

df_lendiff_pairs_duplicated <- bind_rows(
  df_lendiff_pairs,
  df_lendiff_pairs %>% rename(species1 = species2, species2 = species1)
)

head(df_lendiff_pairs_duplicated)
```

We count the number of sites 
```{r rhampho_closer_diffsites_count, fig.width=5, fig.height=5}
df_lendiff_pairs_summary <- df_lendiff_pairs_duplicated %>% 
  filter(len_diff > 0) %>% 
  group_by(species1, species2) %>% 
  summarise(count = n(), median = median(len_diff), mean = mean(len_diff)) %>% 
  bind_rows(data.frame(
      species1 = color_schema_species$species_initials,
      species2 = color_schema_species$species_initials,
      count = 0, median=0, mean=0)
  ) %>% 
  mutate(
    species1 = factor(species1, 
      levels = color_schema_species$species_initials, labels = color_schema_species$species_abrv),
    species2 = factor(species2,
      levels = color_schema_species$species_initials, labels = color_schema_species$species_abrv)
  )

ggplot(df_lendiff_pairs_summary, aes(species1, species2, fill=count/1000)) +
  geom_tile(color='gray10') +
  geom_text(aes(label = format(round(count/1000, 1), nsmall = 1) ), 
            df_lendiff_pairs_summary %>% filter(species1!=species2), color='black') +
  scale_fill_distiller('number of\nsites with\ndifferent alleles\n(x 10^3)', palette = 'Purples') +
  theme(aspect.ratio = 1) +
  xlab(NULL) + ylab(NULL) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), 
        text = element_text(size=15), aspect.ratio = 1,
        legend.position='bottom',
        legend.title=element_text(size=rel(0.6)),
        legend.text=element_text(size=rel(0.5)) )
```

There is much fewer sites between the two Rhamphochromis where the alleles are different, even though the average length difference is quite similar.
```{r}
df_lendiff_pairs_summary %>% reactable()

df_lendiff_pairs_summary %>% 
  filter(!(species1 == 'rhaChi1' & species2 == 'rhaChi2') ) %>% 
  filter(count != 0) %>% 
  summary(digits=10)
```

We also plot the distribution of the allele length difference using rhaChi, rhaChi2, mayZeb (mbuna) and aulStu (benthics) as a point of comparison. We constrain only to places where the alleles are different. We can notice that the pair of Rhamphochromis are much closer to each other genetically than you would expect even when considering species from the same ecological clade.
```{r rhampho_closer_lendiff_hist, fig.width=5, fig.height=3}
df_lendiff_pairs_duplicated %>% 
  filter(species1 %in% c('rc', 'rc2', 'ac'), len_diff>0) %>% 
  ggplot() +
  geom_freqpoly(aes(len_diff + 1, color=species2)) +
  facet_wrap(.~factor(species1, levels = c('rc', 'rc2', 'ac'), 
    labels = c('rhaChi', 'rhaChi2', 'astCal')), ncol=3) +
  scale_color_manual(values=color_schema_species$color, 
    breaks=color_schema_species$species_initials, labels=color_schema_species$species_abrv) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  xlab('length difference between alleles') + labs(color = 'compared species') +
  #ggtitle('Allele length differences between species (excl sites w same allele)') +
  theme(legend.position='top') 
```

## PCA with presence-absence matrix
Constrained to only bi-allelic and tri-allelic sites.
```{r}
source("modules/create_presence_absence_matrix.R")
pa_matrix_pca <- create_pa_matrix(
  df %>% filter(n_paths_bio %in% c(2,3)) %>% select(ends_with('_path')) %>% drop_na()
)
colnames(pa_matrix_pca) <- color_schema_species$species_abrv 
head(pa_matrix_pca)
```

Construct PCA with PAV matrix.
```{r pca_species_pav, fig.width=4, fig.height=3.2}
res.pca2 <- prcomp(t(pa_matrix_pca), scale = TRUE)
fviz_eig(res.pca2, choice='variance')
fviz_pca_ind(res.pca2,
  col.ind = "cos2", repel = TRUE,
  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")
)
```


## Heatmap and phylogeny dendrograms
We construct heatmaps and dendrograms for a few thousand randomly chosen loci (we can't plot all because this dataset is too large)
```{r}
sample_idx = sample(nrow(sv_lengths_subset), 5000)
```

```{r dendrogram, fig.width=10, fig.height=3}
fig_heatmap <- sv_lengths_subset %>% 
  .[sample_idx, ] %>%
  t() %>% scale() %>%  # z score
  Heatmap(
    name="SV length\n(Z-score)",
    col=circlize::colorRamp2(c(-2.5, 0, 2.5), c("#A6611A", "#f7f7f7", "#018571")),
    column_title = NULL, row_title = "Sample",
    height = nrow(.)*unit(8, "mm"),
    show_heatmap_legend = TRUE, show_column_names=FALSE
  )
print(fig_heatmap)

# log2 version, with no scaling
fig_heatmap <- log10(sv_lengths_subset+1) %>% 
  .[sample_idx, ] %>%
  t() %>% 
  Heatmap(
    name="SV length\n(log10)",
    col=circlize::colorRamp2(c(0, 2, 4), c("#A6611A", "#f7f7f7", "#018571")),
    column_title = NULL, row_title = "Sample",
    height = nrow(.)*unit(8, "mm"),
    show_heatmap_legend = TRUE, show_column_names=FALSE
  )
print(fig_heatmap)
```

We plot the same set of loci, but with the lengths subtracted against calliptera, to get a relative allele length. Note that the ordering of the columns is not the same.
```{r dendrogram_relative, fig.width=10, fig.height=3}
sv_lengths_subset_subtracted <- sv_lengths_subset - sv_lengths_subset[['astCal']]
fig_heatmap <- sv_lengths_subset_subtracted %>% 
  select(-astCal) %>% 
  .[sample_idx, ] %>%
  t() %>% scale(center=FALSE) %>% # do not center
  Heatmap(
    name="SV length\n(Z-score)",
    col=circlize::colorRamp2(c(-2, 0, 2), c("#A6611A", "#f7f7f7", "#018571")),
    column_title = NULL, row_title = "Sample",
    height = nrow(.)*unit(8, "mm"),
    show_heatmap_legend = TRUE, show_column_names=FALSE
  )
print(fig_heatmap)
```

No scaling, just log2 transformation.
```{r dendrogram_relative_log10, fig.width=10, fig.height=3}
sv_lengths_subset_subtracted_log10 <- log10(sv_lengths_subset+1) - log10(sv_lengths_subset[['astCal']]+1)
fig_heatmap <- sv_lengths_subset_subtracted_log10 %>% 
  select(-astCal) %>% 
  .[sample_idx, ] %>%
  t() %>% 
  Heatmap(
    name="SV length\n(log10)",
    col=circlize::colorRamp2(c(-4, 0, 4), c("#A6611A", "#f7f7f7", "#018571")),
    column_title = NULL, row_title = "Sample",
    height = nrow(.)*unit(8, "mm"),
    show_heatmap_legend = TRUE, show_column_names=FALSE
  )
print(fig_heatmap)
```

From the heatmaps, we get an early impression of the allele frequency landscape. It seems like for a substantial proportion of SV loci, the minor allele is only present in one species, regardless of whether it's an insertion (orange) or deletion (purple). These are known as *singletons*.

## Type of alleles
### Number of alleles
Counting the number of alleles observed in total, and how many on average there are per bubble.
```{r}
df %>% group_by(n_paths_bio) %>% 
  summarise(count = n()) %>% 
  mutate(n_paths_bio_counted = n_paths_bio * count) %>% 
  summarise(total_bubbles = sum(count), total_alleles = sum(n_paths_bio_counted)) %>% 
  mutate(alleles_per_locus = total_alleles / total_bubbles)
```

### Calculations and data formating
We want to know the type of structural variation (insertion, deletion, substitution/inversion), as well as their frequencies in the five species studied. There is some substantial formating that has to be done here.

Currently, the allele path and path length is in wide format (each species is a column). We'll like to convert this to long format, and add some additional columns about the bubble using `left_join`. Each row would correspond to the path taken  by a species at a particular bubble.

There are rows where the path is `NA`, meaning that the allele is not known. We would discard these rows, keeping only the called paths for that bubble. For these, the value of `n_uncalled` should be greater than 0.
```{r}
# convert into tidy format
# add columns about the reference allele and relative length
tmp_df_allele <- df %>% 
  select(id, n_paths_bio, n_uncalled, ends_with(c('_path', '_len'))) %>% 
  pivot_longer(cols = ends_with(c('_path', '_len')),
    names_to = c('species', '.value'),
    names_pattern = "(.*)_(.*)") %>%
  left_join(select(df, id, `astCal_path`, `astCal_len`), by='id') %>% 
  dplyr::rename(path_ref = `astCal_path`, len_ref = `astCal_len`) %>% 
  mutate(ref = if_else(path == path_ref, TRUE, FALSE),
         len_rel = len - len_ref) %>% 
  drop_na(path)

tmp_df_allele
```

Making use of the `id` (genomic locus) and `path` (allele) in `group_by`, we can tally the allele frequency. Then making use of this information, for each locus, we annotate whether each path is a minor allele and count the number of alleles observed in that locus (`n_alleles` --- note this is different from `n_paths`, which is the number of theoretical paths).
```{r}
# count the number of each allele for each path
# rejoin certain columns needed to determine type of mutation
tmp_df_allele2 <- tmp_df_allele %>% 
  group_by(id, path) %>%  # this is needed because the same allele can show up across multiple loci (e.g. *)
  count(path) %>% 
  group_by(id) %>% 
  mutate(minor = if_else(n == min(n), TRUE, FALSE) ) %>% 
  left_join(select(tmp_df_allele, id, n_paths_bio, n_uncalled, path, len, len_rel, len_ref, ref) %>% unique(), 
            by = c('id', 'path'))

tmp_df_allele2
```

For each allele, we determine its type and add that information as a new column.
```{r}
type_annotations <- rbind(
  tmp_df_allele2 %>% filter(ref==FALSE, len_rel>0, len_ref==0) %>% add_column(allele_type='complete insertion'),
  tmp_df_allele2 %>% filter(ref==FALSE, len_rel>0, len_ref>0)  %>% add_column(allele_type='partial insertion'),
  tmp_df_allele2 %>% filter(ref==FALSE, len_rel<0, len_ref==-len_rel) %>% add_column(allele_type='complete deletion'),
  tmp_df_allele2 %>% filter(ref==FALSE, len_rel<0, len_ref!=-len_rel) %>% add_column(allele_type='partial deletion'),
  tmp_df_allele2 %>% filter(ref==FALSE, len_rel==0) %>% add_column(allele_type='substitution/inversion')
)
tmp_df_allele3 <- left_join(tmp_df_allele2, select(type_annotations, id, path, allele_type), by = c('id', 'path'))
```

Also useful is to summarise the types of structural variation present at a genomic locus. This locus level information is then added into the main dataframe.

* For bi-allelic loci, this simply corresponds to the type of allele, and could be one of either: complete deletion, partial deletion, complete insertion, partial insertion or inversion/substitution.
* For multi-allelic loci, we need to see what types of alleles are found there, and can summarise that information as either all deletions, all insertions or mixed.
```{r}
df_loci_biallele <- tmp_df_allele3 %>% 
  filter(n_paths_bio == 2) %>% filter(ref==FALSE) %>% 
  select(id, allele_type) %>% 
  dplyr::rename(locus_type = allele_type)

df_loci_multiallele <- tmp_df_allele3 %>% 
  filter(n_paths_bio > 2) %>% filter(ref==FALSE) %>% 
  select(id, allele_type) %>% 
  group_by(id) %>% 
  summarise(
    locus_type = if_else(all(allele_type %in% c('partial deletion', 'complete deletion')), 'all deletions',
                         if_else(all(allele_type %in% c('partial insertion', 'complete insertion')), 'all insertions', 
                                 if_else(all(allele_type == 'substitution/inversion'), 'all substitutions/inversions', 'mixed')))
)

df_allele <- left_join(tmp_df_allele3, rbind(df_loci_biallele, df_loci_multiallele), by='id')

rm(df_loci_biallele, df_loci_multiallele)
rm(tmp_df_allele, tmp_df_allele2, tmp_df_allele3)

df_allele
```

We want to refactor certain columns to have a certain ordering of categories for the plots later on.
```{r}
df_allele$allele_type <- factor(df_allele$allele_type, 
  levels = c('complete deletion', 'partial deletion', 'partial insertion', 'complete insertion', 'substitution/inversion')
)

df_allele$locus_type <- factor(df_allele$locus_type, 
  levels = c('complete deletion', 'partial deletion', 'partial insertion', 'complete insertion',
    'substitution/inversion', 'all deletions', 'mixed', 'all insertions')
)
```


### Plots
We plot a frequency histogram of the number of alleles / biological paths observed at the genomic loci. Contrast with the plot from the previous section detailing the number of theoretical paths. Note that a version of this plot has already been generated before in `1-preprocess.Rmd`.
```{r}
df_allele %>% filter(ref == TRUE) %>% group_by(n_paths_bio) %>% 
  summarise(count = n()) %>% 
  mutate(percentage = round(count/sum(count) * 100, 2))
```

Version with only bubbles with no more than 3 uncalled alleles.
```{r}
df_allele %>% filter(ref == TRUE, n_uncalled <= 3) %>% group_by(n_paths_bio) %>% 
  summarise(count = n()) %>% 
  mutate(percentage = round(count/sum(count) * 100, 2))
```

Here's a plot surmising the types of mutations we see at the SV loci. First version is for all bubbles, second version is for only those with no more than 3 uncalled alleles (They are the same though).
```{r sv_type_loci_summary, fig.height=3.5}
df_allele %>% filter(n_paths_bio >= 2, ref==TRUE) %>% 
  mutate(allelism = if_else(n_paths_bio==2, 'bi-allelic', 'multi-allelic')) %>% 
  ggplot(aes(x=locus_type)) + 
  facet_wrap(vars(allelism), scales = 'free_x') +
  geom_bar() +
  scale_y_continuous(expand = expansion(mult = c(.05, .1)) ) +
  geom_text(aes(label = ..count..), vjust = -0.5, 
            position = position_dodge(width = .9), stat = "count") +
  xlab(NULL) +  ylab("number of genomic loci") +
  theme(legend.position = 'none', 
    axis.text.x = element_text(angle = 35, vjust = 1, hjust=1)
  )

df_allele %>% filter(n_paths_bio >= 2, ref==TRUE, n_uncalled <= 3) %>% 
  mutate(allelism = if_else(n_paths_bio==2, 'bi-allelic', 'multi-allelic')) %>% 
  ggplot(aes(x=locus_type)) + 
  facet_wrap(vars(allelism), scales = 'free_x') +
  geom_bar() +
  scale_y_continuous(expand = expansion(mult = c(.05, .1)) ) +
  geom_text(aes(label = ..count..), vjust = -0.5, 
            position = position_dodge(width = .9), stat = "count") +
  xlab(NULL) +  ylab("number of genomic loci") + 
  theme(axis.text.x = element_text(angle = 35, vjust = 1, hjust=1) )
```

Next, for some allelic level analysis. This tabulates the type of alleles observed across all bubbles.
```{r sv_type_loci_summary_color, fig.width=4, fig.height=5}
df_allele %>% filter(n_paths_bio >= 2, ref==FALSE) %>% 
  #mutate(allele_type = if_else(is.na(allele_type), 'inversion/substitution', as.character(allele_type))) %>% 
  mutate(allele_type = factor(allele_type, levels=c('complete deletion', 'partial deletion', 'partial insertion', 'complete insertion', 'substitution/inversion'))) %>% 
  ggplot(aes(x=allele_type, fill=allele_type, color=allele_type)) + 
  geom_bar() +
  scale_y_continuous(expand = expansion(mult = c(.05, .1)) ) +
  scale_fill_manual(values=alpha(c(ggpubr::get_palette('BrBG', k=4), 'gray10'), 0.7)) +
  scale_color_manual(values=c(ggpubr::get_palette('BrBG', k=4), 'gray10')) +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count") +
  xlab(NULL) +
  ylab("number of alleles") +
  theme(axis.text.x = element_text(angle = 35, vjust = 1, hjust=1), legend.position='none')

df_allele %>% filter(n_paths_bio >= 2, n_uncalled == 0, ref==FALSE) %>% 
  #mutate(allele_type = if_else(is.na(allele_type), 'inversion/substitution', as.character(allele_type))) %>% 
  mutate(allele_type = factor(allele_type, levels=c('complete deletion', 'partial deletion', 'partial insertion', 'complete insertion', 'substitution/inversion'))) %>% 
  ggplot(aes(x=allele_type, fill=allele_type, color=allele_type)) + 
  geom_bar() +
  scale_y_continuous(expand = expansion(mult = c(.05, .1)) ) +
  scale_fill_manual(values=alpha(c(ggpubr::get_palette('BrBG', k=4), 'gray10'), 0.7)) +
  scale_color_manual(values=c(ggpubr::get_palette('BrBG', k=4), 'gray10')) +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count") +
  xlab(NULL) +
  ylab("number of alleles") +
  theme(axis.text.x = element_text(angle = 35, vjust = 1, hjust=1), legend.position='none')
```

Poster version (with no substitutions)
```{r sv_type_loci_summary_color-poster, fig.width=2.6, fig.height=2.8}
df_allele %>% 
  filter(n_paths_bio >= 2, n_uncalled == 0, ref==FALSE) %>% 
  filter(allele_type != "substitution/inversion") %>% 
  mutate(allele_type = factor(allele_type, levels=c('complete deletion', 'partial deletion', 'partial insertion', 'complete insertion'))) %>%
  ggplot(aes(x=allele_type, fill=allele_type, color=allele_type)) + 
  geom_bar(width = 0.8) +
  scale_y_continuous(expand = expansion(mult = c(.05, .25)) ) +
  scale_fill_manual(values=alpha(c(ggpubr::get_palette('BrBG', k=4), 'gray10'), 0.8)) +
  scale_color_manual(values=c(ggpubr::get_palette('BrBG', k=4), 'gray10')) +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count") +
  xlab(NULL) +
  ylab("number of alleles") +
  theme_linedraw() +
  theme(axis.text.x = element_text(angle = 35, vjust = 1, hjust=1), legend.position='none')
```

Sized for paper.
```{r sv_type_loci_summary_color-PAPER, fig.width=3.2, fig.height=2.8}
df_allele %>% filter(n_paths_bio >= 2, ref==FALSE) %>% 
  #mutate(allele_type = if_else(is.na(allele_type), 'inversion/substitution', as.character(allele_type))) %>% 
  mutate(allele_type = factor(allele_type, levels=c('complete deletion', 'partial deletion', 'partial insertion', 'complete insertion', 'substitution/inversion'))) %>% 
  ggplot(aes(x=allele_type, fill=allele_type, color=allele_type)) + 
  geom_bar() +
  scale_y_continuous(expand = expansion(mult = c(.05, .25)) ) +
  scale_fill_manual(values=alpha(c(ggpubr::get_palette('BrBG', k=4), 'gray10'), 0.7)) +
  scale_color_manual(values=c(ggpubr::get_palette('BrBG', k=4), 'gray10')) +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count") +
  xlab(NULL) +
  ylab("number of alleles") +
  theme(axis.text.x = element_text(angle = 35, vjust = 1, hjust=1), legend.position='none')
```

Some summary statistics for all loci (regardless of how many uncalled alleles) to quote.
```{r}
# the total number of non-reference alleles
df_allele %>% filter(ref == FALSE) %>% nrow()

# count the number of each type of allele
df_allele %>% filter(ref == FALSE) %>% .$allele_type %>% table()

# total relative length and count
df_allele %>% 
  filter(ref == FALSE) %>% 
  group_by(allele_type) %>% 
  summarise(count = n(), cumulative_lenrel = sum(len_rel), mean_lenrel = round(mean(len_rel), 2),
            min_lenrel = min(len_rel), max_lenrel = max(len_rel)) %>% 
  mutate(percentage = round(count/sum(count) * 100, 2)) %>% 
  reactable()

# recalculate using absolute length
df_allele %>% 
  filter(ref == FALSE) %>% 
  group_by(allele_type) %>% 
  summarise(count = n(), cumulative_len = sum(len), mean_len = round(mean(len), 2),
            min_len = min(len), max_len = max(len)) %>% 
  mutate(percentage = round(count/sum(count) * 100, 2)) %>% 
  reactable()

# insertion and deletions more broadly
df_allele %>% 
  filter(ref == FALSE, allele_type != "substitution/inversion") %>% 
  group_by(stringr::str_detect(allele_type, "insertion")) %>% 
    summarise(count = n(), 
              cumulative_lenrel = sum(len_rel), 
              mean_lenrel = round(mean(len_rel), 2),
              min_lenrel = round(min(len_rel), 2), 
              max_lenrel = round(max(len_rel), 2) ) %>% 
  reactable()
```

Here we plot the relative length of the non-reference allele against the reference.
```{r allele_lenrel_histogram, fig.width=4, fig.height=5.2}
df_allele %>% filter(n_paths_bio >= 2, allele_type != 'substitution/inversion', ref == FALSE) %>% 
  ggplot(aes(abs(len_rel), fill=allele_type, color=allele_type)) +
  facet_grid(rows=vars(allele_type)) +
  scale_fill_brewer("Type of allele", palette = "BrBG") +
  scale_color_brewer("Type of allele", palette = "BrBG") +
  geom_histogram(position='identity', alpha=0.7, center=0, binwidth=0.1) +
  geom_vline(xintercept = 50+1, linetype='dotted') +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  xlab("length difference of nonreference allele vs backbone") +  
  ylab("number of alleles") +
  theme(legend.position = "none")
```

Non-faceted version.
```{r allele_lenrel_histogram_all, fig.width=4.4, fig.height=2.5}
df_allele %>% filter(n_paths_bio >= 2, allele_type != 'substitution/inversion', ref == FALSE) %>% 
  ggplot(aes(abs(len_rel))) +
  geom_histogram(position='identity', fill='gray70', color='black', alpha=0.7, center=0, binwidth=0.1) +
  geom_vline(xintercept = 50+1, linetype='dotted') +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  xlab("length difference of non-reference allele vs backbone") +  
  ylab("number of alleles") +
  theme(legend.position = "none")

# print some stats
df_allele %>%
  filter(n_paths_bio >= 2, allele_type != 'substitution/inversion', ref == FALSE) %>% 
  ungroup() %>% 
  transmute(abs_len_rel = abs(len_rel)) %>% 
  summary()
```

What are the longest SVs?
```{r}
tmp_df_all_svs <- df_allele %>% 
  mutate(allele_type_parsed = stringr::str_extract(allele_type, "insertion|deletion")) %>% 
  filter(ref == FALSE) %>% 
  select(allele_type_parsed, len, len_rel) %>% 
  mutate(abs_len_rel = abs(len_rel)) 

tmp_df_all_svs %>% 
  ungroup() %>% 
  group_by(allele_type_parsed) %>% 
  summarise(
    count = n(), mean = round(mean(abs_len_rel), 2),
    q50 = median(abs_len_rel, 0.5), q90 = quantile(abs_len_rel, 0.9),
    q95 = quantile(abs_len_rel, 0.95), q99 = quantile(abs_len_rel, 0.99),
    q100 = max(abs_len_rel)
  )
```




Length deviation but only for the longest biological allele vs backbone.
```{r allele_len_rel_histogram_all_longest, fig.width=4.2, fig.height=2.2}
tmp_df <- df_allele %>% 
  filter(n_paths_bio >= 2, allele_type != 'substitution/inversion', ref == FALSE) %>%
  group_by(id) %>%
  slice_max(order_by = abs(len_rel)) %>% 
  mutate(allele_type_parsed = stringr::str_extract(allele_type, "insertion|deletion"))
  
tmp_df %>% 
  ggplot(aes(abs(len_rel))) +
  geom_histogram(position='identity', fill='gray70', color='black', alpha=0.7, center=0, binwidth=0.1) +
  geom_vline(xintercept = 50+1, linetype='dotted') +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  xlab("Length difference of longest structural variant vs backbone") +  
  ylab("Number of bubbles") +
  theme(legend.position = "none")
```

```{r allele_len_rel_histogram_all_longest-poster, fig.width=2.9, fig.height=2.8}
tmp_df %>% 
  ggplot(aes(abs(len_rel), color=allele_type_parsed)) +
  geom_freqpoly(position='identity', center=0, binwidth=0.1) +
  geom_vline(xintercept = 50+1, linetype='dotted') +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  xlab("Length deviation of\nlongest structural variant") +  
  ylab("Number of bubbles") +
  scale_color_manual(NULL, values = ggpubr::get_palette('BrBG', k=8)[c(2,7)]) +
  theme_linedraw() +
  theme(
    legend.position = c(.95, .95),
    legend.justification = c("right", "top"),
    legend.box.just = "right",
    legend.margin = margin(0, 0, 0, 0)
  )
```

The allele frequencies of the SVs at the bi-allelic sites reveal what we suspected earlier, that most types of variation are rare and often occur in only one individual. It seems like there is more insertions than deletions, but whether this is biologically true or a technical discovery bias --- I cannot tell.
```{r sv_type_allele_bi, fig.width=8, fig.height=2.5}
df_allele %>% filter(n_paths_bio == 2, allele_type != 'substitution/inversion', ref==FALSE) %>% 
  ggplot(aes(n, fill = allele_type)) + 
  facet_wrap(vars(allele_type), ncol = 4) +
  geom_bar() +
  scale_y_continuous(expand = expansion(mult = c(.05, .2)) ) +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count", size = unit(2.4, "npc")) +
  scale_fill_brewer("Type of allele", palette = "BrBG") +
  #ggtitle("Structural variation at bi-allelic sites") + 
  xlab("sample frequency (out of 8 individuals)") +  ylab("number of alleles") +
  theme(aspect.ratio = 1, legend.position = 'none')
```

The high prevalence of singletons is also seen at multi-allelic loci.
```{r sv_type_allele_multi, fig.width=8, fig.height=2.5}
df_allele %>% filter(n_paths_bio > 2, allele_type != 'substitution/inversion', ref==FALSE) %>% 
  ggplot(aes(n, fill = allele_type)) + 
  facet_wrap(vars(allele_type), ncol = 4) +
  geom_bar() +
  scale_y_continuous(expand = expansion(mult = c(.05, .2)) ) +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count", size = unit(2.4, "npc")) +
  scale_fill_brewer("Type of allele", palette = "BrBG") +
  #ggtitle("Structural variation at multi-allelic sites") + 
  xlab("sample frequency (out of 8 individuals)") +  ylab("number of alleles") +
  theme(aspect.ratio = 1, legend.position = 'none')
```

The same as above but for all the loci.
```{r sv_type_allele_all_facet, fig.width=9, fig.height=3.5}
# this is worth redoing with fill = TE_proportion later on
df_allele %>% filter(n_paths_bio >= 2, allele_type != 'substitution/inversion', ref==FALSE) %>% 
  ggplot() + 
  geom_bar(aes(n, fill = cut(n_paths_bio, breaks=c(2,3,4,num_of_genomes+1), include.lowest = TRUE) )) +
  facet_wrap(vars(allele_type), ncol = 4) +
  scale_fill_manual("Total observed paths at bubble", 
                    values = rev(ggpubr::get_palette('Blues', k=4)[2:4]),
                    labels = c('2', '3', '> 3') ) +
  scale_y_continuous(expand = expansion(mult = c(.05, .1)) ) +
  geom_text(aes(n, label = ..count..), vjust = -0.5, stat = "count", size = unit(2.7, "npc")) +
  #ggtitle("Population prevalence of SVs across all genomic loci") + 
  xlab("sample frequency (out of 8 individuals)") +  ylab("number of alleles") +
  theme(legend.position = 'top')
```

As above but not facetted.
```{r sv_type_allele_all, fig.width=4, fig.height=3.2}
df_allele %>% filter(n_paths_bio >= 2, ref==FALSE) %>% 
  ggplot() + 
  geom_bar(aes(n, fill = cut(n_paths_bio, breaks=c(2,3,4,num_of_genomes+1), include.lowest = TRUE) )) +
  scale_fill_manual("Total observed\npaths at bubble", 
                    values = rev(ggpubr::get_palette('Blues', k=4)[2:4]),
                    labels = c('2', '3', '> 3') ) +
  scale_y_continuous(labels = label_comma(), expand = expansion(mult = c(.05, .1)) ) +
  geom_text(aes(n, label = ..count..), vjust = -0.5, stat = "count", size = rel(3.5)) +
  xlab("sample frequency (out of 8 individuals)") +  ylab("number of alleles") +
  theme(
  legend.position = c(.95, .8),
    legend.justification = c("right", "top"),
    legend.box.just = "right",
    legend.margin = margin(6, 6, 6, 6),
    legend.title = element_text(size = rel(0.8)),
    legend.text = element_text(size = rel(0.8)),
    legend.key.size=unit(0.06, "npc")
  )
```

And some printed statistics for the percentage of alleles by frequency.
```{r}
# count the number of alleles based on their population frequency
df_allele %>%
  filter(ref==FALSE) %>% 
  group_by(n) %>% 
  summarise(count=n()) %>% 
  mutate(percentage = round(count/sum(count) * 100, 2))
```

## Distribution of allele lengths by frequency
Are rarer alleles longer? Maybe... there seems to be enrichment of complete indels at low frequencies, and they also seem to be longer.
```{r allele_lenrel_histogram_facet_freq, fig.width=4, fig.height=5.2}
df_allele %>% filter(n_paths_bio >= 2, allele_type != 'substitution/inversion', ref == FALSE) %>% 
  ggplot(aes(abs(len_rel), fill=allele_type, color=allele_type )) + #fill=allele_type, color=allele_type)) 
  facet_grid(rows=vars(n), scales = 'free_y') +
  geom_freqpoly(position='identity', center=0, binwidth=0.1) +
  scale_fill_brewer("Type of allele", palette = "BrBG") +
  scale_color_brewer("Type of allele", palette = "BrBG") +
  geom_vline(xintercept = 50+1, linetype='dotted') +
  xlab("length difference of non-reference allele vs backbone") +  
  ylab("number of alleles") +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  theme(legend.position = "none")
```

I think this is interesting. It suggests that there is a lot of long insertions at low frequency. I wonder how many of them are attributable to the \textit{Rhamphochromis}? 
```{r}
# mean allele length deviation
df_allele %>% filter(n_paths_bio >= 2, allele_type != 'substitution/inversion', ref == FALSE) %>% 
  group_by(n, allele_type) %>% 
  summarise(avg_len_diff = mean(abs(len_rel))) %>% 
  pivot_wider(id_cols = n, names_from = allele_type, values_from = avg_len_diff) %>% 
  round(digits = 3) %>% 
  reactable()

# median allele length deviation
df_allele %>% filter(n_paths_bio >= 2, allele_type != 'substitution/inversion', ref == FALSE) %>% 
  group_by(n, allele_type) %>% 
  summarise(avg_len_diff = median(abs(len_rel))) %>% 
  pivot_wider(id_cols = n, names_from = allele_type, values_from = avg_len_diff) %>% 
  round(digits = 3) %>% 
  reactable()
```

## Allele sharing between samples
```{r}
# biallelic sites with complete allelic information, same as above
tmp <- df_allele %>% filter(ref == TRUE, n == 1, n_paths_bio == 2, n_uncalled == 0)
glue("{nrow(tmp)} bubbles with reference alleles adding up to {sum(tmp$len)} bp")

# biallelic sites, including places where the allele information is incomplete
tmp <- df_allele %>% filter(ref == TRUE, n == 1, n_paths_bio == 2)  
glue("{nrow(tmp)} bubbles with reference alleles adding up to {sum(tmp$len)} bp")

# multiallelic sites with complete allelic information
tmp <- df_allele %>% filter(ref == TRUE, n == 1, n_paths_bio > 2, n_uncalled == 0)
glue("{nrow(tmp)} bubbles with reference alleles adding up to {sum(tmp$len)} bp")

# multiallelic sites, including places where the allele information is incomplete
tmp <- df_allele %>% filter(ref == TRUE, n == 1, n_paths_bio > 2)  
glue("{nrow(tmp)} bubbles with reference alleles adding up to {sum(tmp$len)} bp")

# everything with complete allelic information (MOST COMPARABLE TO ABOVE)
tmp <- df_allele %>% filter(ref == TRUE, n == 1, n_uncalled == 0)
glue("{nrow(tmp)} bubbles with reference alleles adding up to {sum(tmp$len)} bp")

# everything
tmp <- df_allele %>% filter(ref == TRUE, n == 1)
glue("{nrow(tmp)} bubbles with reference alleles adding up to {sum(tmp$len)} bp")
```
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 sharing calculation for just the non-reference nodes (separate notebook).

## A few more plots (most are archived)
A histogram of the allele lengths is shown below.
```{r allele_len_histogram, fig.height=3}
df_allele %>% filter(n_paths_bio >= 2, allele_type != 'substitution/inversion', ref == FALSE) %>% 
  ggplot(aes(len+1, fill=allele_type, color=allele_type)) +
  scale_fill_brewer("Type of allele", palette = "BrBG") +
  scale_color_brewer("Type of allele", palette = "BrBG") +
  geom_vline(xintercept = 50+1, linetype='dotted') +
  annotate(geom='text', x=70, y=Inf, label="path length = 50", hjust=1, angle=90) +
  geom_area(stat='bin', position='identity', colour=NA, alpha=0.5) +
  geom_line(stat='bin', position='identity', size=1.2, alpha=0.8) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  #ggtitle("Distribution of non-reference allele lengths") + 
  xlab("log(path length + 1)") +  ylab("frequency")
```

Version of the relative lengths of SVs with `ggridges`.
```{r allele_lenrel_ridges}
df_allele %>% filter(n_paths_bio >= 2, allele_type != 'substitution/inversion', ref == FALSE) %>% 
  ggplot(aes(x=abs(len_rel), y=allele_type, fill=allele_type, color=allele_type)) +
  geom_density_ridges2(stat='binline', panel_scaling = FALSE) +
  geom_vline(xintercept = 50+1, linetype='dotted') +
  annotate(geom='text', x=35, y=0.7, label="relative length = 50", hjust=1, angle=0) +
  scale_fill_brewer("Type of allele", palette = "BrBG") +
  scale_color_brewer("Type of allele", palette = "BrBG") +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  #ggtitle("Length difference of non-reference alleles vs backbone") + 
  xlab("log(path length + 1)") + ylab(NULL) +
  theme(legend.position = "none")
```

It is also possible to calculate a ratio of the non-reference allele length against the reference. This gives an idea of the scale of difference of the non-reference alleles compared to the `calliptera`. 
```{r allele_lenratio_histogram, fig.height=3}
df_allele %>% filter(n_paths_bio >= 2, allele_type != 'substitution/inversion', ref == FALSE) %>% 
  mutate(ratio_allele_length = (len+1)/(len_ref+1) ) %>% 
  ggplot(aes(ratio_allele_length, fill=allele_type, color=allele_type)) +
  scale_fill_brewer("Type of allele", palette = "BrBG") +
  scale_color_brewer("Type of allele", palette = "BrBG") +
  geom_vline(xintercept = 50, linetype='dotted') +
  annotate(geom='text', x=75, y=14000, label="50/1", hjust=0) +
  geom_vline(xintercept = 1/50, linetype='dotted') +
  annotate(geom='text', x=1/75, y=14000, label="1/50", hjust=1) +
  geom_histogram(position='identity', alpha=0.5, center=0) +
  geom_line(stat='bin', position='identity', size=1) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  #ggtitle("Ratio of allele lengths: non-reference variant vs backbone") + 
  xlab("ratio of path lengths (variant / reference)") +  ylab("frequency")
```

```{r}
df_allele %>% filter(n_paths_bio >= 2, allele_type != 'substitution/inversion', ref == FALSE) %>% 
  mutate(ratio_allele_length = (len+1)/(len_ref+1) ) %>% 
  ungroup() %>% group_by(allele_type) %>% 
   summarise(mean_length = median(len), mean_ratio = median(ratio_allele_length), mean_ratio_inverse = 1/mean_ratio)
```

A few things I noticed:

- Note that there is a pseudocount of 1 on both the denominator and numerator.
- Complete deletions and insertions are the easiest to interpret, because the reference allele length is 0 + 1, so the position on the x-axis is pretty much the length of the longer allele (reference for deletions, non-reference for insertions).
- Complete deletions and insertions only start showing up at the 50bp mark, because this was the `min_variant_size` parameter in `minigraph`. There is a good number of variants measuring from there until 10,000 bp.
- The one thing that confuses me is that there is more complete insertions than complete deletions around the 1000 bp mark?? Is the biological or technical? ><

Speculating some transposon explanations:

- Complete insertions and deletions should be caused by newly inserted TEs. The fact that many of them are singletons might suggest these are **younger** insertions.
- Smaller partial insertions and deletions might reflect neutral mutations (at possibly older TEs or normal sequences).
- Is there a length difference between TE-rich and TE-poor SVs?

## Saving some data frames
```{r}
write_csv(df_allele, file = paste0("./local/malawi_haplochromines-alleles.csv"))
```