tigmint-parameters.rmd

---
title: Determine sensitivity to parameters of Tigmint
author: Shaun Jackman
output:
  html_document:
    keep_md: true
params:
  input_tsv:
    label: "Input TSV file of parameters and QUAST metrics"
    value: "tigmint-span-parameters.quast.tsv"
    input: text
  output_tsv:
    label: "Output TSV file of parameters and QUAST metrics"
    value: "tigmint-span-parameters.quast.out.tsv"
    input: text
---

```{r setup, message=FALSE}
library(cowplot)
library(dplyr)
library(ggplot2)
library(ggrepel)
library(knitr)
library(magrittr)
library(readr)
library(scales)
library(stringr)
library(tidyr)

knit_print.data.frame <- function(x, ...) kable(x) %>% paste(collapse = "\n") %>% asis_output
input_tsv <- params$input_tsv
output_tsv <- params$output_tsv

input_tsv
output_tsv
```

# Read the data

```{r read-data}
metrics_orig <- read_tsv(input_tsv)

metrics <- metrics_orig %>%
	distinct() %>%
	rename(Misassemblies = `# misassemblies`) %>%
	filter(!str_detect(Assembly, "arcs"), !str_detect(Assembly, "broken")) %>%
	mutate(Assembler = Assembly %>%
			str_replace("abyss2.*", "ABySS") %>%
			str_replace("discovardenovo_abyss.*", "DISCO+ABySS") %>%
			str_replace("discovardenovo_besst.*", "DISCO+BESST")) %>%
	bind_cols(as_tibble(str_match(.$Assembly, "window([0-9]+).*span([0-9]+)")) %>%
		transmute(Window = parse_integer(V2), Span = parse_integer(V3))) %>%
	mutate(Label = str_c("w=", Window, " s=", Span))
```

# Determine the range of the data
```{r range}
ranges <- metrics %>% summarize_at(vars(Window, Span, Misassemblies, NGA50), funs(min, max), na.rm = TRUE)
misassemblies_range <- ranges %$% c(0, Misassemblies_max)
nga50_range <- ranges %$% c(0, NGA50_max)
window_range <- ranges %$% c(Window_min, Window_max)
span_range <- ranges %$% c(Span_min, Span_max)
```

# Plot parameters and assembly metrics
```{r NGA50-misassemblies}
ggplot(drop_na(metrics)) +
	aes(x = Misassemblies, y = NGA50, label = Label) +
	geom_point() +
	geom_text_repel() +
	scale_x_continuous(name = "QUAST Misassemblies", labels = comma) +
	scale_y_continuous(name = "Scaffold NGA50", labels = unit_format(unit = "Mbp", scale = 1e-6)) +
	theme_minimal(base_size = 11)
```

# Plot misassemblies vs window
```{r misassemblies-window}
plot_misassemblies_window <- ggplot(metrics %>% filter(Span == 20)) +
	aes(x = Window, y = Misassemblies,
		colour = Assembler, group = Assembler, shape = Assembler) +
	geom_point() +
	geom_line() +
	scale_y_continuous(name = "Misassemblies", labels = comma) +
	scale_x_log10(name = "Window (bp)", labels = comma,
		breaks = unique(metrics$Window), minor_breaks = NULL) +
	scale_colour_brewer(palette = "Dark2") +
	expand_limits(x = window_range, y = misassemblies_range) +
	theme_minimal(base_size = 11) +
	theme(legend.position="bottom")
plot_misassemblies_window
```

# Plot NGA50 vs window
```{r NGA50-window}
plot_nga50_window <- ggplot(metrics %>% filter(Span == 20)) +
	aes(x = Window, y = NGA50,
		colour = Assembler, group = Assembler, shape = Assembler) +
	geom_point() +
	geom_line() +
	scale_y_continuous(name = "NGA50 (Mbp)", labels = unit_format(unit = "", scale = 1e-6)) +
	scale_x_log10(name = "Window (bp)", labels = comma,
		breaks = unique(metrics$Window), minor_breaks = NULL) +
	scale_colour_brewer(palette = "Dark2") +
	expand_limits(x = window_range, y = nga50_range) +
	theme_minimal(base_size = 11) +
	theme(legend.position="bottom")
plot_nga50_window
```

# Plot misassemblies vs spanning molecules
```{r misassemblies-molecules}
plot_misassemblies_molecules <- ggplot(metrics %>% filter(Window == 2000)) +
	aes(x = Span, y = Misassemblies,
		colour = Assembler, group = Assembler, shape = Assembler) +
	geom_point() +
	geom_line() +
	scale_y_continuous(name = "Misassemblies", labels = comma) +
	scale_x_log10(name = "Spanning molecules",
		breaks = unique(metrics$Span), minor_breaks = NULL) +
	scale_colour_brewer(palette = "Dark2") +
	expand_limits(x = span_range, y = misassemblies_range) +
	theme_minimal(base_size = 11) +
	theme(legend.position="bottom")
plot_misassemblies_molecules
```

# Plot NGA50 vs spanning molecules
```{r NGA50-molecules}
plot_nga50_molecules <- ggplot(metrics %>% filter(Window == 2000)) +
	aes(x = Span, y = NGA50,
		colour = Assembler, group = Assembler, shape = Assembler) +
	geom_point() +
	geom_line() +
	scale_y_continuous(name = "NGA50 (Mbp)", labels = unit_format(unit = "", scale = 1e-6)) +
	scale_x_log10(name = "Spanning molecules",
		breaks = unique(metrics$Span), minor_breaks = NULL) +
	scale_colour_brewer(palette = "Dark2") +
	expand_limits(x = span_range, y = nga50_range) +
	theme_minimal(base_size = 11) +
	theme(legend.position="bottom")
plot_nga50_molecules
```

# Metrics vs parameters
```{r parameters, fig.width = 6, fig.height = 4, dpi = 300}
plot_grid(ncol = 2,
	labels = c("A", "B", "C", "D"),
	plot_misassemblies_molecules + theme(legend.position = "none"),
	plot_misassemblies_window + theme(legend.position = "none"),
	plot_nga50_molecules + theme(legend.position = "none"),
	plot_nga50_window + theme(legend.position = "none")) %>%
	plot_grid(ncol = 1, rel_heights = c(2, 0.1), get_legend(plot_misassemblies_molecules))
```

# Table of parameters and assembly metrics
```{r parameters-table}
parameters_table <- metrics %>%
	arrange(Assembler, !is.na(Window), Window, Span) %>%
	replace_na(list(Window = "Original", Span = "Original")) %>%
	transmute(
		Assembler,
		Window,
		Span,
		`NG50 (Mbp)` = round(NG50 / 1e6, 2),
		`NGA50 (Mbp)` = round(NGA50 / 1e6, 2),
		Misassemblies = comma(Misassemblies))

parameters_table
write_tsv(parameters_table, output_tsv)
```