---
title: "nf-core/genomeassembler report"
author: ""
format: dashboard
editor: source
nav-buttons:
    - icon: github
    - href: https://github.com/nf-core/genomeassembler
params:
    nanoq: false
    busco: false
    quast: false
    jellyfish: false
    merqury: false
---

```{r load libraries and functions}
#| message: false
#| output: false
#| include: false

# Load libraries
library(tidyverse)
library(magrittr)
library(plotly)
# Load functions
list.files("functions", full.names = T, pattern = ".R") %>%
    map(\(x) source(x))

# Set default ggplot theme
theme_set(theme_bw(base_size = 14,
                    base_family = "Arial"))
theme_update(strip.background = element_blank(),
                axis.text.x = element_text(angle = 70, hjust = 1))
## Colors, these come from the khroma package ("muted")
### For <=9 stages:
color_scale_plots <- scale_color_manual(values = c("#CC6677", "#332288", "#DDCC77", "#117733", "#88CCEE", "#882255", "#44AA99", "#999933", "#AA4499"), na.value = "#DDDDDD")
fill_scale_plots <- scale_fill_manual(values = c("#CC6677", "#332288", "#DDCC77", "#117733", "#88CCEE", "#882255", "#44AA99", "#999933", "#AA4499"), na.value = "#DDDDDD")
# Base directory containing reports
data_base = "data/"
```

# About

This report displays the main information gathered from various QC steps.

# nanoq {.tabset}

::: {.content-visible unless-profile="nanoq"}
nanoq was not included in the pipeline run, no ONT reads were included.
:::

```{r nanoq read inputs}
#| eval: !expr params$nanoq
#| include: false
#| message: false
#| output: false

# Parse nanoq reports into table
nanoq_reports <- list.files(paste0(data_base, "nanoq"),
                            pattern = "report.json",
                            full.names = T) %>%
    map_dfr(\(x) read_nanoq(x))
```

```{r}
#| eval: !expr params$nanoq
#| include: false

# For each sample, we create one plot chunk that will be saved into nanoq files
# This is an rmd chunk in plain text.

dir.create("nanoq_files")
for (i in 1:length(unique(nanoq_reports$sample))) {
paste0('```{r}\n
            #| title: "Nanoq read statistics"
            p <- nanoq_reports %>%
            filter(stat %in% c("Median Length", "Longest", "Median Quality","Bases")) %>%
            filter(sample == "', unique(nanoq_reports$sample)[i], '") %>%
            mutate(stat=fct_relevel(stat,c("Bases","Longest","Median Length","Median Quality"))) %>%
            ggplot(aes(x = sample, y = val)) +
            geom_line() +
            geom_point(size = 5, pch=21, aes(fill=stage)) +
            facet_wrap(~stat, scales = "free_y", ncol=2) +
            fill_scale_plots +
            scale_y_continuous(labels = function(x) format(x,scientific=-1,trim=T, digits = 3, drop0trailing=T), n.breaks = 4) +
            theme(axis.title.x = element_blank(),
                    axis.title.y = element_blank(),
                    legend.position = "none",
                    legend.title = element_blank(),
                    panel.grid.minor = element_blank())
        ggplotly(p)\n```') %>%
    write_lines(glue::glue("nanoq_files/_{ unique(nanoq_reports$sample)[i] }_nanoq.Rmd"))
}
```

::: {.content-visible when-profile="nanoq"}

::: {.panel-tabset .flow}

```{r nanoq add subplots}
#| eval: !expr params$nanoq
#| results: asis

# This loop creates one tab per sample.
## Each tab contains 3 valueboxes
## Below the valueboxes, the sample-specific plot code generated above is inserted

for (i in 1:length(unique(nanoq_reports$sample))) {
    cat(paste0('## ', unique(nanoq_reports$sample)[i], '\n\n'),
        paste0('### { width = 30% }', '\n\n'),
        paste0('::: {.valuebox icon="magic" color="primary" title="Total bases sequenced"}','\n'),
        paste0(nanoq_reports %>%
            filter(stat == "Bases") %>%
            filter(sample == unique(nanoq_reports$sample)[i],) %$%
            sum(val) %>%
            format(scientific=-1,trim=T, digits = 3, drop0trailing=T),'\n'),
        paste0(':::', '\n\n'),
        paste0('::: {.valuebox icon="collection" color="secondary" title="Number of reads"}', '\n'),
        paste0(nanoq_reports %>%
                filter(stat == "N Reads") %>%
                filter(sample == unique(nanoq_reports$sample)[i]) %$%
                min(val) %>%
                paste(" bases"), '\n'),
        paste0(':::', '\n\n'),
        paste0('::: {.valuebox icon="chevron-double-up" color="success" title="Longest read"}', '\n'),
        paste0(nanoq_reports %>%
                filter(stat == "Longest") %>%
                filter(sample == unique(nanoq_reports$sample)[i]) %$%
                max(val) %>%
                paste0(" bases"),'\n'),
        paste0(':::', '\n\n'),
        paste0('### ', '\n\n'),
        knitr::knit_child(glue::glue('nanoq_files/_{ unique(nanoq_reports$sample)[i] }_nanoq.Rmd'),
                        envir = globalenv(),
                        quiet = TRUE),
        paste0('\n\n'),
        sep = ""
    )
}
```

```{r}
#| eval: !expr params$nanoq
# Clean up the intermediate files
unlink("nanoq_files", recursive = T)
```

:::

:::

# QUAST {.tabset}

::: {.content-visible unless-profile="quast"}
QUAST was not included in the pipeline run.
:::

::: {.content-visible when-profile="quast"}
QUAST reports assembly statistics, taking into account the reference, if provided.

```{r message = F}
#| eval: !expr params$quast
# This chunk parses the quast reports from data/quast
quast_stats <- list.files(paste0(data_base, "quast"),
                            pattern = "report.tsv",
                            full.names = T) %>%
    map_dfr(\(x) {
    read_quast_report(x) %>%
    mutate(sample = str_extract(x %>% basename(),
                        ".+?(?=_[assembly|links|longstitch|ragtag|medaka|pilon])"),
            stage = case_when(
                str_detect(x, "_ragtag") ~ "RagTag",
                str_detect(x, "_medaka") ~ "medaka",
                str_detect(x, "_pilon") ~ "pilon",
                str_detect(x, "_longstitch") ~ "longstitch",
                str_detect(x, "_links") ~ "LINKS",
                str_detect(x, "assembly") ~ "Assembly",
                TRUE ~ "Unknown")
            )
        }
    )
```

```{r quast write length plots}
#| eval: !expr params$quast
#| include: false
# This creates code that will generate the length plot based on the contents of the quast report.
dir.create("quast_files")
dir.create("quast_files/length")
for (i in 1:length(unique(quast_stats$sample))) {
paste0('```{r}\n
    p <- quast_stats %>%
        filter(sample == "', unique(quast_stats$sample)[i], '") %>%
        filter(str_detect(stat, "[L].*[59]0")) %>%
        mutate(stat = fct_relevel(stat, "L50","L90","LG50","LG90")) %>%
        ggplot(aes(x=stat, y=value)) +
        geom_point(aes(fill = stage),
            size = 5,
            pch = 21,
            alpha = 0.8,
            position = position_dodge(width = 0.4)) +
        facet_wrap(~ sample, scales = "free_y") +
        fill_scale_plots +
        labs(title = "QUAST: L(G) 50 and 90") +
        theme(panel.border = element_rect(fill = NA))
        ggplotly(p) \n```') %>%
    write_lines(glue::glue("quast_files/length/_{ unique(quast_stats$sample)[i] }_quast.Rmd"))
}
```

```{r quast contig plots}
#| eval: !expr params$quast
#| include: false
# This creates code that will generate the contig plots based on the contents of the quast report.

dir.create("quast_files/contigs")
for (i in 1:length(unique(quast_stats$sample))) {
paste0('```{r}\n
    p <- quast_stats %>%
    filter(sample == "', unique(quast_stats$sample)[i], '") %>%
    filter(str_detect(stat, "# contigs \\\\(")) %>%
    filter(!str_detect(stat, ">= 0")) %>%
    mutate(stat = stat %>% str_remove_all("# contigs ") %>% str_remove_all("[()]") %>% fct_inorder()) %>%
    ggplot(aes(x=stat, y=value)) +
    geom_point(aes(fill = stage),
            size = 5,
            pch = 21,
            alpha = 0.8,
            position = position_dodge(width = 0.4)) +
    facet_wrap(~ sample, scales = "free_y") +
    fill_scale_plots +
    theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
    labs(title = "QUAST: Number of contigs by size")
ggplotly(p)
    p <- quast_stats %>%
    filter(sample == "', unique(quast_stats$sample)[i], '") %>%
    filter(str_detect(stat, "Total length")) %>%
    filter(!str_detect(stat, ">= 0")) %>%
    mutate(stat = stat %>% str_remove_all("Total length ") %>% str_remove_all("[()]") %>% fct_inorder()) %>%
    ggplot(aes(x = stat, y = value)) +
    geom_point(
    aes(fill = stage),
    size = 5,
    pch = 21,
    height = 0,
    width = 0.2,
    alpha = 0.8,
    position = position_dodge(width = 0.4)
    ) +
    facet_wrap( ~ sample, scales = "free_y") +
    fill_scale_plots +
    theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
    labs(title = "QUAST: Aggregated length") +
    scale_y_continuous(
    labels = function(x)
        format(
            x,
            scientific = -1,
            trim = T,
            digits = 3,
            drop0trailing = T
        )
    )
ggplotly(p)
    \n```') %>%
    write_lines(glue::glue("quast_files/contigs/_{ unique(quast_stats$sample)[i] }_quast.Rmd"))
}
```

::: {.panel-tabset .flow}
```{r quast add length subplots}
#| eval: !expr params$quast
#| results: asis
# This generates the tab-page for each sample
# Per sample there are 3 value boxes
# Below the value boxes there are two plots, one showing the length and one showing the contig statistics

for (i in 1:length(unique(quast_stats$sample))) {
    cat(paste0('## ', unique(quast_stats$sample)[i], '\n'),
        paste0('### { width=30% }\n\n'),
        paste0('::: {.valuebox icon="arrow-up-right-circle" color="primary" title="Total length"}\n'),
        quast_stats %>%
            filter(sample == unique(quast_stats$sample)[i]) %>%
            filter(stat == "Total length (>= 0 bp)") %$%
            max(value) %>%
            format(
                scientific = -1,
                trim = T,
                digits = 3,
                drop0trailing = T
            ) %>%
            paste("bp"),
        paste0('\n'),
        paste0(':::'),
        paste0('\n\n'),
        paste0('::: {.valuebox icon="percent" color="success" title="GC Content"}\n'),
        quast_stats %>%
            filter(sample == unique(quast_stats$sample)[i]) %>%
            filter(stat == "GC (%)") %$%
            mean(value) %>%
            round(2) %>%
            paste(" %"),
        paste0('\n'),
        paste0(':::'),
        paste0('\n\n'),
        paste0('::: {.valuebox icon="emoji-heart-eyes" color="info" title="Lowest L90"}\n'),
        quast_stats %>%
            filter(sample == unique(quast_stats$sample)[i]) %>%
            filter(stat == "L90") %>%
            filter(value == min(value)) %>%
            unique() %$%
            glue::glue("{unique(value)}, at stage(s): {paste(stage, collapse = ', ')}"),
        paste0('\n'),
        paste0(':::'),
        paste0('\n\n'),
        paste0('### {.tabset}'),
        paste0('\n\n'),
        paste0('#### Tables \n\n'),
        quast_stats %>%
            filter(sample == unique(quast_stats$sample)[i]) %>%
            dplyr::select(sample, stage, stat, value) %>%
            pivot_wider(names_from = "stat", values_from = "value") %>%
            #knitr::kable(format = 'html', caption = glue::glue('QUAST statistics'))
            gt::gt() %>%
            gt::cols_nanoplot(columns = starts_with("# contigs ("),
                    new_col_name = "Contigs_by_size",
                    new_col_label = gt::md("*# Contigs by size*")) %>%
            gt::cols_nanoplot(columns = starts_with("Total length ("),
                    new_col_name = "Total_length",
                    new_col_label = gt::md("*Total length*")) %>%
            gt::tab_footnote(
                    footnote = "Breaks are: contigs >= 0, 1kb, 5kb, 10kb, 25kb, 50kb",
                    locations = gt::cells_column_labels(columns = c(Contigs_by_size, Total_length))) %>%
            gt::cols_align(align = "center",
                    columns = c(Contigs_by_size, Total_length)) %>%
            gt::cols_move(Contigs_by_size, "Largest contig") %>%
            gt::cols_move(Total_length, "Total length") %>%
            gt::as_raw_html(),
        paste0('\n\n'),
        paste0('#### Plots \n\n'),
        knitr::knit_child(glue::glue('quast_files/length/_{ unique(quast_stats$sample)[i] }_quast.Rmd'),
                            envir = globalenv(),
                            quiet = TRUE),
        paste0('\n\n\n'),
        knitr::knit_child(glue::glue('quast_files/contigs/_{ unique(quast_stats$sample)[i] }_quast.Rmd'),
                        envir = globalenv(),
                        quiet = TRUE),
        paste0('\n\n\n'),
    sep = "")
}
```

```{r unlink quast}
#| eval: !expr params$quast
# Remove temporary files, write out collected report (mainly for debugging)
unlink("quast_files/contigs", recursive = T)
unlink("quast_files/length", recursive = T)
write_csv(quast_stats,"quast_files/reports.csv")
```
:::
:::

# BUSCO

::: {.content-visible unless-profile="busco"}
BUSCO was not included in the pipeline run.
:::

```{r}
#| eval: !expr params$busco
#| warning: false
#| message: false
#| echo: false
# Parse the reports from busco
busco_reports <- list.files(paste0(data_base, "busco"),
                            full.names = T,
                            pattern = "batch_summary") %>%
    map_dfr(\(x) read_busco_batch(x))
```


```{r}
#| eval: !expr params$busco
#| include: false
# This creates code that will generate the plots based on BUSCO results

dir.create("busco_files")
dir.create("busco_files/orthologs")
for (i in 1:length(unique(busco_reports$sample))) {
paste0('```{r}\n
            p <- busco_reports %>%
                filter(sample == "', unique(busco_reports$sample)[i], '") %>%
                filter(Var %in% c("Complete","Single","Duplicated","Fragmented")) %>%
                ggplot(aes(y = value, x = Var)) +
                geom_point(
                    aes(fill = stage),
                    size = 6,
                    pch = 21,
                    height = 0,
                    alpha = 0.8,
                    position = position_dodge(width = 0.4)
                ) +
                facet_wrap( ~ sample, nrow = 3) +
                fill_scale_plots +
                labs(   y = "% of Single Copy Orthologs",
                        title = "BUSCO: Conserved Orthologs") +
                coord_cartesian(clip = "on") +
                theme(
                    panel.border = element_rect(fill = NA),
                    legend.position = "bottom",
                    axis.title.y = element_text(angle = 90),
                    axis.title.x = element_blank()
                )
                ggplotly(p)
                \n```') %>%
    write_lines(glue::glue("busco_files/orthologs/_{ unique(busco_reports$sample)[i] }_orthologs.Rmd"))
}
```

::: {.content-visible when-profile="busco"}
BUSCO assess assembly quality based on the presence / absence of expected single-copy orthologs.

::: {.panel-tabset .flow}
```{r busco orthologs add subplots and valueboxes}
#| eval: !expr params$busco
#| results: asis
#
# This generates the tab-page for each sample
# Per sample there are 3 value boxes
# Below the value boxes there are one plots, showing the BUSCO statistics

for (i in 1:length(unique(busco_reports$sample))) {
    cur_sample <-  unique(busco_reports$sample)[i]
    # The BUSCO valueboxes contain information on which stage of the assembly had the highest quality, this requires some variables.
    completeness_val <-  busco_reports %>%
        filter(sample == cur_sample) %>%
        filter(Var == "Complete") %>%
        filter(value == max(value)) %$%
        value %>%
        unique()
    completeness_stage  <-  busco_reports %>%
        filter(sample == cur_sample) %>%
        filter(Var == "Complete") %>%
        filter(value == max(value)) %$%
        stage
    frag_val <-  busco_reports %>%
        filter(sample == cur_sample) %>%
        filter(Var == "Fragmented") %>%
        filter(value == max(value)) %$%
        value %>%
        unique()
    frag_stage  <-  busco_reports %>%
        filter(sample == cur_sample) %>%
        filter(Var == "Fragmented") %>%
        filter(value == max(value)) %$%
        stage
    missing_val <-  busco_reports %>%
        filter(sample == cur_sample) %>%
        filter(Var == "Missing") %>%
        filter(value == max(value)) %$%
    value %>%
        unique()
    missing_stage  <-  busco_reports %>%
        filter(sample == cur_sample) %>%
        filter(Var == "Missing") %>%
        filter(value == max(value)) %$%
    stage
    cat(paste('## ', unique(busco_reports$sample)[i]),
        paste0('\n\n'),
        paste0('### {.fill} \n\n'),
        paste0('::: {.valuebox icon="percent" color="success" title="Max. BUSCO Completeness" }\n'),
        paste0('\n'),
        glue::glue("{unique(completeness_val)}%,\nat stage(s): {paste(unique(completeness_stage), collapse = ', ')}"),
        paste0('\n'),
        paste0(':::'),
        paste0('\n'),
        paste0('::: {.valuebox icon="heartbreak" color="warning" title="Max. BUSCO Fragmented"}\n'),
        glue::glue("{unique(frag_val)}%,\nat stage(s): {paste(unique(frag_stage), collapse = ', ')}"),
        paste0('\n'),
        paste0('\n'),
        paste0(':::'),
        paste0('\n'),
        paste0('::: {.valuebox icon="person-walking" color="danger" title="Max. BUSCOs Missing"}\n'),
        glue::glue("{unique(missing_val)}%, at stage(s): {paste(unique(missing_stage), collapse = ', ')}"),
        paste0('\n'),
        paste0(':::'),
        paste0('\n\n'),
        paste('###'),
        paste0('\n\n'),
        busco_reports %>%
            filter(sample == cur_sample) %>%
            dplyr::select(sample, stage, Var, value) %>%
            mutate(Var = str_replace_all(Var, "_", " ")) %>%
            pivot_wider(names_from = "Var", values_from = "value") %>%
            gt::gt() %>%
            gt::as_raw_html(),
        paste0('\n\n'),
        knitr::knit_child(glue::glue('busco_files/orthologs/_{ unique(busco_reports$sample)[i] }_orthologs.Rmd'),
                            envir = globalenv(),
                            quiet = TRUE),
        paste0('\n\n\n'),
        sep = "")
}
```
:::
:::

```{r}
#| eval: !expr params$busco
# Delete temporary files
unlink("busco_files/orthologs", recursive = T)
```

```{r}
#| eval: !expr params$busco
# Export large report table, mainly for debugging
write_csv(busco_reports,"busco_files/reports.csv")
```

# genomescope

::: {.content-visible unless-profile="jellyfish"}
jellyfish / genomescope was not included in the pipeline run.
:::

```{r}
#| eval: !expr params$jellyfish
#| message: false
#| echo: false
#| output: false
#| warning: false
# Parse the genomescope statistics
genomescope_out <- list.files(paste0(data_base, "genomescope"), full.names = T, pattern = "genomescope.txt") %>%
    map_dfr(\(x) read_genomescope(x))
```

::: {.content-visible when-profile="jellyfish"}
Jellyfish and genomescope are used to infer genome size from the initial ONT reads.

```{r}
#| eval: !expr params$jellyfish
#| output: asis
# Since genomescope produces plots, I am simply including those here instead of recreating them, the proper QC for kmers comes with merqury.
img_files <- list.files(paste0(data_base,"genomescope"), full.names = T, pattern = "plot.png")
dir.create("genomescope_files")
for (file in img_files) {
    file.copy(from = file,
            to   = paste0("genomescope_files/", file %>% basename(), sep =""))

}
img_files <- list.files("genomescope_files", full.names = T, pattern = "plot.png")

cat(":::{.panel-tabset}\n",
    glue::glue('## <<str_extract(img_files %>% basename(), ".+?(?=_plot.png)")>>\n ![](<<img_files>>){width=50% fig-align="centre"}\n\n\n', .open = "<<", .close = ">>"),
    ":::\n",
    sep = ""
)
```
:::

# merqury

::: {.content-visible unless-profile="merqury"}
meryl and merqury were not included in the pipeline run.
:::

```{r}
#| eval: !expr params$merqury
#| include: false
#| message: false
#| output: false

# Here the merqury stats are parsed and the assembly stage is extracted
merqury_stats <- list.files(paste0(data_base, "merqury"), full.names = T, pattern = "stats") %>%
    lapply(\(x) {
    read_tsv(x, col_names = c("sample_stage","all","assembly","total","percent"), show_col_types = FALSE) %>%
    mutate(sample = str_extract(x %>% basename(),
                                ".+?(?=_[assembly|links|longstitch|ragtag|medaka|pilon])"),
            stage = case_when(
                    str_detect(x, "_ragtag") ~ "RagTag",
                    str_detect(x, "_medaka") ~ "medaka",
                    str_detect(x, "_pilon") ~ "pilon",
                    str_detect(x, "_longstitch") ~ "longstitch",
                    str_detect(x, "_links") ~ "LINKS",
                    str_detect(x, "assembly") ~ "Assembly",
                    TRUE ~ "Unknown"
                    )
            )
        }
    ) %>%
    bind_rows()
# This parses the assembly stats
merqury_asm_hists <- list.files(paste0(data_base, "/merqury"), full.names = T, pattern = "asm.hist")  %>%
    lapply(\(x) {
            read_tsv(x, col_names = T, show_col_types = FALSE) %>%
            mutate(
                sample = str_extract(x %>% basename(),
                                ".+?(?=_[assembly|links|longstitch|ragtag|medaka|pilon])"),
                stage = case_when(
                        str_detect(x, "_ragtag") ~ "RagTag",
                        str_detect(x, "_medaka") ~ "medaka",
                        str_detect(x, "_pilon") ~ "pilon",
                        str_detect(x, "_longstitch") ~ "longstitch",
                        str_detect(x, "_links") ~ "LINKS",
                        str_detect(x, "assembly") ~ "Assembly",
                        TRUE ~ "Unknown"),
                Assembly = as.factor(Assembly),
                stage = as.factor(stage),
                sample = as.factor(sample),
                kmer_multiplicity = as.integer(kmer_multiplicity),
                Count = as.integer(Count)
            )
        }
    ) %>%
    bind_rows()
# This parses the copy number file
merqury_cn_hists <- list.files(paste0(data_base, "merqury"), full.names = T, pattern = "cn.hist")  %>%
    lapply(\(x) {
            read_tsv(x, col_names = T, show_col_types = FALSE) %>%
            mutate(
                sample = str_extract(x %>% basename(),
                                ".+?(?=_[assembly|links|longstitch|ragtag|medaka|pilon])"),
                stage = case_when(
                    str_detect(x, "_ragtag") ~ "RagTag",
                    str_detect(x, "_medaka") ~ "medaka",
                    str_detect(x, "_pilon") ~ "pilon",
                    str_detect(x, "_longstitch") ~ "longstitch",
                    str_detect(x, "_links") ~ "LINKS",
                    str_detect(x, "assembly") ~ "Assembly",
                    TRUE ~ "Unknown"),
                Copies = as.factor(Copies),
                stage = as.factor(stage),
                sample = as.factor(sample),
                kmer_multiplicity = as.integer(kmer_multiplicity),
                Count = as.integer(Count)
            )
        }
    ) %>%
    bind_rows()
# This parses the qv file
merqury_qv <-
    list.files(paste0(data_base, "merqury"), full.names = T, pattern = ".qv") %>%
        lapply(\(x) {
                read_tsv(x,
                    col_names = c("Assembly", "kmers_assembly_unique", "kmers_assembly_shared", "QV", "error_rate"),
                    show_col_types = FALSE) %>%
                mutate(
                    sample = str_extract(x %>% basename(),
                                ".+?(?=_[assembly|links|longstitch|ragtag|medaka|pilon])"),
                    stage = case_when(
                        str_detect(x, "_ragtag") ~ "RagTag",
                        str_detect(x, "_medaka") ~ "medaka",
                        str_detect(x, "_pilon") ~ "pilon",
                        str_detect(x, "_longstitch") ~ "longstitch",
                        str_detect(x, "_links") ~ "LINKS",
                        str_detect(x, "assembly") ~ "Assembly",
                        TRUE ~ "Unknown"),
                    stage = as.factor(stage),
                    sample = as.factor(sample),
                    kmers_assembly_shared = as.integer(kmers_assembly_shared),
                    kmers_assembly_unique = as.integer(kmers_assembly_unique),
                    QV = as.double(QV),
                    error_rate = as.double(error_rate)
                )
            }
        ) %>%
    bind_rows()
dir.create("merqury_files")
```

```{r merqury qv}
#| eval: !expr params$merqury
#| include: false
# This generates QV-plots from merqury; the plot function is stuffed into plot_merqury
dir.create("merqury_files/qv_plots/")
for (i in 1:length(unique(merqury_qv$sample))) {
    cur_sample <- unique(merqury_qv$sample)[i]
    paste0('```{r}
p <- merqury_qv %>%
    plot_merqury_qv("', cur_sample,'")
ggplotly(p)\n```') %>%
    write_lines(glue::glue("merqury_files/qv_plots/_{ cur_sample }_qv_plt.Rmd"))
}
```

```{r merqury completeness}
#| eval: !expr params$merqury
#| include: false
# This generates stat-plots from merqury; the plot function is stuffed into plot_merqury

dir.create("merqury_files/stat_plots/")
for (i in 1:length(unique(merqury_stats$sample))) {
    cur_sample <- unique(merqury_stats$sample)[i]
    paste0('```{r}
p <- merqury_stats %>%
    plot_merqury_stats("', cur_sample,'")
ggplotly(p)\n```') %>%
    write_lines(glue::glue("merqury_files/stat_plots/_{ cur_sample }_completeness_plt.Rmd"))
}
```

```{r merqury asm}
#| eval: !expr params$merqury
#| include: false
# This generates assembly plots from merqury; the plot function is stuffed into plot_merqury

dir.create("merqury_files/asm_plots/")
for (i in 1:length(unique(merqury_asm_hists$sample))) {
    cur_sample <- unique(merqury_asm_hists$sample)[i]
    paste0('```{r}
p <- merqury_asm_hists %>%
        plot_merqury_multiplicity("', cur_sample,'")
ggplotly(p)\n```') %>%
    write_lines(glue::glue("merqury_files/asm_plots/_{ cur_sample }_asm_plt.Rmd"))
}
```

```{r merqury cn}
#| eval: !expr params$merqury
#| include: false
# This generates copy-number from merqury; the plot function is stuffed into plot_merqury

dir.create("merqury_files/cn_plots/")
for (i in 1:length(unique(merqury_cn_hists$sample))) {
    cur_sample <- unique(merqury_cn_hists$sample)[i]
    paste0('```{r}
p <- merqury_cn_hists %>%
        plot_merqury_copynumber("', cur_sample,'")
ggplotly(p)\n```') %>%
    write_lines(glue::glue("merqury_files/cn_plots/_{ cur_sample }_cn_plt.Rmd"))
}
```

::: {.content-visible when-profile="merqury"}
merqury compares k-mer spectra between assemblies and short read libraries to assess assembly quality and completeness.

::: {.panel-tabset .flow}
```{r merqury add plots and valueboxes}
#| eval: !expr params$merqury
#| results: asis
# This generates the tab-page for each sample
# Per sample there are 3 value boxes
# Below the value boxes there is a tabset of plots, each tab contains one of the plot-types produced above.
# Those are: Completeness, k-mer specatr, QV and CN

for (i in 1:length(unique(merqury_stats$sample))) {
    cur_sample <-  unique(merqury_stats$sample)[i]
    highest_val <-  merqury_stats %>%
        filter(sample == cur_sample) %>%
        filter(percent == max(percent)) %$%
        percent %>%
        unique()
    highest_stage  <-  merqury_stats %>%
        filter(sample == cur_sample) %>%
        filter(percent == highest_val) %$%
        stage %>%
        unique()
    lowest_val <-  merqury_stats %>%
        filter(sample == cur_sample) %>%
        filter(percent == min(percent)) %$%
        percent %>%
        unique()
    lowest_stage  <-  merqury_stats %>%
        filter(sample == cur_sample) %>%
        filter(percent == lowest_val) %$%
        stage %>%
        unique()
    highest_qv  <- merqury_qv  %>%
        filter(sample == cur_sample)  %>%
        filter(QV == max(QV)) %$%
        QV %>%
        unique()
    highest_qv_stage  <- merqury_qv  %>%
        filter(sample == cur_sample)  %>%
        filter(QV == max(QV)) %$%
        stage %>%
        unique()

    cat(paste('## ', cur_sample),
        paste0('\n\n'),
        paste0('### Valueboxes'),
        paste0('\n\n'),
        paste0('::: {.valuebox icon="exclude" color="primary" title="Merqury QV" }\n'),
        glue::glue("QV: {unique(highest_qv) %>% round(2)}, at stage(s): {paste(highest_qv_stage, collapse = ', ')}"),
        paste0('\n'),
        paste0(':::'),
        paste0('\n\n'),
        paste0('::: {.valuebox icon="percent" color="success" title="Highest k-mer completeness" }\n'),
        glue::glue("{unique(highest_val) %>% round(2)}%, at stage(s): {paste(highest_stage, collapse = ', ')}"),
        paste0('\n'),
        paste0(':::'),
        paste0('\n\n'),
        paste0('::: {.valuebox icon="heartbreak" color="warning" title="Lowest k-mer completeness" }\n'),
        glue::glue("{unique(lowest_val) %>% round(2)}%, at stage(s): {paste(lowest_stage, collapse = ', ')}"),
        paste0('\n'),
        paste0(':::'),
        paste0('\n\n'),
        paste0('\n\n'),
        paste0('### Plots { .tabset }'),
        paste0('\n\n'),
        paste0('#### Completeness \n'),
        paste0('\n'),
        knitr::knit_child(glue::glue('merqury_files/stat_plots/_{ cur_sample }_completeness_plt.Rmd'),
                    envir = globalenv(),
                    quiet = TRUE),
        paste0('\n'),
        paste0('#### QV \n'),
        paste0('\n'),
        paste0('QV is defined as:\n', expression(10*-log10(error_rate))),
        paste0('\n'),
        knitr::knit_child(glue::glue('merqury_files/qv_plots/_{ cur_sample }_qv_plt.Rmd'),
                    envir = globalenv(),
                    quiet = TRUE),
        paste0('\n'),
        paste0('#### Spectra \n'),
        paste0('\n'),
        knitr::knit_child(glue::glue('merqury_files/asm_plots/_{ cur_sample }_asm_plt.Rmd'),
                    envir = globalenv(),
                    quiet = TRUE),
        paste0('\n'),
        paste0('#### Copy Number \n'),
        paste0('\n'),
        knitr::knit_child(glue::glue('merqury_files/cn_plots/_{ cur_sample }_cn_plt.Rmd'),
                    envir = globalenv(),
                    quiet = TRUE),
        paste0('\n\n\n'),
        sep = "")
}
```
:::
:::

```{r}
#| eval: !expr params$merqury
# Delete files.
unlink("merqury_files/cn_plots")
unlink("merqury_files/asm_plots")
unlink("merqury_files")
```

# Software versions

The pipeline was run using the following software versions:

```{r}
versions <- yaml::read_yaml("software_versions.yml")
lapply(1:length(versions), \(process) {
    proc = versions[[process]]
    proc_name = names(versions[process])
    tools <- lapply(1:length(proc), \(tool) {
        tool_name = proc[tool] %>% names
        tool_version = proc[[tool]] %>% as.character()
        return(tibble(Process = proc_name,Tool = tool_name, Version = tool_version))
        }) %>%
    bind_rows()
    }) %>%
bind_rows() %>%
knitr::kable()
```
