metrics.rmd

---
title: Plot and tabulate assembly metrics
author: Shaun Jackman
output:
  html_document:
    keep_md: true
params:
  abyss_fac_tsv:
    label: "Input TSV file of parameters and abyss-fac metrics"
    value: "abyss-fac.tsv"
    input: text
  samtobreak_tsv:
    label: "Input TSV file of parameters and abyss-samtobreak metrics"
    value: "samtobreak.tsv"
    input: text
  output_tsv:
    label: "Output TSV file of parameters and assembly metrics"
    value: "metrics.out.tsv"
    input: text
---

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

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

c(abyss_fac_tsv, samtobreak_tsv, output_tsv)
```

# Read the data

```{r read-data}
assemblers <- c("abyss", "discovardenovo-besst", "supernova")
assembly_names <- c(
	"ABySS",
	"ABySS + Tigmint",
	"ABySS + ARCS",
	"ABySS + Tigmint + ARCS",
	"DISCOVAR",
	"DISCOVAR + Tigmint",
	"DISCOVAR + ARCS",
	"DISCOVAR + Tigmint + ARCS",
	"Supernova",
	"Supernova + Tigmint",
	"Supernova + ARCS",
	"Supernova + Tigmint + ARCS")

abyss_fac_orig <- read_tsv(abyss_fac_tsv)
samtobreak_orig <- read_tsv(samtobreak_tsv)

assemblers_regex <- paste0(assemblers, collapse = "|")
abyss_fac <- abyss_fac_orig %>%
	filter(grepl(assemblers_regex, name), !grepl("scaftigs", name)) %>%
	add_column(Assembly = assembly_names)
samtobreak <- samtobreak_orig %>%
	filter(grepl(assemblers_regex, File)) %>%
	add_column(Assembly = assembly_names)

metrics <- left_join(abyss_fac, samtobreak, by = "Assembly") %>%
	separate(Assembly, "Assembler", extra = "drop", remove = FALSE)
```

# Plot NGA50 and breakpoints
```{r nga50-breakpoints}
ggplot(metrics) +
	aes(x = Total_breakpoints, y = Scaffold_NGA50, label = Assembly) +
	geom_point() +
	geom_text_repel(aes(colour = Assembler), segment.alpha = 0.5) +
	scale_x_continuous(name = "Breakpoints", labels = comma) +
	scale_y_continuous(name = "Scaffold NGA50", labels = unit_format(unit = "Mbp", scale = 1e-6)) +
	scale_colour_brewer(palette = "Dark2", guide = FALSE) +
	expand_limits(x = c(4000, 8000), y = c(0, 10e6)) +
	theme_minimal(base_size = 20)
```

# Plot NG50, NGA50, and breakpoints
```{r nga50-ng50-breakpoints}
ggplot(metrics) +
	aes(x = Total_breakpoints, y = NG50, label = Assembly) +
	geom_errorbar(aes(ymin = Scaffold_NGA50, ymax = NG50), width = 50) +
	geom_text_repel(nudge_x = -250, nudge_y = 500e3, segment.alpha = 0.5) +
	scale_x_continuous(name = "Breakpoints", labels = comma) +
	scale_y_continuous(name = "Scaffold NGA50-NG50", labels = unit_format(unit = "Mbp", scale = 1e-6)) +
	expand_limits(x = c(2000, 4000), y = c(0, 25e6)) +
	theme_minimal(base_size = 20)
```

# Table of assembly metrics
```{r metrics-table}
metrics_table <- metrics %>%
	transmute(
		Assembly,
		`NG50 (Mbp)` = round(Scaffold_NG50 / 1e6, 2),
		`NGA50 (Mbp)` = round(Scaffold_NGA50 / 1e6, 2),
		Breakpoints = comma(Total_breakpoints),
		Reduction = lag(Total_breakpoints) - Total_breakpoints,
		Reduction = str_c(Reduction, " (", percent(Reduction / lag(Total_breakpoints)), ")"),
		Reduction = ifelse(str_detect(Assembly, "Tigmint"), Reduction, NA))

metrics_table
write_tsv(metrics_table, output_tsv)
```