Skip to content

dev(compare-asvs): --per-sample mode with within-sample flip pairing #384

dev(compare-asvs): --per-sample mode with within-sample flip pairing

dev(compare-asvs): --per-sample mode with within-sample flip pairing #384

Workflow file for this run

name: Concordance (vs R DADA2)
# Sanity guardrail: run short dada2-rs pipelines on small fixtures and compare
# the post-chimera ASVs + per-ASV counts against a STATIC R DADA2 reference CSV
# (committed under dev/concordance/reference/; CI never runs R).
#
# GATING: the compare steps run with --gate, so a metric below its threshold
# fails the build. Thresholds are deliberately generous (the calibrated baseline
# is recall/precision = 1.000, count_corr ~0.994 on both platforms) so the gate
# catches real regressions without flaking on cross-implementation noise.
#
# Triggers: pull requests, pushes to non-main branches, and manual dispatch
# (the intended way to run it against main).
on:
pull_request:
push:
branches-ignore: [main]
workflow_dispatch:
concurrency:
group: concordance-${{ github.ref }}
cancel-in-progress: true
env:
CARGO_TERM_COLOR: always
# Gating thresholds (generous vs the 1.000/1.000/~0.994 baseline → catch real
# regressions, not cross-implementation noise). Tune here.
MIN_RECALL: "0.85"
MIN_PRECISION: "0.80"
MIN_COUNT_CORR: "0.95"
MIN_ABUNDANCE: "2"
# Pseudo arm only. Measured on the committed fixture: recall 0.974, precision
# 0.993 against the R pseudo reference, the gap being R's between-round
# error-model re-fit (issue #100), which we deliberately do not do. Set below
# measured so the known gap passes but a real regression does not.
PSEUDO_MIN_RECALL: "0.93"
PSEUDO_MIN_PRECISION: "0.95"
jobs:
concordance:
name: Concordance vs R reference
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v5
- name: Install Rust stable
uses: dtolnay/rust-toolchain@stable
- name: Cache cargo
uses: actions/cache@v5
with:
path: |
~/.cargo/registry/index
~/.cargo/registry/cache
~/.cargo/git/db
target
key: ${{ runner.os }}-cargo-concordance-${{ hashFiles('**/Cargo.lock') }}
restore-keys: ${{ runner.os }}-cargo-
# Built with `--features wfa` because the WFA concordance steps below drive
# the experimental backend via `ALIGN_BACKEND=wfa2`. The feature is
# off-by-default (issue #63), so without it those steps would hit the
# NW-only stub. The resulting binary still defaults to NW for the plain
# concordance runs.
- name: Build (release, with WFA backend)
run: cargo build --release --features wfa
# ---------------- Illumina (in-repo MiSeq tutorial pair) ----------------
- name: Illumina — run dada2-rs pipeline
run: |
bash dev/concordance/run_illumina.sh \
./target/release/dada2-rs \
dev/concordance/data/illumina work/illumina 2
- name: Illumina — compare to R reference
run: |
ref=dev/concordance/reference/illumina_seqtab_nochim.csv
if [ ! -f "$ref" ]; then
echo "::notice::Illumina reference CSV not committed yet ($ref); skipping comparison. Generate with write_reference.R."
echo "Illumina reference not yet committed — pipeline ran but no comparison." >> "$GITHUB_STEP_SUMMARY"
exit 0
fi
python3 dev/concordance/compare_to_reference.py \
--rs work/illumina/seqtab.nochim.json --reference "$ref" \
--label illumina --min-abundance "$MIN_ABUNDANCE" \
--min-recall "$MIN_RECALL" --min-precision "$MIN_PRECISION" \
--min-count-corr "$MIN_COUNT_CORR" --gate \
--summary "$GITHUB_STEP_SUMMARY" --json work/illumina_report.json
# -------- Illumina with the experimental WFA backend (issue #49) ---------
# Same fixtures + same R reference + same gate, but the whole pipeline
# (learn-errors, dada, remove-bimera-denovo) runs with --align-backend
# wfa2. WFA is ASV-equivalent to NW on this data, so it must clear the
# same thresholds; this guards against the WFA path regressing.
#
# CAUTION — this gate is REASSURING BUT NOT SUFFICIENT. WFA's ends-free
# alignment is genuinely suboptimal (upstream WFA2-lib #102; see
# nwalign.rs::wfa_endsfree_known_divergence). On these small fixtures the
# divergence is below the count thresholds so the gate passes, but at
# scale WFA over-calls low-abundance ASVs (jaccard ~0.92-0.95 unbounded,
# ~0.997-0.999 capped) — see issue #51 scale validation. A green WFA gate
# means "did not regress on the fixture", NOT "WFA == NW". WFA stays
# EXPERIMENTAL; NW is the default + error-model backend.
- name: Illumina (WFA) — run dada2-rs pipeline
run: |
ALIGN_BACKEND=wfa2 bash dev/concordance/run_illumina.sh \
./target/release/dada2-rs \
dev/concordance/data/illumina work/illumina_wfa 2
- name: Illumina (WFA) — compare to R reference
run: |
ref=dev/concordance/reference/illumina_seqtab_nochim.csv
if [ ! -f "$ref" ]; then
echo "::notice::Illumina reference CSV not committed yet ($ref); skipping WFA comparison."
exit 0
fi
echo "::warning title=WFA gate is fixture-scale only::A green Illumina WFA gate means the experimental WFA backend did not regress on this small fixture, NOT that WFA == NW. WFA ends-free is suboptimal (WFA2-lib #102) and over-calls low-abundance ASVs at scale (issue #51). NW remains the default + error-model backend."
python3 dev/concordance/compare_to_reference.py \
--rs work/illumina_wfa/seqtab.nochim.json --reference "$ref" \
--label illumina-wfa --min-abundance "$MIN_ABUNDANCE" \
--min-recall "$MIN_RECALL" --min-precision "$MIN_PRECISION" \
--min-count-corr "$MIN_COUNT_CORR" --gate \
--summary "$GITHUB_STEP_SUMMARY" --json work/illumina_wfa_report.json
# -------- Illumina, pseudo-pooling (issue #100) --------------------------
# dada-pseudo vs an R `dada(pool="pseudo")` reference. Separate fixture:
# the 2-sample illumina fixture selects priors that change nothing, so
# pseudo there is indistinguishable from a per-sample run and the check
# would pass while testing nothing. run_illumina.sh POOL=pseudo enforces
# this with a positive control that fails if the priors turn out inert.
#
# A KNOWN, UNDERSTOOD GAP is expected here. R re-fits the error model
# between its two pseudo rounds; we deliberately do not (see
# docs/findings/pseudo-pooling-priors-vs-error-model.md). On this fixture
# that costs recall 0.974 / precision 0.993 -- 4 low-abundance ASVs R calls
# and we do not. Running our opt-in emulation
# (--reestimate-err-between-rounds) instead scores 1.000/1.000 against the
# same reference, which is how we know the residual is entirely the re-fit
# and not drift. Thresholds are set below the measured default-arm values,
# so a regression in the pseudo path still trips the gate.
- name: Illumina pseudo — run pipeline
run: |
POOL=pseudo bash dev/concordance/run_illumina.sh \
./target/release/dada2-rs \
dev/concordance/data/illumina_pseudo work/illumina_pseudo 2
- name: Illumina pseudo — compare to R reference
run: |
ref=dev/concordance/reference/illumina_pseudo_seqtab_nochim.csv
if [ ! -f "$ref" ]; then
echo "::notice::Illumina pseudo reference CSV not committed yet ($ref); skipping comparison. Generate with: Rscript dev/concordance/write_reference.R illumina dev/concordance/data/illumina_pseudo \"$ref\" --pool=pseudo"
echo "Illumina pseudo reference not yet committed — pipeline ran but no comparison." >> "$GITHUB_STEP_SUMMARY"
exit 0
fi
python3 dev/concordance/compare_to_reference.py \
--rs work/illumina_pseudo/seqtab.nochim.json --reference "$ref" \
--label illumina-pseudo --min-abundance "$MIN_ABUNDANCE" \
--min-recall "$PSEUDO_MIN_RECALL" \
--min-precision "$PSEUDO_MIN_PRECISION" \
--min-count-corr "$MIN_COUNT_CORR" --gate \
--summary "$GITHUB_STEP_SUMMARY" --json work/illumina_pseudo_report.json
# ---------------- PacBio (committed subsampled fixture) -----------------
# Skipped until both the fixture and the reference CSV are committed.
- name: PacBio — run + compare (if fixture present)
run: |
data=dev/concordance/data/pacbio
ref=dev/concordance/reference/pacbio_seqtab_nochim.csv
if [ ! -d "$data" ] || [ -z "$(ls -A "$data"/*.fastq.gz 2>/dev/null)" ]; then
echo "::notice::PacBio fixture not committed yet ($data); skipping."
echo "PacBio fixture not yet committed — skipped." >> "$GITHUB_STEP_SUMMARY"
exit 0
fi
# k=5 matches R's fixed KMER_SIZE (apples-to-apples vs the reference).
bash dev/concordance/run_pacbio.sh \
./target/release/dada2-rs "$data" work/pacbio \
AGRGTTYGATYMTGGCTCAG RGYTACCTTGTTACGACTT 2
# k=7 is the recommended PacBio setting; the k-mer screen is a prefilter,
# so it should match the same R(k=5) reference. Asserts the production
# config stays concordant.
PACBIO_KMER=7 bash dev/concordance/run_pacbio.sh \
./target/release/dada2-rs "$data" work/pacbio_k7 \
AGRGTTYGATYMTGGCTCAG RGYTACCTTGTTACGACTT 2
if [ ! -f "$ref" ]; then
echo "::notice::PacBio reference CSV not committed yet ($ref); pipeline ran, no comparison."
echo "PacBio reference not yet committed — pipeline ran but no comparison." >> "$GITHUB_STEP_SUMMARY"
exit 0
fi
python3 dev/concordance/compare_to_reference.py \
--rs work/pacbio/seqtab.nochim.json --reference "$ref" \
--label pacbio-k5 --min-abundance "$MIN_ABUNDANCE" \
--min-recall "$MIN_RECALL" --min-precision "$MIN_PRECISION" \
--min-count-corr "$MIN_COUNT_CORR" --gate \
--summary "$GITHUB_STEP_SUMMARY" --json work/pacbio_k5_report.json
python3 dev/concordance/compare_to_reference.py \
--rs work/pacbio_k7/seqtab.nochim.json --reference "$ref" \
--label pacbio-k7 --min-abundance "$MIN_ABUNDANCE" \
--min-recall "$MIN_RECALL" --min-precision "$MIN_PRECISION" \
--min-count-corr "$MIN_COUNT_CORR" --gate \
--summary "$GITHUB_STEP_SUMMARY" --json work/pacbio_k7_report.json
# -------- PacBio with the experimental WFA backend (issue #51) ----------
# WFA on PacBio was originally ~2.3x slower than NW at k=5 (the looser
# screen lets divergent non-error-copy pairs through, and WFA pays its full
# O(n·s) extend on them). The default WFA edit-budget cap (--wfa-max-edits,
# issue #51) bounds those pairs — capped pairs fall back to NW — which
# erases the slowdown (locally faster than NW) while staying ASV-identical
# to the R reference (93=93). k=5 matches R's KMER_SIZE. WFA is
# experimental; this gate guards the WFA+cap path against regressing.
#
# CAUTION — same caveat as the Illumina WFA gate: PASSING HERE IS
# FIXTURE-SCALE ONLY. The 93=93 match holds because the committed fixture
# is small; at 95-sample scale capped WFA is NW parity-to-slower (NOT the
# local speedup) and unbounded WFA is a 3.6x slowdown that over-calls ~160
# ASVs (issue #51 scale validation). The cap does not make WFA == NW — it
# bounds the divergence. Green = "did not regress on the fixture".
- name: PacBio (WFA) — run + compare (if fixture present)
run: |
data=dev/concordance/data/pacbio
ref=dev/concordance/reference/pacbio_seqtab_nochim.csv
if [ ! -d "$data" ] || [ -z "$(ls -A "$data"/*.fastq.gz 2>/dev/null)" ]; then
echo "::notice::PacBio fixture not committed yet ($data); skipping WFA gate."
exit 0
fi
ALIGN_BACKEND=wfa2 bash dev/concordance/run_pacbio.sh \
./target/release/dada2-rs "$data" work/pacbio_wfa \
AGRGTTYGATYMTGGCTCAG RGYTACCTTGTTACGACTT 2
if [ ! -f "$ref" ]; then
echo "::notice::PacBio reference CSV not committed yet ($ref); WFA pipeline ran, no comparison."
exit 0
fi
echo "::warning title=WFA gate is fixture-scale only::A green PacBio WFA gate means the experimental WFA+cap backend did not regress on this small fixture, NOT that WFA == NW. At scale capped WFA is NW parity-to-slower and unbounded WFA over-calls ~160 ASVs (issue #51). NW remains the default + error-model backend."
python3 dev/concordance/compare_to_reference.py \
--rs work/pacbio_wfa/seqtab.nochim.json --reference "$ref" \
--label pacbio-k5-wfa --min-abundance "$MIN_ABUNDANCE" \
--min-recall "$MIN_RECALL" --min-precision "$MIN_PRECISION" \
--min-count-corr "$MIN_COUNT_CORR" --gate \
--summary "$GITHUB_STEP_SUMMARY" --json work/pacbio_wfa_report.json
- name: Upload reports
if: always()
uses: actions/upload-artifact@v4
with:
name: concordance-reports
path: work/*_report.json
if-no-files-found: ignore