---
title: "Benchmarking BamScale: reader throughput and end-to-end workflow impact"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Benchmarking BamScale: reader throughput and end-to-end workflow impact}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE,
                      fig.width = 7, fig.height = 4, dpi = 120)
has_ggplot <- requireNamespace("ggplot2", quietly = TRUE)
```

## Overview

`BamScale` is a multithreaded BAM reader, built on the `ompBAM` OpenMP engine, that
returns the same native Bioconductor objects as the standard readers and is a
drop-in replacement for `Rsamtools::scanBam` and
`GenomicAlignments::readGAlignments`.

Reading a BAM is only one step of an analysis, so a faster reader speeds a workflow
end-to-end **only in proportion to the time that workflow spends reading**
(Amdahl's law). This vignette summarises benchmarks that quantify *both* the raw
read speedup and its end-to-end effect on two standard workflows — coverage/bigWig
track generation and ATAC-seq fragment-size QC — with BamScale output verified
**byte-identical** to the standard tools.

The numbers below are a concise summary of a representative benchmark run
(Intel Xeon Gold 6252, 96 cores; warm page cache; median of 5 iterations). The
benchmark harness that produces them is in `inst/benchmarks/`
(`run_server_benchmark.R`, `run_workflow_benchmark.R`); see the last section to
reproduce the full results.

```{r results-data}
## Summary results (from the benchmark runs; see inst/benchmarks/ to reproduce).
read_tbl <- data.frame(
  Workload    = c("Core fields", "GAlignments", "Sequence + quality"),
  Reads       = c("qname/flag/rname/pos/mapq/cigar", "-> GAlignments object", "seq + base quality"),
  `Standard (s)` = c(15.6, 10.9, 22.3),
  `BamScale (s)` = c(6.8, 3.4, 8.1),
  Threads     = c(48L, 48L, 24L),
  Speedup     = c(2.29, 3.22, 2.76),
  check.names = FALSE
)

spectrum_tbl <- data.frame(
  Endpoint      = c("ATAC fragment-size QC", "GAlignments (read)", "Sequence + quality (read)",
                    "Coverage -> RleList", "Core fields (read)", "Coverage -> bigWig"),
  Layer         = c("workflow", "read", "read", "workflow", "read", "workflow"),
  `Read fraction` = c(0.93, 1.00, 1.00, 0.76, 1.00, 0.21),
  Speedup       = c(3.66, 3.22, 2.76, 2.50, 2.29, 1.23),
  check.names = FALSE
)
```

## Read throughput on a single large BAM

On a single large BAM, BamScale reads **2.3--3.2x faster** than the standard reader
across three representative access patterns, because it spreads the BGZF decode
across cores:

```{r read-table}
knitr::kable(
  transform(read_tbl, Speedup = sprintf("%.2fx", Speedup)),
  caption = "Single-file read throughput: best BamScale configuration vs the single-threaded standard reader (median of 5 iterations)."
)
```

```{r read-plot, eval=has_ggplot}
library(ggplot2)
ggplot(read_tbl, aes(stats::reorder(Workload, Speedup), Speedup)) +
  geom_col(fill = "#2a78d6", width = 0.65) +
  geom_hline(yintercept = 1, linetype = "dashed", colour = "grey55") +
  geom_text(aes(label = sprintf("%.2fx", Speedup)), hjust = -0.15, size = 3.5) +
  coord_flip() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(x = NULL, y = "speedup over standard reader (dashed = parity)",
       title = "Single-file read speedup") +
  theme_minimal(base_size = 12)
```

Two honest points. At **one thread** BamScale is at or near the single-threaded
readers (the win comes from threading, not a faster single-core path), and the
scaling is **strongly sublinear** -- roughly 2--3x on 48 threads, saturating by
about 24 threads. The value is not near-linear scaling but the ability to use cores
that a single-threaded reader cannot use at all -- as on one large BAM in an
interactive session. Object-free command-line decoders such as `samtools` reach
several-fold higher *raw* throughput, but return no Bioconductor object; BamScale's
role is fast, object-faithful decoding *inside* R.

## End-to-end impact tracks the read fraction

Embedding BamScale as the read step of a complete workflow, the end-to-end speedup
follows how read-dominated that workflow is (Amdahl's law):

```{r spectrum-table}
knitr::kable(
  transform(spectrum_tbl,
            `Read fraction` = paste0(round(100 * spectrum_tbl$`Read fraction`), "%"),
            Speedup = sprintf("%.2fx", Speedup)),
  caption = "End-to-end speedup by workload endpoint. Read-pattern rows are single-file read throughput; workflow rows are end-to-end. Read fraction is the share of the standard pipeline spent reading."
)
```

```{r spectrum-plot, eval=has_ggplot, fig.height=4.4}
ggplot(spectrum_tbl, aes(stats::reorder(Endpoint, Speedup), Speedup, fill = `Read fraction`)) +
  geom_col(width = 0.65) +
  geom_hline(yintercept = 1, linetype = "dashed", colour = "grey55") +
  geom_text(aes(label = sprintf("%.2fx", Speedup)), hjust = -0.15, size = 3.4) +
  coord_flip() +
  scale_fill_gradient(low = "#cfe0f3", high = "#12355b", labels = scales::percent) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(x = NULL, y = "end-to-end speedup (dashed = parity)", fill = "read\nfraction",
       title = "The gain tracks the read fraction") +
  theme_minimal(base_size = 12)
```

- **ATAC fragment-size QC** is nearly pure read, so BamScale's read speedup carries
  through almost intact: **3.66x** end-to-end.
- **Coverage -> bigWig** is dominated by the shared, single-threaded bigWig writer
  (identical work on both arms), so the end-to-end gain is bounded at **1.23x** even
  though the read phase itself is much faster. Stopping at the in-memory coverage
  `RleList` -- what many analyses consume next -- recovers a **2.50x** gain.

## A note on cores and fairness

A multithreaded reader can occupy cores a single-threaded one cannot. At **matched
core counts** (both arms given the same number of cores) the two readers are
approximately at parity for multi-file processing (coverage ~1.10x, ATAC ~0.96x):
the single-threaded reader already saturates a machine by parallelising *across*
files, so BamScale's multi-file advantage comes specifically from also threading
*within* a file -- which pays off when there are fewer files than cores.

## Output is identical to the standard tools

Every speedup above is against output verified equal to the standard tool. The
coverage `RleList` is `identical()` to `GenomicAlignments::coverage()` on
`readGAlignments` output at every thread count, and the ATAC fragment-size table is
`identical()` to `Rsamtools::scanBam` and byte-identical (maximum absolute
difference 0) to `ATACseqQC::fragSizeDist` over 49.8 million reads. This is what
makes BamScale a genuine drop-in rather than an approximation.

## Try it

```{r demo, echo=TRUE, eval=TRUE}
library(BamScale)
bam <- ompBAM::example_BAM("Unsorted")

## Alignment fields as a GAlignments object (drop-in for readGAlignments)
ga <- bam_read(bam, what = c("rname", "pos", "cigar", "strand"),
               as = "GAlignments", threads = 2)
ga

## ...or core fields as a data.frame; `threads` controls within-file parallelism
df <- bam_read(bam, what = c("qname", "flag", "mapq"),
               as = "data.frame", threads = 2)
head(df)
```

## Reproducing the full benchmarks

The complete benchmark harness ships under `inst/benchmarks/`:

```{r repro, echo=TRUE, eval=FALSE}
dir(system.file("benchmarks", package = "BamScale"))
# run_server_benchmark.R      : read-pattern micro-benchmarks (step1 / GAlignments / seq+qual)
# run_workflow_benchmark.R    : end-to-end coverage and ATAC fragment-size QC workflows
# download_atac_data.R        : fetch + index the ENCODE GM12878 ATAC BAMs used here
```

Each writes machine-readable result tables (`summary.csv`) plus host and
correctness metadata, from which the summary figures above are derived.

## Session information

```{r session-info}
sessionInfo()
```