---
title: "Analysis of the minigraph parameter sweep"
output: html_document
---

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

In this document, we are analysing the results of sweeping through various `minigraph` parameters. We want to ensure that the built graphs are relatively robust across various metrics. The version of minigraph used is `v0.18`.

## Setup
```{r, cache=FALSE}
library(dplyr)
library(ggplot2)
library(tibble)
library(readr)
library(scales)
library(glue)
library(reactable)

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")
percent_around_default <- 0.05   # error bars canonical value
```

## Species order / backbone
Out of the 5 species, there are a total of 5! = 120 possible in which the genomes can be added to build the graph. With the default ordering being *astCal, mayZeb, troMau, aulStu, otoArg, copChr (utaka), rhaChi, rhaChi2* (*cmtaourR*), we tried running minigraph with all other possible orderings for comparison. We plot various metrics about the graph topology to see how robust they are.

First we need to preprocess the data.
```{r}
dfA1 <- read_csv(paste0(dir_sweep, "sweep_species_order.summary"), col_types='ciiiiiiiiidddddddd') %>% 
  mutate(species_order = gsub('S', 'R', species_order))
dfA2 <- read_csv(paste0(dir_sweep, "sweep_species_order.benchmark"), col_types='ciiddd') %>% 
  mutate(species_order = gsub('S', 'R', species_order))

df_species_order <- inner_join(dfA1, dfA2, by = c("species_order", "min_variant_size")) %>% 
  select(-bb_size) %>% 
  rename(mean_segment_length = mean, peak_memory = max_pss) %>% 
  mutate(peak_memory_gb = 1.04858 * peak_memory / 1000, 
         wallclock_time_hr = wallclock_time / 3600) %>% 
  select(species_order, n_segments:mean_segment_length, wallclock_time_hr, peak_memory_gb) %>% 
  tidyr::pivot_longer(-species_order, names_to='key') %>% 
  mutate(backbone = substr(species_order, 1, 1) %>% 
           forcats::fct_relevel(., 'c', 'm', 't', 'a', 'o', 'u', 'r', 'R')) %>% 
  mutate(
    key = forcats::fct_relevel(key,
    'n_segments', 'n_edges', 'n_bubbles',
    'linear_count', 'variable_count', 'mean_segment_length',
    'linear_size', 'variable_size', 'graph_size',
    'peak_memory_gb', 'wallclock_time_hr')) 

df_species_order_default <- df_species_order %>% 
  filter(species_order == 'cmtaourR') %>% 
  mutate(value_upper = value*(1 + percent_around_default), value_lower = value*(1 - percent_around_default))
```

We see that for most metrics about the graph (excluding benchmark parameters), the values fall within plus minus `r percent_around_default * 100`% of their respective value for the default species ordering (gray lines). This indicates that the graph generation is relatively consistent regardless of backbone, probably because the Lake Malawi cichlids are quite similar genetically to begin with.

```{r speciesorder_graphmetrics, fig.width = 9, fig.height = 5.5}
df_subset <- df_species_order %>% filter(!key %in% c('peak_memory_gb', 'wallclock_time_hr'))
df_species_order_default <- df_species_order_default %>% filter(!key %in% c('peak_memory_gb', 'wallclock_time_hr'))

ggplot(df_subset, aes(x = backbone, y = value, fill = backbone)) +
  facet_wrap(vars(key), ncol = 3, scales = 'free') +
  geom_hline(aes(yintercept = value_upper), df_species_order_default, color = 'gray70', linetype="dotted") +
  geom_hline(aes(yintercept = value), df_species_order_default, color = 'gray70') +
  geom_hline(aes(yintercept = value_lower), df_species_order_default, color = 'gray70', linetype="dotted") +
  geom_boxplot() +
  ylab(NULL) +
  scale_fill_manual(
    values = color_schema_species$color,
    breaks = color_schema_species$single_letter,
    labels = color_schema_species$species_abrv
  ) +
  scale_x_discrete(breaks = color_schema_species$single_letter, labels = color_schema_species$species_initials)
```

```{r speciesorder_graphmetrics-reordered, fig.width = 9, fig.height = 5.5}
color_schema_species_reordered <- color_schema_species %>% 
  mutate(single_letter = forcats::fct_relevel(single_letter, 'c', 'm', 't', 'a', 'r', 'o', 'u', 'R')) %>% 
  arrange(single_letter)

df_subset %>% 
  mutate(backbone = forcats::fct_relevel(backbone, 'c', 'm', 't', 'a', 'r', 'o', 'u', 'R') ) %>% 
  ggplot(aes(x = backbone, y = value, fill = backbone)) +
  facet_wrap(vars(key), ncol = 3, scales = 'free') +
  geom_hline(aes(yintercept = value_upper), df_species_order_default, color = 'gray70', linetype="dotted") +
  geom_hline(aes(yintercept = value), df_species_order_default, color = 'gray70') +
  geom_hline(aes(yintercept = value_lower), df_species_order_default, color = 'gray70', linetype="dotted") +
  geom_boxplot() +
  ylab(NULL) +
  scale_fill_manual(
    values = color_schema_species_reordered$color,
    breaks = color_schema_species_reordered$single_letter,
    labels = color_schema_species_reordered$species_abrv
  ) +
  scale_x_discrete(breaks = color_schema_species_reordered$single_letter, labels = color_schema_species_reordered$species_initials)
```

Wallclock runtime and memory benchmarks.
```{r speciesorder_benchmarks, fig.width = 6, fig.height = 2.5}
df_subset <- df_species_order %>% filter(key %in% c('peak_memory_gb', 'wallclock_time_hr'))
df_subset_default <- df_species_order_default %>% filter(key %in% c('peak_memory_gb', 'wallclock_time_hr'))

ggplot(df_subset, aes(x = backbone, y = value, fill = backbone)) +
  facet_wrap(vars(key), ncol = 3, scales = 'free') +
  geom_hline(aes(yintercept = value_upper), df_subset_default, color = 'gray70', linetype="dotted") +
  geom_hline(aes(yintercept = value), df_subset_default, color = 'gray70') +
  geom_hline(aes(yintercept = value_lower), df_subset_default, color = 'gray70', linetype="dotted") +
  geom_boxplot() + 
  ylab(NULL) +
  scale_fill_manual(
    values = color_schema_species$color,
    breaks = color_schema_species$single_letter,
    labels = color_schema_species$species_abrv
  ) +
  scale_x_discrete(
    breaks = color_schema_species$single_letter, labels = color_schema_species$species_initials
  )
```

## Minimum variant size
Now's let's fix the order to be the canonical one. We then vary the `-L` parameter, which constraints the minimum variant length allowed in the minigraph genome graph.

```{r minvariantsize_graphmetrics}
dfB1 <- read_csv(paste0(dir_sweep, "sweep_variant_size.summary"), col_types='ciiiiiiiiidddddddd')

df_vs_metrics <- dfB1 %>% 
  select(-bb_size) %>% 
  rename(mean_segment_length = mean) %>% 
  select(min_variant_size, n_segments:mean_segment_length) %>% 
  tidyr::pivot_longer(-min_variant_size, names_to='key') %>%
  mutate(key = forcats::fct_relevel(key,
    'n_segments', 'n_edges', 'n_bubbles',
    'linear_count', 'variable_count', 'mean_segment_length',
    'linear_size', 'variable_size', 'graph_size'))

df_vs_metrics_50 <- df_vs_metrics %>% 
  filter(min_variant_size == 50) %>% 
  mutate(value_upper = value*(1 + percent_around_default), value_lower = value*(1 - percent_around_default),
         xmin=1, xmax=1000)
df_vs_metrics_50 <- rbind(
  mutate(df_vs_metrics_50, min_variant_size=1),
  mutate(df_vs_metrics_50, min_variant_size=1000)
)

df_vs_metrics %>% 
  ggplot(aes(x = min_variant_size, y = value)) +
  facet_wrap(vars(key), ncol = 3, scales = 'free') + 
  scale_x_continuous("minimum variant size (-L)", trans='log10') +
  geom_ribbon(aes(x=min_variant_size, ymin=value_lower, ymax=value_upper), 
              df_vs_metrics_50, fill = 'gray90') +
  geom_hline(aes(yintercept = value), df_vs_metrics_50, linetype="dotted") +
  geom_line(color='black') + geom_point() + ylab(NULL)

df_vs_metrics %>% filter(min_variant_size != 1) %>% 
  ggplot(aes(x = min_variant_size, y = value)) +
  facet_wrap(vars(key), ncol = 3, scales = 'free') + 
  scale_x_continuous("minimum variant size (-L)", trans='log10') +
  geom_ribbon(aes(x=min_variant_size, ymin=value_lower, ymax=value_upper), 
              df_vs_metrics_50, fill = 'gray90') +
  geom_hline(aes(yintercept = value), df_vs_metrics_50, linetype="dotted") +
  geom_line(color='black') + geom_point() + ylab(NULL)
```

As before, plus minus `r percent_around_default * 100`% intervals from the default `L = 50` parameter are indicated. We note that the topology of the graph is noticeably different with decreasing counts across segments, edges, bubbles, linear nodes and variable nodes. However, it would seem that the amount of **information** encapsulated within the genome graph remains unchanged. The total size of the graph is relatively constant, as are those of linear and variable counterparts.

However, when the `-L` value becomes higher than a certain point (above 100), the linear part of the graph start to expand in size and the variable shrinks, as smaller scale interspecies variation are excluded. This means that at `-L` values below this point, any differences we notice in the graph topology might be caused by alignments behaving capriciously. The author did note that `minigraph` is not optimised to deal with smaller scale variation, so this might be what is at play here.

```{r minvariantsize_segmentlength, fig.width = 2.8, fig.height = 3.5}
dfB1 %>% select(min_variant_size, q0:q100) %>% 
ggplot(aes(x = min_variant_size)) +
  geom_boxplot(aes(group = min_variant_size, ymin = q0, lower = q25, middle = q50, upper = q75, ymax = q100), 
               stat = "identity", fill = "gray90") +
  coord_flip() +
  scale_x_continuous('minimum variant size (-L)', trans='log10', breaks = dfB1$min_variant_size ) +
  scale_y_log10("segment length (bp)",
        breaks = trans_breaks("log10", function(x) 10^x),
        labels = trans_format("log10", math_format(10^.x)))
```

Benchmark performances indicate relatively constant memory use, but an asymptotically decreasing wallclock time with increasing `-L` sizes.

```{r minvariantsize_benchmark, fig.width=2.5, fig.height=3.5}
dfB2 <- read_csv(paste0(dir_sweep, "sweep_variant_size.benchmark"), col_types='ciiddd')

df_vs_bmk_mean <- dfB2 %>% 
  mutate(peak_memory_gb = 1.04858 * max_pss / 1000, wallclock_time_hr = wallclock_time / 3600) %>% 
  select(min_variant_size, wallclock_time_hr, peak_memory_gb) %>% 
  group_by(min_variant_size) %>% 
  summarise(peak_memory_gb = mean(peak_memory_gb), wallclock_time_hr=mean(wallclock_time_hr)) %>% 
  tidyr::pivot_longer(-min_variant_size, names_to='key')

df_vs_bmk <- dfB2 %>% 
  mutate(peak_memory_gb = 1.04858 * max_pss / 1000, wallclock_time_hr = wallclock_time / 3600) %>% 
  select(min_variant_size, wallclock_time_hr, peak_memory_gb) %>% 
  tidyr::pivot_longer(-min_variant_size, names_to='key')

ggplot() +
  facet_wrap(vars(key), ncol = 1, scales = 'free') +
  scale_x_continuous('minimum variant size (-L)', trans='log10') +
  ylab(NULL) +
  geom_line(data=df_vs_bmk_mean, aes(x = min_variant_size, y = value), color='gray') + 
  geom_point(data=df_vs_bmk, aes(x = min_variant_size, y = value)) 
```

## Final printout of statistics for the canonical graph
Graph metrics.
```{r}
canonical_stats <- df_vs_metrics %>% 
  filter(min_variant_size == 50) %>% 
  select(key, value) %>% 
  column_to_rownames(var="key")

options(scipen=10)
print(canonical_stats, scientific = F)
options(scipen=0)
```

Print out of the boxplot values of the mean segment length for canonical `-L` parameter.
```{r}
dfB1 %>% filter(min_variant_size == 50) %>% select(q0, q25, q50, q75, q100)
```

```{r}
glue(
  "Backbone genome size               : {calliptera_genome_size}",
  "\nBackbone genome size - linear genome : {calliptera_flexible_size} ",
  " (backbone's variable part: {calliptera_flexible_percent}%)",
  "\nNOTE: These numbers do not account for coverage.",
  "\n",
  "\nTotal extra information      : {total_extra_info}",
  "\nExtra information per genome :  {extra_info_per_genome}",
  calliptera_genome_size = 880428986,
  calliptera_flexible_size = calliptera_genome_size - canonical_stats['linear_size', ],
  calliptera_flexible_percent = round(calliptera_flexible_size / calliptera_genome_size * 100, 2),
  total_extra_info = canonical_stats['graph_size',] - calliptera_genome_size,
  extra_info_per_genome = total_extra_info / (num_of_genomes - 1)
)
```

Print out graph stats for the different values of L.
```{r}
df_vs_metrics %>% 
  tidyr::pivot_wider(id_cols=min_variant_size, names_from=key, values_from=value) %>% 
  mutate(
    linear_percent = linear_size / graph_size * 100,
    variable_percent = variable_size / graph_size * 100,
    bubble_density = n_bubbles / 880428986 * 1000,
    additional_bases = (graph_size - 880428986) / 880428986 * 100,
    flex_on_bb = (880428986 - linear_size) / 880428986 * 100
  ) %>% 
  relocate(
    min_variant_size, n_segments, n_edges, mean_segment_length, 
    graph_size, linear_percent, variable_percent, 
    additional_bases, flex_on_bb, n_bubbles, bubble_density
  ) %>% 
  reactable(columns = list(min_variant_size = colDef(sticky = "left")))
```

## Other minigraph versions

### minigraph without the -C option (no base alignment)
**NOTE: this section uses the incorrect formulae to calculate linear and variable sizes (called core and variable here). Has yet to be corrected!** It only counts core genome using core segments that are in bubbles, which underestimates the size.

The statistics from altering the minimum variant size are shown.
```{r minvariantsize_graphmetrics_noC}
dfC1 <- read_csv(paste0(dir_sweep_archive, "sweep_variant_size_noC.summary"), col_types='ciiiiiiiiidddddddd')

df_vs_metrics <- dfC1 %>% 
  rename(mean_segment_length = mean) %>% 
  select(min_variant_size, n_segments:mean_segment_length) %>% 
  tidyr::pivot_longer(-min_variant_size, names_to='key') %>%
  mutate(key = forcats::fct_relevel(key,
    'n_segments', 'n_edges', 'n_bubbles',
    'core_count', 'flexible_count', 'mean_segment_length',
    'core_size', 'flexible_size', 'total_pangenome_size'))

df_vs_metrics_50 <- df_vs_metrics %>% 
  filter(min_variant_size == 50) %>% 
  mutate(value_upper = value*(1 + percent_around_default), value_lower = value*(1 - percent_around_default),
         xmin=1, xmax=1000)
df_vs_metrics_50 <- rbind(
  mutate(df_vs_metrics_50, min_variant_size=1),
  mutate(df_vs_metrics_50, min_variant_size=1000)
)

df_vs_metrics %>% 
  ggplot(aes(x = min_variant_size, y = value)) +
  facet_wrap(vars(key), ncol = 3, scales = 'free') + 
  scale_x_continuous("minimum variant size (-L)", trans='log10') +
  geom_ribbon(aes(x=min_variant_size, ymin=value_lower, ymax=value_upper), 
              df_vs_metrics_50, fill = 'gray90') +
  geom_hline(aes(yintercept = value), df_vs_metrics_50, linetype="dotted") +
  geom_line(color='black') + geom_point() + ylab(NULL)
```

### minigraph version 0.15
**NOTE: this section uses the incorrect formulae to calculate core and flexible sizes. Has yet to be corrected!** It only counts core genome using core segments that are in bubbles, which underestimates the size.

An older version of minigraph I initially used. It seems like the graph topology is a bit less robust depending on the specified species ordering.

```{r speciesorder_graphmetrics_0.15, fig.width = 9, fig.height = 6}
dfD1 <- read_csv(paste0(dir_sweep_0.15, "sweep_species_order.summary"), col_types='ciiiiiiiiidddddddd')
dfD2 <- read_csv(paste0(dir_sweep_0.15, "sweep_species_order.benchmark"), col_types='ciiddd')

df_species_order <- inner_join(dfD1, dfD2, by = c("species_order", "min_variant_size")) %>% 
  rename(mean_segment_length = mean, peak_memory = max_pss) %>% 
  mutate(peak_memory_gb = 1.04858 * peak_memory / 1000, 
         wallclock_time_hr = wallclock_time / 3600) %>% 
  select(species_order, n_segments:mean_segment_length, wallclock_time_hr, peak_memory_gb) %>% 
  tidyr::pivot_longer(-species_order, names_to='key') %>% 
  mutate(backbone = substr(species_order, 1, 1) %>% 
           forcats::fct_relevel(., 'c', 'm', 't', 'a', 'o', 'u', 'r', 'R')) %>% 
  mutate(
    key = forcats::fct_relevel(key,
    'n_segments', 'n_edges', 'n_bubbles',
    'core_count', 'flexible_count', 'mean_segment_length',
    'core_size', 'flexible_size', 'total_pangenome_size',
    'peak_memory_gb', 'wallclock_time_hr')) 

df_species_order_default <- df_species_order %>% 
  filter(species_order == 'cmtaourR') %>% 
  mutate(value_upper = value*(1 + percent_around_default), value_lower = value*(1 - percent_around_default))

df_subset <- df_species_order %>% filter(!key %in% c('peak_memory_gb', 'wallclock_time_hr'))
df_species_order_default <- df_species_order_default %>% filter(!key %in% c('peak_memory_gb', 'wallclock_time_hr'))

ggplot(df_subset, aes(x = backbone, y = value, fill = backbone)) +
  facet_wrap(vars(key), ncol = 3, scales = 'free') +
  geom_hline(aes(yintercept = value_upper), df_species_order_default, color = 'gray70', linetype="dotted") +
  geom_hline(aes(yintercept = value), df_species_order_default, color = 'gray70') +
  geom_hline(aes(yintercept = value_lower), df_species_order_default, color = 'gray70', linetype="dotted") +
  geom_boxplot() +
  ylab(NULL) +
  scale_fill_manual(
    values = color_schema_species$color,
    breaks = color_schema_species$single_letter,
    labels = color_schema_species$species_abrv
  ) +
  scale_x_discrete(breaks = color_schema_species$single_letter, labels = color_schema_species$species_initials)
```

The results of sweeping the minimum variant size,  however, is closer to the latest version with base alignment than without.
```{r minvariantsize_graphmetrics_0.15}
dfE1 <- read_csv(paste0(dir_sweep_0.15, "sweep_variant_size.summary"), col_types='ciiiiiiiiidddddddd')

df_vs_metrics <- dfE1 %>% 
  rename(mean_segment_length = mean) %>% 
  select(min_variant_size, n_segments:mean_segment_length) %>% 
  tidyr::pivot_longer(-min_variant_size, names_to='key') %>%
  mutate(key = forcats::fct_relevel(key,
    'n_segments', 'n_edges', 'n_bubbles',
    'core_count', 'flexible_count', 'mean_segment_length',
    'core_size', 'flexible_size', 'total_pangenome_size'))

df_vs_metrics_50 <- df_vs_metrics %>% 
  filter(min_variant_size == 50) %>% 
  mutate(value_upper = value*(1 + percent_around_default), value_lower = value*(1 - percent_around_default),
         xmin=1, xmax=1000)
df_vs_metrics_50 <- rbind(
  mutate(df_vs_metrics_50, min_variant_size=1),
  mutate(df_vs_metrics_50, min_variant_size=1000)
)

df_vs_metrics %>% 
  ggplot(aes(x = min_variant_size, y = value)) +
  facet_wrap(vars(key), ncol = 3, scales = 'free') + 
  scale_x_continuous("minimum variant size (-L)", trans='log10') +
  geom_ribbon(aes(x=min_variant_size, ymin=value_lower, ymax=value_upper), 
              df_vs_metrics_50, fill = 'gray90') +
  geom_hline(aes(yintercept = value), df_vs_metrics_50, linetype="dotted") +
  geom_line(color='black') + geom_point() + ylab(NULL)
```