Add benchmark runner script and regenerate all benchmark results

Add benchmark/run_benchmarks.py — a Python script that automates running
all upstream tools (dupRadar, featureCounts, 7 RSeQC tools) via Docker/Wave
containers and RustQC natively, with psrecord CPU/memory profiling.

Regenerated all benchmark results on large (10 GB BAM) and small datasets:
- RustQC (all tools, 10 threads): 16m 11s
- Traditional pipeline total: ~1h 25m (dupRadar 27m, featureCounts 4m,
  7 RSeQC tools 54m running sequentially)

Updated all markdown files (README.md, benchmark/README.md, docs benchmark
pages) with fresh timing numbers. Regenerated SVG bar charts with 10 bars
(RustQC + featureCounts + dupRadar + 7 individual RSeQC scripts).

New profiling output in benchmark/profiling/ with per-tool psrecord logs,
CPU/memory profile plots, and a results.json summary.

Author: ewels

README.md

@@ -19,7 +19,7 @@ All tools accept SAM/BAM/CRAM input and support processing multiple files in a s
 <p align="center">
 <picture>
    <source media="(prefers-color-scheme: dark)" srcset="docs/public/benchmarks/benchmark_dark.svg">
-   <img src="docs/public/benchmarks/benchmark_light.svg" alt="Benchmark: RustQC 54s vs featureCounts 16m vs dupRadar 30m" width="600">
+   <img src="docs/public/benchmarks/benchmark_light.svg" alt="Benchmark: RustQC 16m 11s vs traditional tools ~1h 25m (dupRadar + featureCounts + 7 RSeQC tools)" width="600">
 </picture>
 </p>
 
@@ -36,7 +36,7 @@ Performs dupRadar-equivalent duplicate rate analysis, featureCounts-compatible r
 | Language | R | Rust |
 | Dependencies | R, Bioconductor, Rsubread | None (static binary) |
 | Read counting | 4 separate featureCounts calls | Single-pass alignment reading |
-| Speed | ~30 min for 10 GB BAM | <1 min for 10 GB BAM |
+| Speed | ~30 min for 10 GB BAM | ~16 min for 10 GB BAM (all tools) |
 | Memory | High (R overhead) | Low |
 | Output format | Identical | Identical |
 

benchmark/README.md

@@ -8,9 +8,33 @@ Benchmark data and scripts for comparing RustQC against
 For detailed results, tables, and side-by-side plot comparisons, see the
 documentation:
 
+- [Combined benchmarks](https://ewels.github.io/RustQC/benchmarks/combined/)
 - [dupRadar benchmarks](https://ewels.github.io/RustQC/benchmarks/dupradar/)
 - [featureCounts benchmarks](https://ewels.github.io/RustQC/benchmarks/featurecounts/)
 
+## Latest results (large dataset)
+
+Run via `python3 benchmark/run_benchmarks.py --all` on a 10-core Apple Silicon
+Mac. Upstream tools run via Docker with x86 emulation; RustQC runs natively.
+
+| Tool | Wall time |
+|------|----------:|
+| dupRadar (R) | 27m 21s |
+| featureCounts (Subread) | 3m 39s |
+| bam_stat (RSeQC) | 6m 07s |
+| infer_experiment (RSeQC) | 7s |
+| read_duplication (RSeQC) | 29m 43s |
+| read_distribution (RSeQC) | 6m 00s |
+| junction_annotation (RSeQC) | 4m 37s |
+| junction_saturation (RSeQC) | 6m 32s |
+| inner_distance (RSeQC) | 1m 09s |
+| **Traditional total** | **1h 25m** |
+| **RustQC (10 threads)** | **16m 11s** |
+
+> **Note:** Docker x86 emulation on ARM inflates the upstream tool timings.
+> The key takeaway is that RustQC replaces 9 separate tool invocations
+> (each requiring a full BAM pass) with a single command and single BAM pass.
+
 ## Directory structure
 
 ```
@@ -86,7 +110,29 @@ pending):
 | `junction_saturation` | 256,466 junctions at 100% (163,710 known, 92,756 novel) |
 | `inner_distance` | 1,000,000 pairs sampled (mean 29.43, median 27.5, SD 32.80) |
 
-## Reproducing benchmarks
+## Running benchmarks
+
+The simplest way to run benchmarks is with the included script:
+
+```bash
+# Run everything (requires Docker for upstream tools)
+python3 benchmark/run_benchmarks.py --all
+
+# Run only specific tools
+python3 benchmark/run_benchmarks.py --rustqc --dupradar
+
+# Run only the small dataset
+python3 benchmark/run_benchmarks.py --all --small-only
+```
+
+The script profiles CPU/memory usage with `psrecord`, saves timing results to
+`benchmark/profiling/results.json`, per-tool profile logs and plots under
+`benchmark/profiling/<tool>/<dataset>/`, and regenerates the SVG bar charts
+in `docs/public/benchmarks/`.
+
+Dependencies: `pip install psrecord matplotlib psutil`
+
+## Manual reproduction
 
 ### 1. Download large benchmark input files
 

benchmark/RSeQC/large/infer_experiment.txt

@@ -1,4 +1,4 @@
-Reading reference gene model /data/genes.bed ... Done
+Reading reference gene model /data/input.bed ... Done
 Loading SAM/BAM file ...  Total 200000 usable reads were sampled
 
 

benchmark/RSeQC/large/inner_distance.inner_distance_plot.pdf

@@ -1,5 +1,5 @@
-out_file = '/out/inner_distance'
-pdf('/out/inner_distance.inner_distance_plot.pdf')
+out_file = '/data/output/inner_distance'
+pdf('/data/output/inner_distance.inner_distance_plot.pdf')
 fragsize=rep(c(-247.5,-242.5,-237.5,-232.5,-227.5,-222.5,-217.5,-212.5,-207.5,-202.5,-197.5,-192.5,-187.5,-182.5,-177.5,-172.5,-167.5,-162.5,-157.5,-152.5,-147.5,-142.5,-137.5,-132.5,-127.5,-122.5,-117.5,-112.5,-107.5,-102.5,-97.5,-92.5,-87.5,-82.5,-77.5,-72.5,-67.5,-62.5,-57.5,-52.5,-47.5,-42.5,-37.5,-32.5,-27.5,-22.5,-17.5,-12.5,-7.5,-2.5,2.5,7.5,12.5,17.5,22.5,27.5,32.5,37.5,42.5,47.5,52.5,57.5,62.5,67.5,72.5,77.5,82.5,87.5,92.5,97.5,102.5,107.5,112.5,117.5,122.5,127.5,132.5,137.5,142.5,147.5,152.5,157.5,162.5,167.5,172.5,177.5,182.5,187.5,192.5,197.5,202.5,207.5,212.5,217.5,222.5,227.5,232.5,237.5,242.5,247.5),times=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,72,398,420,458,532,567,737,957,1235,1537,1821,2075,2516,3238,4295,6187,9181,13629,19899,27200,37423,45241,55599,64240,70990,74193,75210,73148,67875,61125,53144,44885,36500,28801,22231,15885,10801,7035,4620,2968,2132,1513,1296,1154,1089,1018,1007,921,922,874,871,813,757,702,749,621,630,582,531,579,540,521,534,453,375,325,339,308,259,287,292))
 frag_sd = sd(fragsize)
 frag_mean = mean(fragsize)

benchmark/RSeQC/large/inner_distance.inner_distance_plot.r

@@ -0,0 +1,2 @@
+Get exon regions from /data/input.bed ...
+Load BAM file ...  Total read pairs  used 1000000

benchmark/RSeQC/large/inner_distance.txt

@@ -1,8 +1,8 @@
-pdf("/out/junction_annotation.splice_events.pdf")
+pdf("/data/output/junction_annotation.splice_events.pdf")
 events=c(1.1224763530979862,0.7011813022911578,98.012362937516)
 pie(events,col=c(2,3,4),init.angle=30,angle=c(60,120,150),density=c(70,70,70),main="splicing events",labels=c("partial_novel 1%","complete_novel 1%","known 98%"))
 dev.off()
-pdf("/out/junction_annotation.splice_junction.pdf")
+pdf("/data/output/junction_annotation.splice_junction.pdf")
 junction=c(19.86072227897655,10.42360390851029,69.71567381251316)
 pie(junction,col=c(2,3,4),init.angle=30,angle=c(60,120,150),density=c(70,70,70),main="splicing junctions",labels=c("partial_novel 20%","complete_novel 10%","known 70%"))
 dev.off()

benchmark/RSeQC/large/junction_annotation.junction_plot.r

@@ -0,0 +1,18 @@
+Reading reference bed file:  /data/input.bed  ...  Done
+Load BAM file ...  Done
+
+===================================================================
+Total splicing  Events:	13065665
+Known Splicing Events:	12805967
+Partial Novel Splicing Events:	146659
+Novel Splicing Events:	91614
+Filtered Splicing Events:	21425
+
+Total splicing  Junctions:	256466
+Known Splicing Junctions:	178797
+Partial Novel Splicing Junctions:	50936
+Novel Splicing Junctions:	26733
+
+===================================================================
+Create BED file ...
+Create Interact file ...