---
title: "Preprocessing SVs"
output: html_document
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  cache = TRUE, cache.path = '.cache_knitr/preprocess/',
  fig.width = 7, fig.height = 4, dpi = 350,
  fig.path='figure/1-preprocess/',
  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 variation they are and how widespread/private alleles are.

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

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")
setwd(paste0(Sys.getenv("HOME"), "/code/malawi_transposon/script/R_SCRIPT/"))
```

## Data import and preprocessing
```{r}
df_raw <- read_csv(
  paste0(dir_sv, "malawi_haplochromines-summary.csv"),
    col_types = paste0('cciiiiiic', paste0(rep('ciicic', num_of_genomes), collapse=''), collapse=''),
    na = c("", "NA", ".")
)

# some chromosome columns need to have their entries renamed
df_raw$`troMau2.0_chr` <- gsub('troMau_', '', df_raw$`troMau2.0_chr`)
df_raw$`rhaChi1.0_chr` <- gsub('rhaChi_', '', df_raw$`rhaChi1.0_chr`)
df_raw$`aulStu5.0_chr` <- gsub('aulStu_', '', df_raw$`aulStu5.0_chr`)
df_raw$`copChr1.0_chr` <- gsub('copChr_', '', df_raw$`copChr1.0_chr`)
df_raw$`otoArg1.0_chr` <- gsub('otoArg_', '', df_raw$`otoArg1.0_chr`)
df_raw$`rhaChi2.0_chr` <- gsub('rhaChiTwo_', '', df_raw$`rhaChi2.0_chr`)

# add identifier column for the SV coordinates in the backbone
df_raw <- mutate(df_raw, id=paste0(chr, ':', start, '-', end)) %>% relocate(id)

# rename columns so species names are more concise (i.e. "astCal_start")
first_species_idx <- ncol(df_raw) - nrow(color_schema_species)*6 + 1
colnames(df_raw)[first_species_idx:ncol(df_raw)] <- paste(
  rep(color_schema_species$species_abrv, each=6), c('chr', 'start', 'end', 'path', 'len', 'strand'), sep='_'
)

glue("Read in a dataframe with {nrow(df_raw)} rows")
```

Some bubbles are a single path where two backbone segments are joined together.
These shouldn't be identified bubbles to start with, so we remove them.
```{r}
filter(df_raw, n_paths == 1) %>%  select(id, n_paths, segment_list)

nrow_before <- nrow(df_raw)
df_raw <- filter(df_raw, n_paths > 1)
glue("removed {nrow_before - nrow(df_raw)} bubbles with n_paths == 1, {nrow(df_raw)} remaining")
```

## Filtering bubbles --- Checkpoint 1
### SV calling success rate
To perform variant calling, `minigraph` realigns the genomes onto the genome graph. It is not always possible to call an allele at each polymorphic site or bubble. This happens because certain samples will not have read/assembly coverage for that particular region. 
```{r}
df_raw['n_uncalled'] <- df_raw %>% 
  select(ends_with('_path')) %>% 
  is.na() %>% 
  rowSums() 

glue('Out of the {sv_count} SVs from the raw dataframe,',
  ' {all_count} ({round(all_count/sv_count*100, 2)}%)',
  ' had alleles successfully called in all samples.',
  n_uncalled = summary(as.factor(df_raw$n_uncalled)),
  sv_count = sum(n_uncalled),
  all_count = n_uncalled[['0']]
)

glue("There are {nrow(filter(df_raw, n_uncalled > 0))} SVs where alleles were not called in all species.")
```

For a large majority of bubbles, alleles were called for all the studied species. The bar chart below depicts the success rate of calling variants. At this stage, it is wise to not discard any loci just because the variant calling is not complete, because there is still useful biological information in there.
```{r sv_calling_success, fig.width=4.5, fig.height=4.5}
ggplot(df_raw, aes(as.factor(n_uncalled))) + 
  geom_histogram(stat='count') +
  #ggtitle("No. of SVs where alleles are not called") + 
  xlab("uncalled samples (lower is better)") +  
  ylab("frequency") +
  scale_x_discrete(expand = expansion(mult = c(.1, .1))) +
  scale_y_continuous(expand = expansion(mult = c(.05, .1))) +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count") +
  theme(aspect.ratio=1)
```

However, the sites where the uncalled samples was `r num_of_genomes - 1` is quite contradictory, because you wouldn't expect them to be bubbles to start with. I suspect the reason is because variant calling is done through assembly mapping to graph, there might be regions where it is difficult to align properly and call variants. This seems to happen a lot in a repetitive region in chromosome 3. These variants we are likely to exclude from many analyses later on.
```{r}
# where are these weird alleles on?
filter(df_raw, n_uncalled == 7) %>% select(chr) %>% table()

# how many paths do they have?
filter(df_raw, n_uncalled == 7) %>% select(n_paths) %>% table()

# if we exclude them, how many true bubbles left?
glue("If these weird backbone only bubbles are excluded, we have {nrow(df_raw %>% filter(n_uncalled < 7))} left.")
```


### Same allele and/or same length?
There are some bubbles where the same allele is called for all the species, meaning that only one of the possible paths is used. This could indicate that that particular bubble might be an artifact caused by errors during multiple sequence alignment. Some of these are the bubbles where the number of uncalled samples is `r num_of_genomes - 1` (because only one allele is called.)

Here we add columns indicating whether the SVs in all `r num_of_genomes` species are of the same allele and of the same length.
```{r}
df_raw$allSameLengthBool <- rowSums(select(df_raw, ends_with('_len') ) == df_raw$astCal_len, na.rm = TRUE ) == (num_of_genomes - df_raw$n_uncalled)
df_raw$allSameAlleleBool <- rowSums(select(df_raw, ends_with('_path')) == df_raw$astCal_path, na.rm = TRUE) == (num_of_genomes - df_raw$n_uncalled)

glue("Same length : {sum(df_raw$allSameLengthBool)}\n", "Same path   : {sum(df_raw$allSameAlleleBool)}")
```

The same allele being called for all species could be an error caused by assembly-to-graph mapping or an artefactual bubble. These should probably be excluded for most analysis. Sites where the alleles are different but they have the same length could be substitutions or inversions. These should still be included for most of the analyses.

Here we regenerate the plot above to illustrate SV calling success. However, we drop the bubbles where all the paths are the same across the samples. This also inadvertently drops the bubbles where the number of uncalled samples was `r num_of_genomes - 1`.
```{r sv_calling_success_filtered, fig.width=3.8, fig.height=3.2}
df_raw %>% 
  filter(allSameAlleleBool == FALSE) %>% 
  ggplot(aes(as.factor(n_uncalled))) + 
  geom_histogram(stat='count') +
  #ggtitle(glue("No. of SVs where alleles are not called (excl fake bubbles)")) + 
  xlab("uncalled samples (lower is better)") +  
  ylab("frequency") +
  scale_x_discrete(expand = expansion(mult = c(.1, .1))) +
  scale_y_continuous(expand = expansion(mult = c(.05, .1))) +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count")
```

```{r}
glue("We started with {nrow(df_raw)} bubbles.")
glue("Excluding n_uncalled==7 bubbles, {nrow(df_raw %>% filter(n_uncalled < 7))} left.")
glue("Excluding the same path bubbles (which includes the previous category), {nrow(df_raw) - sum(df_raw$allSameAlleleBool)} left.")
glue("Excluding both, {nrow(df_raw %>% filter(n_uncalled < 7, allSameAlleleBool == FALSE))} left. This is the number I would report as true bubbles.")
```

### Bottleneck
For further analysis, it is probably more meaningful to constrain to bubbles where there is a decent number of individuals with known called variants, and that there are no fake bubbles. I decide on the arbitrary thresholds.
```{r}
df_raw2 <- df_raw %>% 
  filter(
    n_uncalled <= 3,            # at least 5 samples
    allSameAlleleBool == FALSE  # alleles must be different
  )

glue("removed {nrow(df_raw)-nrow(df_raw2)} loci, {nrow(df_raw2)} remain.")
```

## Filtering bubbles --- Checkpoint 2
### Number of theoretical paths
```{r path_histogram_theoretical, fig.width=6, fig.height=2.6}
ggplot(df_raw2, aes(n_paths, fill = n_paths > 8)) + 
  geom_histogram(bins=60) + 
  scale_x_continuous(trans='log2', 
    breaks = sapply(c(1, seq(4, 32, 4)), function(x) 2^x),
    labels = trans_format("log2", label_math(2^.x)) ) +
  scale_fill_manual(NULL,
    values = c("#3a6e9c", "#cc8227"), labels = c("fewer or equal 8", "more than 8")
  ) +
  xlab("number of theoretical paths (log2 scale)") +
  ylab("number of bubbles") +
  theme(
    legend.position = c(.96, .99),
    legend.justification = c("right", "top"),
    legend.box.just = "right",
    legend.margin = margin(6, 6, 6, 6),
    legend.text = element_text(size = rel(0.6)),
    legend.key.size=unit(0.035, "npc")
  )
```

It seems like most bubbles are relatively simple. Let's zoom into the leftmost section of the distribution.
```{r path_histogram_theoretical_zoom, fig.width=3.2, fig.height=2.8}
ggplot(filter(df_raw2, n_paths <= 8), aes(n_paths)) + 
  geom_histogram(stat='count', fill = "#3a6e9c") +
  xlab("number of theoretical paths") + 
  ylab(NULL) +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count", size = rel(2.9)) +
  scale_y_continuous(expand = expansion(mult = c(.05, .1)))
```

```{r path_histogram_theoretical_zoom_complex, fig.width=3.2, fig.height=2.8}
df_raw2 %>% 
  filter(n_paths>8) %>% 
  ggplot(aes(n_paths)) + 
  geom_histogram(bins=55, fill = "#cc8227") + 
  scale_x_continuous(trans='log2', 
    breaks = sapply(c(1, seq(4, 32, 4)), function(x) 2^x),
    labels = trans_format("log2", label_math(2^.x)) ) +
  xlab("number of theoretical paths (log2 scale)") +
  ylab(NULL)

df_raw2 %>% 
  filter(n_paths > 8) %>% 
  ggplot(aes(cut(n_paths, breaks=sapply(c(3,seq(4,32,4)), function(x) 2**x), include.lowest=TRUE) )) + 
  geom_histogram(stat='count', fill = "#cc8227") +
  scale_x_discrete(labels=expression(
    2^{3~"-"~4}, 2^{4~"-"~8}, 2^{8~"-"~12}, 2^{12~"-"~16}, 
    2^{16~"-"~20}, 2^{20~"-"~24}, 2^{24~"-"~28}, 2^{28~"-"~32})) +
  xlab("number of theoretical paths") + 
  ylab(NULL) +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count", size = rel(2.9)) +
  scale_y_continuous(expand = expansion(mult = c(.05, .1)))
```

```{r}
df_raw2 %>% 
  filter(n_paths > 8) %>% .$n_paths %>% log2() %>% 
  cut(breaks=c(3,seq(4,32,1)), include.lowest=TRUE) %>% 
  table()
```

The left-skewed distribution implies that bubbles are generally simple.
```{r, echo=FALSE}
glue('{count_under_8} ({perc_under_8}%) variants have up to 8 theoretical paths.',
  '\nBreaking down further, there are',
  ' {count2} ({perc2}%) 2-path and',
  ' {count3} ({perc3}%) 3-path variants.',
  df_filt = filter(df_raw2, n_paths <= 8),
  total_count = nrow(df_raw2),
  count_under_8 = nrow(df_filt),
  perc_under_8  = round(count_under_8 / total_count * 100, 2),
  count2 = sum(df_filt$n_paths == 2),
  perc2  = round(count2 / total_count * 100, 2),
  count3 = sum(df_filt$n_paths == 3),
  perc3  = round(count3 / total_count * 100, 2)
  )
```


### Number of biological paths / alleles
We add a column that counts the unique biological alleles present at each bubble.
```{r}
df_npathbio <- df_raw2 %>% 
  select(id, n_paths, ends_with(c('_path', '_len'))) %>% 
  pivot_longer(cols = ends_with(c('_path', '_len')),
    names_to = c('species', '.value'),
    names_pattern = "(.*)_(.*)") %>% 
  drop_na(path, len) %>%   # some species have no variants called
  group_by(id, path) %>% 
  count(path) %>% 
  group_by(id) %>% 
  summarise(n_paths_bio = n())

df_raw2 <- left_join(df_raw2, df_npathbio)
rm(df_npathbio)
```

This is visualised as a bar chart. In lower capacity are for bubbles where not all species have alleles successfully called.
```{r sv_calling_success_rawstage, fig.width=4, fig.height=4}
# all bubbles
df_raw2 %>%
  ggplot(aes(n_paths_bio)) +
  geom_histogram(stat='count') +
  #ggtitle("Allele diversity of SV bubbles") + 
  xlab("number of observed paths (alleles)") +  ylab("frequency") +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count") +
  scale_y_continuous(breaks = seq(0, 140000, 20000), expand = expansion(mult = c(.05, .1)) ) +
  theme(aspect.ratio = 1, legend.position = 'none')

# all bubbles, but those with unknown paths are highlighted
df_raw2 %>%
  ggplot(aes(n_paths_bio)) +
  geom_histogram(aes(alpha = n_uncalled == 0), stat='count') +
  #ggtitle("Allele diversity of SV bubbles") + 
  xlab("number of observed paths (alleles)") +  ylab("frequency") +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count") +
  scale_alpha_manual(values=c(1, 0.35), breaks=c(TRUE, FALSE)) +
  scale_y_continuous(breaks = seq(0, 140000, 20000), expand = expansion(mult = c(.05, .1)) ) +
  theme(aspect.ratio = 1, legend.position = 'none')

# only for bubbles where all species have known paths
df_raw2 %>%
  filter(n_uncalled == 0) %>% 
  ggplot(aes(n_paths_bio)) +
  geom_histogram(stat='count') +
  #ggtitle("Allele diversity of SV bubbles (where all species have paths)") + 
  xlab("number of observed paths (alleles)") +  ylab("frequency") +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count") +
  scale_y_continuous(breaks = seq(0, 140000, 20000), expand = expansion(mult = c(.05, .1)) ) +
  theme(aspect.ratio = 1, legend.position = 'none')
```

In tabular format.
```{r}
df_raw2 %>% group_by(n_paths_bio, n_uncalled) %>% summarise(n()) %>% 
  pivot_wider(id_cols=n_paths_bio, names_from=n_uncalled, 
              values_from=`n()`, names_prefix = "uncalled_") %>% 
  reactable()
```

### Quantifying bubble complexity
There is a strong correlation between the number of flexible segments and the log number of theoretical paths
```{r nsegment_vs_npaths}
ggplot(df_raw2, aes(x=n_segments-2, y=n_paths)) + 
  geom_point() +
  xlab("number of flexible segments") +
  ylab("number of theoretical paths") + 
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  ggpubr::stat_cor(label.x = 100, label.y.npc = "middle",)
```

It is expected that the theoretical complexity of the bubble will increase with the biological complexity. To attempt to quantify this, we want to first bin these values into categories to make them easier to visualise.

Theoretical complexity:
```{r}
vec_npaths_theoretical <- cut(df_raw2$n_paths, breaks=c(2,3,4,5,6,11,max(df_raw2$n_paths)+1), right=FALSE)
mapping_dict <- c('2', '3', '4', '5', '6-10', 'more than 10')
names(mapping_dict) <- levels(vec_npaths_theoretical)
df_raw2$n_paths_categorised <- vec_npaths_theoretical %>% dplyr::recode_factor(!!!mapping_dict)

table(df_raw2$n_paths_categorised)
```

Biological / observed complexity:
```{r}
vec_npaths_observed <- cut(df_raw2$n_paths_bio, breaks=c(2,3,4,5,6,7,8+1), right=FALSE)
mapping_dict <- c('2', '3', '4', '5', '6', 'more than 6')
names(mapping_dict) <- levels(vec_npaths_observed)
df_raw2$n_paths_bio_categorised <- vec_npaths_observed %>% dplyr::recode_factor(!!!mapping_dict)

table(df_raw2$n_paths_bio_categorised)
```

This plot shows how the number of observed paths (facets) relate to the number of theoretical paths and number of segments. We are particularly more interested in the simpler bubbles as they are more straightforward to interpret and are less likely to be artifacts.
```{r complexity_relations_obs_facet, fig.width=9, fig.height=5}
ggplot(df_raw2, aes(x=n_segments-2, fill=n_paths_categorised )) +
  facet_wrap(.~paste0('Number of observed paths: ', n_paths_bio_categorised), scales='free_y', ncol=3) +
  geom_histogram(position = 'identity', color='black', bins = 40) +
  scale_fill_brewer('Number of theoretical paths', palette = 'Pastel2') +
  scale_x_continuous("Number of flexible segments",
    trans='log2', breaks = sapply(seq(0, floor(log2(max(df_raw2$n_segments)-2))), function(x) 2^x) ) +
  ylab("Number of bubbles") +
  theme_bw() +
  theme(legend.position='top', legend.key.size=unit(0.015, "npc"),
        legend.title = element_text(size=9), legend.text = element_text(size=8)) 
```

The same information but with number of theoretical paths as the facets.
```{r complexity_relations_theo_facet, fig.width=9, fig.height=5}
ggplot(df_raw2, aes(x=n_segments-2, fill=n_paths_bio_categorised )) +
  facet_wrap(.~paste0('Number of theoretical paths: ', n_paths_categorised), scales='free_y') +
  geom_histogram(position = 'identity', color='black', bins = 40) +
  scale_fill_brewer('Number of observed paths', palette = 'Pastel2') +
  scale_x_continuous("Number of flexible segments",
    trans='log2', breaks = sapply(seq(0, floor(log2(max(df_raw2$n_segments)-2))), function(x) 2^x) ) +
  ylab("Number of bubbles") +
  theme_bw() +
  theme(legend.position='top', legend.key.size=unit(0.015, "npc"),
        legend.title = element_text(size=9), legend.text = element_text(size=8)) 
```

This suggests that the number of observed and theoretical paths give a good indication about how many segments there are in the bubble. It could be a good metric to decide which bubbles are simpler to look at.

### Length of alleles
Next, for each bubble, we obtain the shortest and longest allele length observed. We also calculate the ratio min / max allele length.

The argument is that the smaller the ratio between min / max allele length, it would be more interesting to follow up, as it is closer to being a complete insertion or deletion. If the minimum allele length is 0, this ratio would be 0 automatically.
```{r}
tmp_df <- df_raw2 %>% 
  select(id, ends_with("_len")) %>% 
  pivot_longer(cols = ends_with('_len'),
    names_to = c('species', '.value'),
    names_pattern = "(.*)_(.*)") %>% 
  select(id, len) %>% 
  group_by(id) %>% 
  drop_na() %>%    # remove species at each bubble with no variant called
  distinct() %>% 
  summarise(
    lenbio_max = max(len), 
    lenbio_min = min(len),
    lenbio_minmaxratio = min(len) / max(len)
  )
  
df_raw2 <- left_join(df_raw2, tmp_df, by='id') 
```

The min / max allele ratio shows a polynomial descending pattern. The rate of decrease is faster for simpler bubbles.
```{r minmaxratio_histogram, fig.width=9, fig.height=5.2}
ggplot(df_raw2, aes(x=lenbio_minmaxratio)) +
  facet_wrap(.~paste0('Number of observed paths: ', n_paths_bio_categorised), scales='free_y') +
  xlab("Ratio of min / max allele length") +
  ylab("Number of bubbles") +
  geom_histogram()
```

Like done previously, we bin the min / max allele length ratio into categories to aid visualisation. The categories are not even because this is a very skewed distribution. Note that there are `r filter(df_raw2, lenbio_minmaxratio == 1) %>% nrow()` inversion / substitution loci which are not included due to the `right = FALSE` parameter. This is intentional.
```{r}
# those which have a value of 1 are substitutions and inversions
df_raw2$lenbio_minmaxratio_categorised <- cut(df_raw2$lenbio_minmaxratio, 
  breaks=c(0,0.05,0.1,0.2,0.5,1.0), right=FALSE
)
table(df_raw2$lenbio_minmaxratio_categorised)
```

This plot shows how the number of observed paths (facets) relate to the minimum allele length and the min / max ratio. Notable observations:

- The leftmost peak is where the shortest allele is 0. There are more of these in simpler bubbles, as the relative height of that peak decreases as complexity increases.
- We can't safely exclude bubbles with a high number of observed or theoretical paths, because there might be some meaningful high quality indels (although the proportion is significantly lower).
- It does seem that number of observed path might be a better filtering mechanism, because the distribution falls off faster.
- Instead, we could consider a **cutoff threshold for the minimum allele length** of somewhere between 100 to 1000.

```{r complexity_vs_minmaxratio_minX, fig.width=9, fig.height=5}
# biological complexity as facets
df_raw2 %>% 
  filter(lenbio_minmaxratio!=1) %>%  # do not plot substitutions and inversions 
  ggplot(aes(x=lenbio_min+1, fill=lenbio_minmaxratio_categorised)) +
  facet_wrap(.~paste0('Number of observed paths: ', n_paths_bio_categorised), scales='free_y') +
  geom_histogram(bins=25, color='black', alpha = 0.75) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  scale_fill_brewer('Ratio of min / max allele length', palette = 'PuOr') +
  xlab("Minimum allele length") +
  ylab("Number of bubbles") +
  theme_bw() +
  theme(legend.position='top', legend.key.size=unit(0.015, "npc"),
        legend.title = element_text(size=9), legend.text = element_text(size=8))

# theoretical complexity as facets
df_raw2 %>% 
  filter(lenbio_minmaxratio!=1) %>%  # do not plot substitutions and inversions 
  ggplot(aes(x=lenbio_min+1, fill=lenbio_minmaxratio_categorised)) +
  facet_wrap(.~paste0('Number of theoretical paths: ', n_paths_categorised), scales='free_y') +
  geom_histogram(bins=25, color='black', alpha = 0.75) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  scale_fill_brewer('Ratio of min / max allele length', palette = 'RdGy') +
  xlab("Minimum allele length") +
  ylab("Number of bubbles") +
  theme_bw() +
  theme(legend.position='top', legend.key.size=unit(0.015, "npc"),
        legend.title = element_text(size=9), legend.text = element_text(size=8)) 
```

This plot shows how the number of observed paths (facets) relate to the maximum allele length and the min / max ratio.

- The maximum allele length distribution shifts rightwards as the bubble complexity increases.
- In the facet where `n_path_bio = 2` and `n_path = 2`, the leftmost peak where maximum allele length is 50 is caused by the minigraph `-L` parameter. There is a chance that many of these are not going to contain TEs, but they might still hold transposable elements. Or could be just noise.
- It does seem like the maximum allele length centers around 1000 to 10000, which is the expected TE size.

```{r complexity_vs_minmaxratio_maxX, fig.width=9, fig.height=5}
# biological complexity as facets
df_raw2 %>% 
  filter(lenbio_minmaxratio!=1) %>%  # do not plot substitutions and inversions 
  ggplot(aes(x=lenbio_max, fill=lenbio_minmaxratio_categorised)) +
  facet_wrap(.~paste0('Number of observed paths: ', n_paths_bio_categorised), scales='free_y') +
  geom_histogram(bins=25, color='black', alpha = 0.75) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  scale_fill_brewer('Ratio of min / max allele length', palette = 'PuOr') +
  xlab("Maximum allele length") +
  ylab("Number of bubbles") +
  theme_bw() +
  theme(legend.position='top', legend.key.size=unit(0.015, "npc"),
        legend.title = element_text(size=9), legend.text = element_text(size=8)) 

# theoretical complexity as facets
df_raw2 %>% 
  filter(lenbio_minmaxratio!=1) %>%  # do not plot substitutions and inversions 
  ggplot(aes(x=lenbio_max, fill=lenbio_minmaxratio_categorised)) +
  facet_wrap(.~paste0('Number of theoretical paths: ', n_paths_categorised), scales='free_y') +
  geom_histogram(bins=25, color='black', alpha = 0.75) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  scale_fill_brewer('Ratio of min / max allele length', palette = 'RdGy') +
  xlab("Maximum allele length") +
  ylab("Number of bubbles") +
  theme_bw() +
  theme(legend.position='top', legend.key.size=unit(0.015, "npc"),
        legend.title = element_text(size=9), legend.text = element_text(size=8)) 
```

### (Soft) Bottleneck to identify interesting SVs for PCR
Based on the above, some good soft filters would be:
```{r}
thresh_obspath <- 5
thresh_lenbiomin <- 1000
thresh_lenbio_minmaxratio <- 0.2
```

* Number of observed paths must be **`r thresh_obspath` or below**. At higher numbers (6 onwards), there is already a large proportion of bubbles with more than 10 theoretical paths: `r sum(df_raw2$n_paths_bio <= thresh_obspath)` bubbles (`r round(sum(df_raw2$n_paths_bio <= thresh_obspath) / nrow(df_raw2) * 100, 2)`%).
* The minimum allele length must be **smaller than `r thresh_lenbiomin`**, or 100 if prefer to be stricter: `r sum(df_raw2$lenbio_min < thresh_lenbiomin)` bubbles (`r round(sum(df_raw2$lenbio_min < 1000) / nrow(df_raw2) * 100, 2)`%).
* The ratio of min / max allele length should be **below `r thresh_lenbio_minmaxratio`**: `r sum(df_raw2$lenbio_minmaxratio < thresh_lenbio_minmaxratio)` bubbles (`r round(sum(df_raw2$lenbio_minmaxratio < thresh_lenbio_minmaxratio) / nrow(df_raw2) * 100, 2)`%).

At this stage, I would not filter anything out, but simply add a column stating whether a bubble is "simple".
```{r}
df_raw2$simple <- df_raw2$n_paths_bio <= thresh_obspath & df_raw2$lenbio_min < thresh_lenbiomin & df_raw2$lenbio_minmaxratio < thresh_lenbio_minmaxratio

glue("There are {count} bubbles classified as simple, which is {perc}%. This is the INTERSECT (not union) of the three filters above.",
  count = sum(df_raw2$simple),
  perc = round(count / nrow(df_raw2) * 100, 2)
)
```

Rearrange columns for saving.
```{r}
df_raw2 %>% select(id, chr, start, end, n_segments,
  n_paths, n_paths_bio, n_paths_categorised, n_paths_bio_categorised,
  inv, len_min, len_max, segment_list,
  n_uncalled, allSameLengthBool, simple,
  lenbio_min, lenbio_max, lenbio_minmaxratio, lenbio_minmaxratio_categorised,
  astCal_chr:rhaChi2_strand) %>% 
  write_csv(file = paste0(
    Sys.getenv("HOME"),
    "/code/malawi_transposon/local/script/malawi_haplochromines-summary-preprocessed.csv")
  )
```
