---
title: "Benchmarking BamScaleR: reader throughput and end-to-end workflow impact"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Benchmarking BamScaleR: 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

`BamScaleR` 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 BamScaleR output verified
**byte-identical** to the standard tools.

The numbers below are a concise summary of a representative benchmark run
(a deeply sequenced public ENCODE ATAC-seq BAM, ENCFF191LNL, 151,057,932
alignments; Intel Xeon Gold 6252, 96 logical 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 (benchmark run 2026-10-08, BamScaleR 0.99.1, on the 25-file
## public ENCODE ATAC-seq corpus; 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(263.3, 290.1, 409.8),
  `BamScaleR (s)` = c(136.0, 115.9, 179.3),
  Threads     = c(36L, 72L, 24L),
  Speedup     = c(1.94, 2.50, 2.29),
  check.names = FALSE
)

spectrum_tbl <- data.frame(
  Endpoint      = c("ATAC fragment-size QC", "GAlignments (read)", "Coverage -> RleList",
                    "Sequence + quality (read)", "Core fields (read)", "Coverage -> bigWig"),
  Layer         = c("workflow", "read", "workflow", "read", "read", "workflow"),
  `Read fraction` = c(0.94, 1.00, 0.72, 1.00, 1.00, 0.41),
  Speedup       = c(4.39, 2.50, 2.49, 2.29, 1.94, 1.52),
  check.names = FALSE
)
```

## Read throughput on a single large BAM

On a single large BAM, BamScaleR reads **1.9--2.5x faster** than the standard reader
across three representative access patterns, because it spreads the BGZF decode
across cores and assembles the output objects in native code:

```{r read-table}
knitr::kable(
  transform(read_tbl, Speedup = sprintf("%.2fx", Speedup)),
  caption = "Single-file read throughput: best BamScaleR 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** BamScaleR is modestly faster than the
single-threaded readers (about 1.1--1.3x, from native-code object assembly); most
of the win comes from threading, and the scaling is **strongly sublinear** --
roughly 2.5--4x on 48 threads, largely saturating between 24 and 48 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; BamScaleR's role is
fast, object-faithful decoding *inside* R.

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

Embedding BamScaleR 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 BamScaleR's read speedup carries
  through almost intact: **4.39x** end-to-end.
- **Coverage -> bigWig** is dominated by the shared, single-threaded
  `rtracklayer::export.bw` step (identical work on both arms), so the end-to-end
  gain is bounded at **1.52x** even though the read phase itself is 5.9x faster.
  Stopping at the in-memory coverage `RleList` -- what many analyses consume
  next -- recovers a **2.49x** gain. (`bam_coverage()` and
  `bam_coverage_bigwig()` compute these endpoints directly inside the reader,
  removing the R-side read phase entirely; see the function documentation.)

## 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.20x, ATAC ~0.97x):
the single-threaded reader already saturates a machine by parallelising *across*
files, so BamScaleR's multi-file advantage comes specifically from also threading
*within* a file -- which pays off when there are fewer files than cores.

## Beyond reading: in-reader aggregation

For workflows whose endpoint is a summary rather than the alignments themselves,
BamScaleR provides four functions that fold the computation inside the
multithreaded reader, so no per-read R objects are ever materialised:
`fragment_sizes()` (paired-end fragment-size distribution),
`mapq_dist()` (mapping-quality distribution), `bam_coverage()` (per-base coverage
`RleList`), and `bam_coverage_bigwig()` (single-pass BAM to bigWig with parallel
block compression). Each is verified byte-identical to its standard equivalent
(`table(abs(scanBam(isize)))`, `table(scanBam(mapq))`,
`coverage(readGAlignments())`, and the serial bigWig writer respectively) -- see
their reference pages. These remove the read phase from the Amdahl decomposition
above entirely; comprehensive benchmarks on public ENCODE data accompany the
manuscript.

## 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 all 151,057,932 alignments. This is what
makes BamScaleR a genuine drop-in rather than an approximation.

## Try it

```{r demo, echo=TRUE, eval=TRUE}
library(BamScaleR)
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 = "BamScaleR"))
# 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()
```