---
title: "Pangenome growth and wider context"
output: html_document
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  cache = TRUE, cache.path = '.cache_knitr/pangenome_growth/',
  fig.width = 7, fig.height = 4, dpi = 350,
  fig.path='figure/5-pangenome_growth/',
  dev = c('png', 'pdf')
)
```

In this document, we make some plots to illustrate changes in graph size and pangenome for a wider context of species.

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

setwd(paste0(Sys.getenv("HOME"), "/code/malawi_transposon/script/R_SCRIPT/"))

# color schema for species
color_schema_species <- readr::read_csv(
  paste0(Sys.getenv("HOME"), 
         "/code/malawi_transposon/metadata/species_colorscheme.csv"), comment='#') %>% 
  mutate(color = paste0('#', color))

# most useful stats
useful_stats <- unlist(strsplit(glue(
  "n_segments,n_edges,graph_size,linear_percent,variable_percent,",
  "bb_extra_perc,bb_cov_perc,n_bubbles,bubble_density,mean_segment_length"), split = ','))
```

## Data import
Import final graph stats and individual sample contributions + coverage by sample.
```{r}
df_stats <- readr::read_csv(
  paste0(Sys.getenv("HOME"),
    "/code/malawi_transposon/cloud/workflow/wider_context_comparisons/results/",
    "allgraphs.stat.finalvalue.summary") )

df_stats_bysample <- readr::read_csv(
  paste0(Sys.getenv("HOME"),
    "/code/malawi_transposon/cloud/workflow/wider_context_comparisons/results/",
    "allgraphs.stat.bysample.finalvalue.summary") ) %>% 
  mutate(
    species = if_else(species == 'rhaChiTwo', 'rhaChi2', species), 
    rank = rank + 1
  )

df_cov_bysample <- readr::read_csv(
  paste0(Sys.getenv("HOME"),
    "/code/malawi_transposon/cloud/workflow/wider_context_comparisons/results/",
    "allgraphs.cov.bysample.finalvalue.summary") ) %>% 
  mutate(
    species = if_else(species == 'rhaChiTwo', 'rhaChi2', species), 
    rank = rank + 1
  )
```

Import dataframes related to pangenome stat and coverage growth.
```{r}
tmp_df1 <- readr::read_csv(
  paste0(Sys.getenv("HOME"),
    "/code/malawi_transposon/cloud/workflow/wider_context_comparisons/results/",
    "allgraphs.stat.growth.summary"))

tmp_df2 <- readr::read_csv(
  paste0(Sys.getenv("HOME"), 
    "/code/malawi_transposon/cloud/workflow/wider_context_comparisons/results/",
    "allgraphs.covstat.growth.summary")) 

df_growth_pangenome_raw <- inner_join(tmp_df1, tmp_df2, by=c('graph', 'sample')) %>% 
  mutate(
    sample = sample + 1, # correct zero-indexed rank
    bubble_density = n_bubbles / bb_size * 1000,
    bubble_density_corr = n_bubbles / bb_cov * 1000,
    linear_percent = linear_size / graph_size * 100,
    variable_percent = variable_size / graph_size * 100,
    variable_size_bb_rel = variable_size_bb / bb_size * 100,
    variable_size_nonref_rel = variable_size_nonref / bb_size * 100,
    variable_size_bb_rel_corr = variable_size_bb / bb_ubiq * 100,
    variable_size_nonref_rel_corr = variable_size_nonref / bb_ubiq * 100
  )

rm(tmp_df1, tmp_df2)
```

Add some default values to some features when only the backbone is present (i.e. when rank = 1).
```{r}
tmp_df <- df_growth_pangenome_raw %>% 
  filter(sample == 2) %>% 
  mutate(
    sample = 1,
    n_segments = 0,
    n_edges = 0,
    n_bubbles = 0, 
    mean_segment_length = NA,
    graph_size = bb_size,
    linear_size = bb_size,
    linear_count = 0,
    variable_size = 0,
    variable_count = 0,
    variable_size_bb = 0,
    variable_count_bb = 0,
    variable_size_nonref = 0,
    variable_count_nonref = 0,
    bb_size = bb_size,
    bb_extra_perc = 0,
    bb_cov = NULL,  # 0
    bb_cov_perc = NULL, # 0
    bb_extra_perc_corr = 0,
    bb_ubiq = NULL, # bb_size
    bb_ubiq_perc = NULL, # 100
    bb_conserv_perc = 100,
    core_size = bb_size,
    flex_size = 0,
    pangenome_size = bb_size,
    core_perc = 100,
    flex_perc = 0,
    bubble_density = 0,
    bubble_density_corr = 0,
    linear_percent = 100,
    variable_percent = 0,
    variable_size_bb_rel = 0,
    variable_size_nonref_rel = 0,
    variable_size_bb_rel_corr = 0,
    variable_size_nonref_rel_corr = 0
  )

df_growth_pangenome_wide <- df_growth_pangenome_raw %>% 
  bind_rows(tmp_df) %>% 
  arrange(graph, sample) %>% 
  mutate(graph = forcats::fct_relevel(graph, "malawi", "african_cichlid", "great_ape", "chicken"))
rm(tmp_df)

df_growth_pangenome_wide
```

Convert to long format needed for some of the plots.
```{r}
df_growth_pangenome <- df_growth_pangenome_wide %>% gather(key, value, -graph, -sample)
```

Extract a subset of the dataframes focused just on Lake Malawi graphs, but with different species as the backbone. Do some preprocessing to extract the backbone from the graph name.
```{r}
extract_backbone_name <- function(input_df) {
  output_df <- input_df %>% 
    filter(stringr::str_detect(graph, "malawi")) %>% 
    mutate(backbone = str_replace(graph, "malawi_", "")) %>% 
    mutate(backbone = if_else(backbone=="malawi", "astCal", backbone)) %>% 
    mutate(backbone = if_else(backbone=="rhaChiTwo", "rhaChi2", backbone)) %>% 
    mutate(backbone = factor(backbone, levels = color_schema_species$species_abrv))
  return(output_df)
}

df_growth_pangenome_malawi <- extract_backbone_name(df_growth_pangenome)
df_stats_bysample_malawi <- extract_backbone_name(df_stats_bysample)
df_cov_bysample_malawi <- extract_backbone_name(df_cov_bysample)
```


## Plots for the default Lake Malawi graph by itself
Individual stat growth plots for select empirical features: number of segments, number of edges, graph size, number of bubbles, bubble density, mean segment length
```{r malawi-solo_panel, fig.width=2.5, fig.height=2.2}
# number of segments
df_growth_pangenome %>% 
  filter(key %in% c("n_segments"), graph == "malawi") %>% 
  ggplot(aes(x=sample, y=value)) +
  geom_line() + geom_point() +
  xlab("rank of sample") +
  ylab("number of segments") +
  theme(aspect.ratio=1)

# number of edges
df_growth_pangenome %>% 
  filter(key %in% c("n_edges"), graph == "malawi") %>% 
  ggplot(aes(x=sample, y=value)) +
  geom_line() + geom_point() +
  xlab("rank of sample") +
  ylab("number of edges") +
  theme(aspect.ratio=1)

# graph size
df_growth_pangenome %>% 
  filter(key %in% c("graph_size"), graph == "malawi") %>% 
  ggplot(aes(x=sample, y=value/10**9)) +
  geom_line() + geom_point() +
  xlab("rank of sample") +
  ylab("total size of graph (Gbp)") +
  theme(aspect.ratio=1)

# number of bubbles
df_growth_pangenome %>% 
  filter(key %in% c("n_bubbles"), graph == "malawi") %>% 
  ggplot(aes(x=sample, y=value)) +
  geom_line() + geom_point() +
  xlab("rank of sample") +
  ylab("number of bubbles") +
  theme(aspect.ratio=1)

# bubble density
df_growth_pangenome %>% 
  filter(key %in% c("bubble_density"), graph == "malawi") %>% 
  ggplot(aes(x=sample, y=value)) +
  geom_line() + geom_point() +
  xlab("rank of sample") +
  ylab("bubbles detected per kbp") +
  theme(aspect.ratio=1)

# cumulative coverage
df_growth_pangenome %>% 
  filter(key %in% c("bb_cov_perc"), graph == "malawi") %>% 
  ggplot(aes(x=sample, y=value)) +
  geom_line() + geom_point() +
  xlab("rank of sample") +
  ylab("backbone coverage (%)") +
  theme(aspect.ratio=1)

# mean segment length
df_growth_pangenome %>% 
  filter(key %in% c("mean_segment_length"), graph == "malawi") %>% 
  ggplot(aes(x=sample, y=value)) +
  geom_line() + geom_point() +
  xlab("rank of sample") +
  ylab("mean segment length") +
  theme(aspect.ratio=1)
```

```{r malawi-raw_prop, fig.width=3.5, fig.height=2.5}
# linear and variable segments
df_growth_pangenome %>% 
  filter(key %in% c("linear_percent", "variable_percent"), graph == "malawi") %>% 
  ggplot(aes(x=sample, y=value, linetype=key, shape=key)) +
  geom_line() + geom_point() +
  coord_cartesian(ylim=c(0,100)) +
  scale_linetype_manual(values = c(2,3), labels=c('linear', 'variable')) +
  scale_shape_manual(values = c(19,15), labels=c('linear', 'variable')) +
  xlab("rank of sample") +
  ylab("percent sequence of graph") +
  theme(legend.position = 'top', legend.title=element_blank())
```

```{r malawi-pangenome, fig.width=3.5, fig.height=2.5}
# core and flexible (NOT USED: TOO CONFUSING)
df_growth_pangenome %>% 
  filter(key %in% c("core_perc", "flex_perc"), graph == "malawi") %>% 
  ggplot(aes(x=sample, y=value, linetype=key, shape=key)) +
  geom_line() + geom_point() +
  coord_cartesian(ylim=c(0,100)) +
  scale_linetype_manual(values = c(2,3), labels=c('core', 'flexible')) +
  scale_shape_manual(values = c(19,15), labels=c('core', 'flexible')) +
  xlab("rank of sample") +
  ylab("estimated percent of pangenome") +
  theme(legend.position = 'top', legend.title=element_blank())
```

```{r malawi-coverage, fig.width=3.3, fig.height=2.8}
# coverage (NOT USED: TOO CONFUSING)
df_growth_pangenome %>% 
  filter(key %in% c("bb_cov_perc", "bb_ubiq_perc"), graph == "malawi") %>% 
  drop_na() %>% 
  ggplot(aes(x=sample, y=value, color=key, shape=key)) +
  geom_line() + geom_point() +
  coord_cartesian(ylim=c(75,100)) +
  scale_color_manual(values=c("#92D050", "#018571"), labels=c('at least one', 'all')) +
  scale_shape_manual(values = c(19,15), labels=c('at least one', 'all')) +
  xlab("rank of sample") +
  ylab("spanning coverage (%)") +
  theme(legend.position = 'top', legend.title=element_blank())
```

**SUGGESTED FROM RICHARD DISCUSSIONS:** show how the nonreference variable sequences increase indefinitely, whereas the backbone variability starts to plateau.
```{r malawi-variable_v_bb, fig.width=4, fig.height=3.2}
stats_to_plot <- c("variable_size", "variable_size_nonref", "variable_size_bb")
df_growth_pangenome %>% 
  filter(graph == "malawi", key %in% stats_to_plot) %>% 
  ggplot(aes(x=sample, y=value/10**6, color=key)) +
  geom_point() +
  geom_line() + 
  xlab("rank of sample") +
  ylab("length of variable sequence (Mbp)") +
  scale_color_manual(
    NULL, breaks=stats_to_plot, 
    values=c('gray30', color_schema_species$color[3], color_schema_species$color[1]), 
    labels=c('all segments', 'nonreference', 'backbone')
  ) +
  scale_y_continuous(sec.axis = sec_axis(~ . * 10**6 / 880428986 * 100, "% of backbone length"), ) +
  theme(legend.position='top')
```

```{r malawi-variable_v_bb-poster, fig.width=4, fig.height=2.8}
stats_to_plot <- c("variable_size", "variable_size_nonref", "variable_size_bb")
df_growth_pangenome %>% 
  filter(graph == "malawi", key %in% stats_to_plot) %>% 
  ggplot(aes(x=sample, y=value/10**6, color=key)) +
  geom_point() +
  geom_line() + 
  xlab("Number of assemblies") +
  ylab("Sequence detected\ninside graph bubbles (Mbp)") +
  scale_color_manual(
    NULL, breaks=stats_to_plot, 
    values=c('gray30', color_schema_species$color[3], color_schema_species$color[1]), 
    labels=c('all', 'non-backbone', 'backbone')
  ) +
  scale_y_continuous(sec.axis = sec_axis(~ . * 10**6 / 880428986 * 100, "% of backbone length"), ) +
  theme_linedraw() +
  theme(legend.position='top')
```

Quantifying the numbers a bit (for paper):
```{r}
tmp_df <- df_stats_bysample %>% 
  filter(graph == 'malawi', species != 'astCal') %>% 
  select(species, n_bases) %>% 
  mutate(n_bases_perc = n_bases / 880428986 * 100)

tmp_df

glue(
  "PERCENTAGE mean: ", round(mean(tmp_df$n_bases_perc), 3), ", median: ", round(median(tmp_df$n_bases_perc), 3),
  "\nN_BASE mean: ", round(mean(tmp_df$n_bases), 3), ", median: ", round(median(tmp_df$n_bases), 3)
)
```

```{r wider-variable_v_bb, fig.width=5.5, fig.height=3.5}
stats_to_plot <- c("variable_size_nonref_rel", "variable_size_bb_rel")
df_growth_pangenome_malawi %>% 
  filter(key %in% stats_to_plot) %>% 
  ggplot(aes(x=sample, y=value, color=backbone, shape=backbone)) +
  geom_line(alpha=0.8, linetype=2) +
  geom_point() +
  facet_grid(.~factor(key, levels=stats_to_plot, labels=c("nonreference", "backbone"))) +
  scale_color_manual(values=color_schema_species$color, breaks=color_schema_species$species_abrv) +
  scale_shape_manual(values = c(19,19,19,17,17,17,15,15)) +
  xlab("rank of sample") +
  ylab("variable length as a\n% of backbone length") +
  theme(legend.position = 'top')

df_growth_pangenome %>% 
  filter(!stringr::str_detect(graph, "^malawi_")) %>% 
  filter(key %in% stats_to_plot) %>% 
  ggplot(aes(x=sample, y=value, color=graph, shape=graph)) +
  geom_line(alpha=0.8, linetype=2) +
  geom_point() +
  facet_grid(.~factor(key, levels=stats_to_plot, labels=c("nonreference", "backbone"))) +
  scale_shape_manual(values = c(19,19,17,15)) +
  xlab("rank of sample") +
  ylab("variable length as a\n% of backbone length") +
  theme(legend.position = 'top')
```

```{r wider-variable_v_bb_corr, fig.width=5.5, fig.height=3.2}
stats_to_plot <- c("variable_size_nonref_rel_corr", "variable_size_bb_rel_corr")
df_growth_pangenome_malawi %>% 
  filter(key %in% stats_to_plot) %>% 
  ggplot(aes(x=sample, y=value, color=backbone, shape=backbone)) +
  geom_line(alpha=0.8, linetype=2) +
  geom_point() +
  facet_grid(.~factor(key, levels=stats_to_plot, labels=c("nonreference", "backbone"))) +
  scale_color_manual(values=color_schema_species$color, breaks=color_schema_species$species_abrv) +
  scale_shape_manual(values = c(19,19,19,17,17,17,15,15)) +
  xlab("rank of sample") +
  ylab("variable length as a\n% of covered backbone") +
  theme(legend.position = 'top')

df_growth_pangenome %>% 
  filter(!stringr::str_detect(graph, "^malawi_")) %>% 
  filter(key %in% stats_to_plot) %>% 
  ggplot(aes(x=sample, y=value, color=graph, shape=graph)) +
  geom_line(alpha=0.8, linetype=2) +
  geom_point() +
  facet_grid(.~factor(key, levels=stats_to_plot, labels=c("nonreference", "backbone"))) +
  scale_shape_manual(values = c(19,19,17,15)) +
  xlab("rank of sample") +
  ylab("variable length as a\n% of covered backbone") +
  theme(legend.position = 'top')
```


## Stat growth plots across diverse graphs + different backbones for Lake Malawi
We make these growth plots across the diverse graphs, for a selected number of features at a time (for presentation purposes).

Linear and variable percent in graph.
```{r wider-raw_prop, fig.width=3.8, fig.height=3.5}
df_growth_pangenome_malawi %>% 
  filter(key %in% c("linear_percent", "variable_percent")) %>% 
  ggplot(aes(x=sample, y=value, linetype=key, color=backbone, shape=backbone)) +
  geom_line() + geom_point() +
  coord_cartesian(ylim=c(0,100)) +
  scale_linetype_manual(values = c(2,3), labels=c('linear', 'variable')) +
  scale_color_manual(values=color_schema_species$color, breaks=color_schema_species$species_abrv) +
  scale_shape_manual(values = c(19,19,19,17,17,17,15,15)) +
  xlab("rank of sample") +
  ylab("percent sequence of graph") +
  theme(legend.position = 'right')

df_growth_pangenome %>% 
  filter(!stringr::str_detect(graph, "^malawi_")) %>% 
  filter(key %in% c("linear_percent", "variable_percent")) %>% 
  ggplot(aes(x=sample, y=value, linetype=key, color=graph, shape=graph)) +
  geom_line() + geom_point() +
  coord_cartesian(ylim=c(0,100)) +
  scale_linetype_manual(values = c(2,3), labels=c('linear', 'variable')) +
  scale_shape_manual(values = c(19,19,17,15)) +
  xlab("rank of sample") +
  ylab("percent sequence of graph") +
  theme(legend.position = 'right')
```

Core and flexible estimation (NOT USED, TOO CONFUSING).
```{r wider-pangenome, fig.width=3.8, fig.height=3.5}
# core and flexible (NOT USED: TOO CONFUSING)
df_growth_pangenome %>% 
  filter(!stringr::str_detect(graph, "^malawi_")) %>% 
  filter(key %in% c("core_perc", "flex_perc")) %>% 
  ggplot(aes(x=sample, y=value, linetype=key, color=graph, shape=graph)) +
  geom_line() + geom_point() +
  coord_cartesian(ylim=c(0,100)) +
  scale_linetype_manual(values = c(2,3), labels=c('core_percent', 'flex_percent')) +
  scale_shape_manual(values = c(19,19,17,15)) +
  xlab("rank of sample") +
  ylab("estimated percent of pangenome") +
  theme(legend.position = 'right')

df_growth_pangenome_malawi %>% 
  filter(key %in% c("core_perc", "flex_perc")) %>% 
  ggplot(aes(x=sample, y=value, linetype=key, color=backbone, shape=backbone)) +
  geom_line() + geom_point() +
  coord_cartesian(ylim=c(0,100)) +
  scale_linetype_manual(values = c(2,3), labels=c('core_percent', 'flex_percent')) +
  scale_color_manual(values=color_schema_species$color, breaks=color_schema_species$species_abrv) +
  scale_shape_manual(values = c(19,19,19,17,17,17,15,15)) +
  xlab("rank of sample") +
  ylab("estimated percent of pangenome") +
  theme(legend.position = 'right')
```

Panel of all the useful stats.
```{r wider-panel_useful_stat, fig.width=9.5, fig.height=4}
df_growth_pangenome_malawi %>% 
  filter(key %in% useful_stats) %>% 
  ggplot(aes(x=sample, y=value, color=backbone, linetype=backbone)) +
  geom_line() +
  xlab("rank of sample") + ylab(NULL) +
  scale_color_manual(values=color_schema_species$color, breaks=color_schema_species$species_abrv) +
  scale_linetype_manual(values = c(1,2,2,3,3,3,4,4), breaks=color_schema_species$species_abrv) +
  facet_wrap(.~factor(key, levels=useful_stats), scales="free_y", nrow=2) +
  theme(legend.position = 'top', aspect.ratio=4/5) +
  guides(color = guide_legend(nrow = 1))

df_growth_pangenome %>%
  filter(!stringr::str_detect(graph, "^malawi_")) %>% 
  filter(key %in% useful_stats ) %>% 
  ggplot(aes(x=sample, y=value, shape=graph, color=graph, linetype=graph)) +
  geom_line() +
  geom_point() +
  scale_color_discrete('graph') +
  scale_shape_manual(values = c(19,19,17,15)) +
  scale_linetype_manual(values = c(1,1,2,2)) +
  xlab("rank of sample") + ylab(NULL) +
  facet_wrap(.~factor(key, levels=useful_stats), scales = 'free_y', nrow=2) +
  theme(legend.position = 'top', aspect.ratio=4/5)
```

Select stats for presentation.
```{r wider-panel_cov_stat, fig.width=8, fig.height=4.4}
variables_to_plot <- paste0(
  "linear_size,variable_size,bb_extra_perc,bb_cov_perc,n_bubbles,bubble_density,bubble_density_corr,") %>% 
  strsplit(split=',') %>%  
  unlist()

df_growth_pangenome_malawi %>% 
  filter(key %in% variables_to_plot) %>% 
  ggplot(aes(x=sample, y=value, color=backbone, linetype=backbone)) +
  geom_line() +
  xlab("rank of sample") +
  scale_color_manual(values=color_schema_species$color, breaks=color_schema_species$species_abrv) +
  scale_linetype_manual(values = c(1,2,2,3,3,3,4,4), breaks=color_schema_species$species_abrv) +
  facet_wrap(.~factor(key, levels=variables_to_plot), scales="free_y", ncol=4) +
  theme(legend.position = 'top', aspect.ratio=4/5)

df_growth_pangenome %>% 
  filter(!stringr::str_detect(graph, "^malawi_")) %>% 
  filter(key %in% variables_to_plot) %>% 
  ggplot(aes(x=sample, y=value, shape=graph, color=graph, linetype=graph)) +
  geom_line() +
  geom_point() +
  xlab("rank of sample") +
  scale_color_discrete('graph') +
  scale_shape_manual(values = c(19,19,17,15)) +
  scale_linetype_manual(values = c(1,1,2,2)) +
  facet_wrap(.~factor(key, levels=variables_to_plot), scales="free_y", ncol=4) +
  theme(legend.position = 'top', aspect.ratio=4/5)

df_growth_pangenome %>% 
  filter(graph == "malawi", key %in% variables_to_plot) %>%
  ggplot(aes(x=sample, y=value, color='malawi')) +
  geom_line() +
  geom_point() +
  xlab("rank of sample") +
  facet_wrap(.~factor(key, levels=variables_to_plot), scales="free_y", ncol=4) +
  theme(legend.position = 'top', aspect.ratio=4/5)
```

## Contribution of individual samples to overall graph size
Some preprocessing and sorting the factor levels.
```{r}
df_stats_bysample_2 <- df_stats_bysample %>% 
  filter(!stringr::str_detect(graph, "malawi_")) %>% 
  select(-sample) %>% 
  rename(sample = rank) %>% 
  inner_join(df_stats %>% select(graph, total_size=graph_size, total_segments=n_segments), by = 'graph') %>% 
  mutate(graph = forcats::fct_relevel(graph, "malawi", "african_cichlid", "great_ape", "chicken"))

df_stats_bysample_malawi_2 <- df_stats_bysample_malawi %>% 
  select(-sample) %>% 
  rename(sample = rank) %>% 
  inner_join(df_stats %>% select(graph, total_size=graph_size, total_segments=n_segments), by = 'graph')
```

Lake Malawi alone.
```{r malawi-bar_contribution, fig.width=3.3, fig.height=2.8}
df_stats_bysample_2 %>% 
  filter(graph=="malawi") %>% 
  mutate(species = factor(species, levels=color_schema_species$species_abrv)) %>% 
  ggplot(aes(x=species, y=n_bases/10**6, fill=species)) +
  scale_x_discrete(NULL) +
  scale_y_continuous("nucleotides contributed (Mbp)", expand = expansion(mult = c(.05, .1)) ) +
  scale_fill_manual(
    values = color_schema_species$color,
    breaks = color_schema_species$species_abrv,
    labels = color_schema_species$species_abrv
  ) +
  geom_bar(stat='identity') +
  geom_text(aes(label = paste0(format(round(n_bases/total_size*100,2), n_small=2), "%") ), 
            vjust = -0.5, position = position_dodge(width = .9), size=rel(2.4)) +
  theme(legend.position = 'none', axis.text.x = element_text(angle = 30, hjust = 1))

df_stats_bysample_2 %>% 
  filter(graph=="malawi") %>% 
  mutate(species = factor(species, levels=color_schema_species$species_abrv)) %>% 
  ggplot(aes(x=species, y=n_segment/1000, fill=species)) +
  scale_x_discrete(NULL) +
  scale_y_continuous("segments contributed", expand = expansion(mult = c(.05, .1)) ) +
  scale_fill_manual(
    values = color_schema_species$color,
    breaks = color_schema_species$species_abrv,
    labels = color_schema_species$species_abrv
  ) +
  geom_bar(stat='identity') +
  geom_text(aes(label = paste0(format(round(n_segment/total_segments*100,2), n_small=2), "%") ), 
            vjust = -0.5, position = position_dodge(width = .9), size=rel(2.4)) +
  theme(legend.position = 'none', axis.text.x = element_text(angle = 30, hjust = 1))
```

```{r wider-bar_contribution_robustness, fig.width=6, fig.height=5}
df_stats_bysample_2 %>% 
  filter(!stringr::str_detect(graph, "^malawi_")) %>% 
  ggplot(aes(x=factor(sample), y=n_bases)) +
  xlab("rank of sample") +
  scale_y_continuous("nucleotides contributed to graph", expand = expansion(mult = c(.05, .25)) ) +
  geom_bar(stat='identity') +
  geom_text(aes(label = paste0(format(round(n_bases/total_size*100,1), n_small=1), "%") ), 
            vjust = -0.5, size=rel(2.5), position = position_dodge(width = .9)) +
  facet_wrap(.~graph, scales='free', nrow=2) +
  theme(legend.position = 'top') +
  theme(aspect.ratio = 4/5)
```

```{r malawi-bar_contribution_robustness, fig.width=11, fig.height=5}
df_stats_bysample_malawi_2 %>% 
  ggplot(aes(x=factor(sample), y=n_bases, fill=species)) +
  xlab("rank of sample") +
  scale_y_continuous("nucleotides contributed to graph", expand = expansion(mult = c(.05, .25)) ) +
  scale_fill_manual(
    values = color_schema_species$color,
    breaks = color_schema_species$species_abrv,
    labels = color_schema_species$species_abrv
  ) +
  geom_bar(stat='identity') +
  geom_text(aes(label = paste0(format(round(n_bases/total_size*100,1), n_small=1), "%") ), 
            vjust = -0.5, size=rel(2.5), position = position_dodge(width = .9)) +
  facet_wrap(.~backbone, scales='free', nrow=2) +
  theme(legend.position = 'top') +
  guides(fill = guide_legend(nrow = 1))
```

## Coverage of individual samples
Some preprocessing and sorting the factor levels.
```{r}
df_cov_bysample_2 <- df_cov_bysample %>% 
  select(-sample) %>% 
  rename(sample = rank) %>% 
  mutate(graph = forcats::fct_relevel(graph, "malawi", "african_cichlid", "great_ape", "chicken"))

df_cov_bysample_malawi_2 <- df_cov_bysample_malawi %>% 
  select(-sample) %>% 
  rename(sample = rank)
```

```{r malawi-bar_coverage_bysample, fig.width=3.3, fig.height=2.8}
df_cov_bysample_2 %>% 
  filter(graph=="malawi") %>% 
  mutate(species = factor(species, levels=color_schema_species$species_abrv)) %>% 
  ggplot(aes(x=species, y=bb_cov_perc_bysample, fill=species)) +
  geom_bar(stat='identity') +
  geom_text(aes(y=if_else(bb_cov_perc_bysample > 90, bb_cov_perc_bysample, 90),
                label = format(round(bb_cov_perc_bysample,2), n_small=2) ), 
            size=rel(2.8), vjust = -0.5, position = position_dodge(width = .9)) +
  geom_hline(yintercept=90, linetype=2) +
  scale_x_discrete(NULL) +
  scale_y_continuous("spanning coverage (%)", breaks=seq(0,100,20),
                     expand = expansion(mult = c(.05, .1)) ) +
  scale_fill_manual(
    values = color_schema_species$color,
    breaks = color_schema_species$species_abrv,
    labels = color_schema_species$species_abrv
  ) +
  theme(legend.position = 'none', axis.text.x = element_text(angle = 30, hjust = 1))
```

```{r malawi-bar_coverage_bysample_robustness, fig.width=9.2, fig.height=4.2}
df_cov_bysample_malawi_2 %>% 
  ggplot(aes(x=factor(sample), y=bb_cov_perc_bysample, fill=species)) +
  geom_bar(stat='identity') +
  geom_text(aes(y=if_else(bb_cov_perc_bysample > 90, bb_cov_perc_bysample, 90),
                label = paste0(format(round(bb_cov_perc_bysample,1), n_small=1)) ), 
            vjust = -0.8, size=2.8, position = position_dodge(width = .9)) +
  geom_hline(yintercept=90, linetype=2) +
  xlab("rank of sample") +
  scale_y_continuous("spanning coverage on backbone", breaks=seq(0,100,20),
                     expand = expansion(mult = c(.05, .25)) ) +
  scale_fill_manual(
    "added assembly",
    values = color_schema_species$color,
    breaks = color_schema_species$species_abrv,
    labels = color_schema_species$species_abrv
  ) +
  facet_wrap(.~backbone, scales='free_x', nrow=2) +
  theme(legend.position='top') +
  guides(fill = guide_legend(nrow = 1))
```

## Some more coverage stuff
```{r malawi-bar_misc, fig.width=9.2, fig.height=2.4}
variables_to_plot <- paste0(
  "bb_extra_perc_corr,bb_cov_perc,bb_ubiq_perc,bb_conserv_perc") %>% 
  strsplit(split=',') %>%  
  unlist()
labels_to_plot <- c(
  'extra bases vs backbone', 
  'cumulative coverage',
  'ubiquitous coverage',
  'detectable in all other samples'
  )

df_growth_pangenome_malawi %>% 
  filter(sample == 8, key %in% variables_to_plot) %>% 
  mutate(key = factor(key, levels=variables_to_plot, labels=labels_to_plot)) %>% 
  ggplot(aes(x=backbone, y=value, fill=backbone)) +
  geom_bar(stat='identity') +
  facet_wrap(.~key, nrow=1)  +
  scale_y_continuous("percentage sequence", expand = expansion(mult = c(.05, .15)) ) +
  scale_fill_manual(
    values = color_schema_species$color,
    breaks = color_schema_species$species_abrv,
    labels = color_schema_species$species_abrv
  ) +
  geom_text(aes(label = paste0(format(round(value,1), n_small=1)) ), 
            vjust = -0.5, size=2.8, position = position_dodge(width = .9)) +
  theme(legend.position = 'none', axis.text.x = element_text(angle = 30, hjust = 1))
```

## Printout of final statistics
Final stats for all the graphs can be retrieved for the input file for `df_stats`, except for these:
* `bb_ubiq_perc` refers to the part of the backbone where alleles can be detectable in all species (includes backbone flexible alleles).
* `bb_conserv_perc` is the proportion of linear segments in `bb_ubiq_perc`
* `bubble density` is per kbp and accounts for covered regions only
```{r}
tmp <- df_growth_pangenome %>% 
  filter(!stringr::str_detect(graph, "^malawi_")) %>% 
  #filter(key %in% c('bubble_density', 'bb_conserv_perc', 'bb_ubiq_perc') ) %>% 
  group_by(graph) %>% 
  filter(sample == max(sample)) %>% 
  pivot_wider(id_cols=c('graph','sample'), names_from='key', values_from='value')

# the most useful stats that would go into a table
tmp %>% reactable(columns = list(graph = colDef(sticky = "left")))

# all stats
tmp %>% 
  select(c(useful_stats, "variable_size_nonref", "bubble_density_corr", "bb_size")) %>% 
  reactable(columns = list(graph = colDef(sticky = "left")))
```

```{r}
# coverage final stats
df_growth_pangenome_malawi %>% 
  filter(key %in% c('bubble_density', 'bb_conserv_perc', 'bb_ubiq_perc') ) %>% 
  group_by(graph) %>% 
  filter(sample == max(sample)) %>% 
  pivot_wider(id_cols=c('graph','sample'), names_from='key', values_from='value') %>% 
  as.data.frame()

# individual coverage on astCal backbone
df_cov_bysample_malawi %>% 
  filter(graph == "malawi") %>% 
  select(backbone, species, bb_cov_bysample, bb_cov_perc_bysample) %>%
  reactable(columns = list(backbone = colDef(sticky = "left")))

# self coverage on own backbone
df_cov_bysample %>% 
  filter(rank==1, str_detect(graph, '^malawi')) %>%
  select(graph, species, bb_cov_bysample, bb_cov_perc_bysample) %>% 
  reactable(columns = list(graph = colDef(sticky = "left")))

# all stats, bb_cov and bb_ubiq: absolute values here
df_growth_pangenome_malawi %>% 
  filter(sample == 8) %>% select(-graph, -sample) %>% 
  pivot_wider(id_cols = key, names_from = backbone, values_from = value) %>% 
  reactable(columns = list(key = colDef(sticky = "left")))

# the most useful stats that would go into a table
df_growth_pangenome_malawi %>% 
  filter(sample == 8) %>% 
  filter(key %in% c(useful_stats, "variable_size_nonref", "bubble_density_corr", "bb_size")) %>% 
  select(-graph, -sample) %>% 
  pivot_wider(id_cols = key, names_from = backbone, values_from = value) %>% 
  reactable(columns = list(key = colDef(sticky = "left")))
```

Note that the empirical values may not necesarily add up. Might be caused by some corrections.

This is the stats for the graph after 5 samples (backbone + 4 non-reference). We probably only need the `astCal` one.
```{r}
# all stats, bb_cov and bb_ubiq: absolute values here
df_growth_pangenome_malawi %>% 
  filter(sample == 5) %>% select(-graph, -sample) %>% 
  pivot_wider(id_cols = key, names_from = backbone, values_from = value) %>% 
  reactable(columns = list(key = colDef(sticky = "left")))

# the most useful stats that would go into a table
df_growth_pangenome_malawi %>% 
  filter(sample == 5) %>% 
  filter(key %in% c(useful_stats, "variable_size_nonref", "bubble_density_corr", "bb_size")) %>% 
  select(-graph, -sample) %>% 
  pivot_wider(id_cols = key, names_from = backbone, values_from = value) %>% 
  reactable(columns = list(key = colDef(sticky = "left")))
```