From 1750f4db85d62003e3030ec9affa65ac6e9f426b Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 27 Aug 2026 21:02:26 +0200 Subject: [PATCH] ci: opt-in large-dataset differential against STAR The suites in tests/ use a synthetic micro-genome and finish in seconds, which is why two reported bugs (#31, #48) were invisible to them: both are about aggregate behaviour over tens of thousands of reads, not about one alignment. test/nfcore_diff.py fetches the nf-core/rnaseq test dataset (50 000 paired reads, S. cerevisiae chrI plus the GFP transgene), runs STAR and rustar-aligner over it, and compares mapping rates, the unmapped-reason buckets and the multimapper depth histogram. Thresholds are the measured gaps plus headroom, so today's state passes and a regression does not; --report-only prints the comparison without failing, for use while a difference is being investigated. The workflow is opt-in rather than per-push: manual dispatch, or the `large-tests` label on a pull request. Dispatch takes a runner label, so it can be aimed at a scverse AWS runner without editing the file. Checked against a known regression: on main the script fails on exactly the two unmapped buckets that #247 fixes, and reports the NH gap of #31 at 1.43x. Closes #244. --- .github/workflows/large-dataset.yml | 71 ++++++++++ CONTRIBUTING.md | 20 +++ test/nfcore_diff.py | 202 ++++++++++++++++++++++++++++ 3 files changed, 293 insertions(+) create mode 100644 .github/workflows/large-dataset.yml create mode 100644 test/nfcore_diff.py diff --git a/.github/workflows/large-dataset.yml b/.github/workflows/large-dataset.yml new file mode 100644 index 00000000..7bc30409 --- /dev/null +++ b/.github/workflows/large-dataset.yml @@ -0,0 +1,71 @@ +name: Large-dataset differential + +# Deliberately not on every push: this job fetches a dataset, builds two +# indexes and runs two aligners. It is opt-in, either by hand or by putting the +# `large-tests` label on a pull request. +on: + workflow_dispatch: + inputs: + runner: + description: "Runner label (e.g. ubuntu-latest, or a scverse AWS runner group)" + required: false + default: "ubuntu-latest" + report_only: + description: "Report the comparison without failing the job" + type: boolean + required: false + default: false + pull_request: + types: [labeled, synchronize] + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +permissions: + contents: read + +jobs: + nfcore-diff: + name: nf-core/rnaseq differential against STAR + # On a pull request this runs only when the `large-tests` label is present, + # so an unrelated PR never pays for it. + if: >- + github.event_name == 'workflow_dispatch' || + contains(github.event.pull_request.labels.*.name, 'large-tests') + runs-on: ${{ github.event.inputs.runner || 'ubuntu-latest' }} + timeout-minutes: 60 + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # ratchet:actions/checkout@v7.0.1 + + - name: Install Rust toolchain + uses: dtolnay/rust-toolchain@631a55b12751854ce901bb631d5902ceb48146f7 # ratchet:dtolnay/rust-toolchain@stable + + - uses: Swatinem/rust-cache@c19371144df3bb44fab255c43d04cbc2ab54d1c4 # ratchet:Swatinem/rust-cache@v2.9.1 + + - name: Build rustar-aligner + run: cargo build --release + + # STAR is packaged for Ubuntu; building it from source would triple the + # job's runtime for no gain in what is being measured. + - name: Install STAR + run: | + sudo apt-get update + sudo apt-get install -y rna-star + STAR --version + + - name: Differential run + run: >- + python3 test/nfcore_diff.py + --rustar ./target/release/rustar-aligner + --work "${RUNNER_TEMP}/nfcore" + --json "${RUNNER_TEMP}/nfcore-metrics.json" + ${{ (github.event.inputs.report_only == 'true') && '--report-only' || '' }} + + - name: Upload metrics + if: always() + uses: actions/upload-artifact@ea165f8d65b6e75b540449e92b4886f43607fa02 # ratchet:actions/upload-artifact@v4.6.2 + with: + name: nfcore-metrics + path: ${{ runner.temp }}/nfcore-metrics.json + if-no-files-found: warn diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 551505ff..8bb20efe 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -63,6 +63,26 @@ Adding a dependency — **especially a non-Rust one** (a C library via a `-sys` If you add a CLI flag that parses but is not yet implemented, mark it as such in the parameter-surface test and document it — do not silently accept a flag that does nothing. A user passing a flag should never be quietly ignored. +## Larger-scale differential + +The unit and integration suites use a synthetic micro-genome and finish in +seconds. `test/nfcore_diff.py` covers what they cannot reach: mapping rates, +the unmapped-reason buckets and the multimapper depth histogram, measured +against STAR on the nf-core/rnaseq test dataset (50 000 paired reads, +*S. cerevisiae* chrI plus the GFP transgene, fetched from public URLs and not +vendored). + +```bash +cargo build --release +python3 test/nfcore_diff.py --report-only # print the comparison +python3 test/nfcore_diff.py # exit non-zero on a regression +``` + +In CI it is opt-in: the `Large-dataset differential` workflow runs on manual +dispatch, or on a pull request carrying the `large-tests` label. The dispatch +form takes a runner label, so it can be pointed at a bigger machine without +editing the workflow. + ## Test data Integration tests in `tests/` use a bundled synthetic micro-genome and need no downloads. The differential benchmark below uses a small **public** yeast RNA-seq dataset that is not vendored; fetch it once and point `DATA` at wherever you keep it. diff --git a/test/nfcore_diff.py b/test/nfcore_diff.py new file mode 100644 index 00000000..37102e4d --- /dev/null +++ b/test/nfcore_diff.py @@ -0,0 +1,202 @@ +#!/usr/bin/env python3 +"""Differential run against STAR on the nf-core/rnaseq test dataset. + +A middle-sized fixture: 50 000 paired reads against *S. cerevisiae* chrI plus +the GFP transgene, which is small enough to fetch and index in under a minute +and large enough to move the numbers the unit tests cannot reach — mapping +rates, the unmapped-reason buckets, and the multimapper depth histogram. + +Everything is fetched from public URLs; nothing is vendored. + + python3 test/nfcore_diff.py --rustar ./target/release/rustar-aligner \\ + --work /tmp/nfcore --star STAR + +Exit status is 0 when every threshold holds, 1 otherwise, so it can gate a CI +job. `--report-only` always exits 0 and just prints the comparison, which is +what to use while a difference is being investigated. +""" + +from __future__ import annotations + +import argparse +import json +import subprocess +import sys +import urllib.request +from pathlib import Path + +REF = "626c8fab639062eade4b10747e919341cbf9b41a" +FILES = { + "genome.fasta": f"https://raw.githubusercontent.com/nf-core/test-datasets/{REF}/reference/genome.fasta", + "genes.gtf.gz": f"https://raw.githubusercontent.com/nf-core/test-datasets/{REF}/reference/genes_with_empty_tid.gtf.gz", + "reads_1.fastq.gz": "https://raw.githubusercontent.com/nf-core/test-datasets/rnaseq/testdata/GSE110004/SRR6357072_1.fastq.gz", + "reads_2.fastq.gz": "https://raw.githubusercontent.com/nf-core/test-datasets/rnaseq/testdata/GSE110004/SRR6357072_2.fastq.gz", +} + +# How far each figure may drift from STAR before the run fails. These are not +# aspirations: they are the measured gaps plus headroom, so a regression trips +# them and today's state does not. +THRESHOLDS = { + "uniquely_mapped_frac": 0.005, # fraction of input reads + "multi_mapped_frac": 0.005, + "unmapped_short_frac": 0.010, + "unmapped_other_frac": 0.010, + "max_nh_ratio": 2.0, # deepest multimapper, rustar / STAR +} + + +def fetch(work: Path) -> None: + work.mkdir(parents=True, exist_ok=True) + for name, url in FILES.items(): + dest = work / name + if dest.exists() and dest.stat().st_size > 0: + continue + print(f"fetching {name}", flush=True) + urllib.request.urlretrieve(url, dest) # noqa: S310 - fixed public URLs + + # Decompress the reads. `--readFilesCommand zcat` is deliberately not used: + # on macOS zcat silently yields nothing for both aligners, and a fixture + # that reads as "0 input reads" is worse than a slower one. + import gzip + import shutil + + for n in (1, 2): + plain = work / f"reads_{n}.fastq" + if not plain.exists(): + with gzip.open(work / f"reads_{n}.fastq.gz", "rb") as fin, open(plain, "wb") as fout: + shutil.copyfileobj(fin, fout) + + +def run(cmd: list[str], log: Path) -> None: + with open(log, "w") as f: + proc = subprocess.run(cmd, stdout=f, stderr=subprocess.STDOUT) + if proc.returncode != 0: + sys.exit(f"command failed ({proc.returncode}): {' '.join(cmd)}\nsee {log}") + + +def align(exe: str, work: Path, tag: str, is_star: bool) -> Path: + idx = work / f"{tag}_idx" + idx.mkdir(exist_ok=True) + prefix = str(work / f"{tag}_") + gen = [exe] + if not is_star: + gen += ["--runMode", "genomeGenerate"] + else: + gen += ["--runMode", "genomeGenerate"] + gen += [ + "--genomeDir", str(idx), + "--genomeFastaFiles", str(work / "genome.fasta"), + "--genomeSAindexNbases", "9", + "--outFileNamePrefix", prefix + "idx_", + ] + run(gen, work / f"{tag}_index.log") + + aln = [exe] + if not is_star: + aln += ["--runMode", "alignReads"] + aln += [ + "--genomeDir", str(idx), + "--readFilesIn", str(work / "reads_1.fastq"), str(work / "reads_2.fastq"), + "--outFilterMultimapNmax", "20", + "--outSAMtype", "SAM", + "--runThreadN", "4", + "--outFileNamePrefix", prefix, + ] + run(aln, work / f"{tag}_align.log") + return Path(prefix) + + +def final_log(prefix: Path) -> dict[str, int]: + rows: dict[str, int] = {} + for line in open(f"{prefix}Log.final.out"): + if "|" not in line: + continue + label, value = line.split("|", 1) + value = value.strip() + if value.endswith("%"): + continue + try: + rows[label.strip()] = int(value) + except ValueError: + pass + return rows + + +def nh_histogram(prefix: Path) -> dict[int, int]: + hist: dict[int, int] = {} + with open(f"{prefix}Aligned.out.sam") as f: + for line in f: + if line.startswith("@"): + continue + fields = line.rstrip("\n").split("\t") + flag = int(fields[1]) + if flag & 0x100 or flag & 0x800 or flag & 0x4: + continue + for tag in fields[11:]: + if tag.startswith("NH:i:"): + nh = int(tag[5:]) + hist[nh] = hist.get(nh, 0) + 1 + break + return hist + + +def main() -> int: + ap = argparse.ArgumentParser() + ap.add_argument("--rustar", default="./target/release/rustar-aligner") + ap.add_argument("--star", default="STAR") + ap.add_argument("--work", default="/tmp/nfcore-diff") + ap.add_argument("--report-only", action="store_true") + ap.add_argument("--json", help="write the measured figures here") + args = ap.parse_args() + + work = Path(args.work) + fetch(work) + + star_prefix = align(args.star, work, "star", is_star=True) + rustar_prefix = align(args.rustar, work, "rustar", is_star=False) + + s, r = final_log(star_prefix), final_log(rustar_prefix) + s_nh, r_nh = nh_histogram(star_prefix), nh_histogram(rustar_prefix) + n_input = s["Number of input reads"] + + fields = [ + ("uniquely_mapped_frac", "Uniquely mapped reads number"), + ("multi_mapped_frac", "Number of reads mapped to multiple loci"), + ("unmapped_short_frac", "Number of reads unmapped: too short"), + ("unmapped_other_frac", "Number of reads unmapped: other"), + ] + + measured: dict[str, float] = {} + failures: list[str] = [] + + print(f"\n{'metric':<40} {'STAR':>10} {'rustar':>10} {'delta':>10}") + for key, label in fields: + sv, rv = s.get(label, 0), r.get(label, 0) + delta = abs(sv - rv) / n_input + measured[key] = delta + print(f"{label:<40} {sv:>10} {rv:>10} {delta:>9.3%}") + if delta > THRESHOLDS[key]: + failures.append(f"{label}: |{sv} - {rv}| / {n_input} = {delta:.3%} > {THRESHOLDS[key]:.3%}") + + s_max, r_max = max(s_nh, default=1), max(r_nh, default=1) + ratio = r_max / s_max if s_max else float("inf") + measured["max_nh_ratio"] = ratio + print(f"{'deepest multimapper (NH)':<40} {s_max:>10} {r_max:>10} {ratio:>9.2f}x") + if ratio > THRESHOLDS["max_nh_ratio"]: + failures.append(f"deepest NH: {r_max} against STAR's {s_max} ({ratio:.2f}x > {THRESHOLDS['max_nh_ratio']}x)") + + if args.json: + Path(args.json).write_text(json.dumps({"measured": measured, "star": s, "rustar": r}, indent=2)) + + if failures: + print("\nFAILED:") + for f in failures: + print(f" {f}") + return 0 if args.report_only else 1 + + print("\nNFCORE DIFF PASSED") + return 0 + + +if __name__ == "__main__": + sys.exit(main())