From eeed9e2bdf35e47bac415628c602b06c08b21325 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 27 Aug 2026 21:30:40 +0200 Subject: [PATCH] test: yeast-tier differential in a minute, and a SAM agreement tool CONTRIBUTING.md's differential wants the whole ERR12389696 run, which is why "measure it on the yeast tier first" has been easy to skip. This fetches a byte range instead: the first ~120 MB of each mate is 10 000 pairs, the reads are ordered so the mates still pair up, and the whole thing takes about a minute including both index builds. test/sam_agreement.py is the comparison itself, usable on any two SAM files. It reports position agreement and NH agreement separately on purpose: a tie broken differently moves the first and not the second, while a real regression usually moves both, so one number cannot tell them apart. Measured on main at 10 000 pairs: 7860 uniquely mapped against STAR's 7861, 98.48% of mates at the same chromosome, position and CIGAR, 99.96% with the same NH. Written while investigating #31, where it is what showed that a candidate fix regressed faithfulness rather than improving it. --- CONTRIBUTING.md | 18 +++++++++++ test/sam_agreement.py | 46 +++++++++++++++++++++++++++++ test/yeast_tier.sh | 69 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 133 insertions(+) create mode 100644 test/sam_agreement.py create mode 100755 test/yeast_tier.sh diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 551505ff..5da8c872 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -63,6 +63,24 @@ 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. +## Yeast tier in a minute + +The full differential in the next section wants the whole ERR12389696 run. +`test/yeast_tier.sh` fetches the first ~120 MB of each mate instead — 10 000 +pairs, enough to see a faithfulness regression in the numbers that matter — +builds both indexes, runs both aligners and prints the per-read agreement. + +```bash +cargo build --release +test/yeast_tier.sh # ~1 minute, needs STAR on PATH +DATA=~/yeast-tier PAIRS=50000 test/yeast_tier.sh +``` + +`test/sam_agreement.py` is the comparison on its own, if you already have two +SAM files. It prints position agreement and NH agreement separately, because a +tie broken differently shows up in the first and not the second, while a real +regression usually moves both. + ## 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/sam_agreement.py b/test/sam_agreement.py new file mode 100644 index 00000000..ebfce0f9 --- /dev/null +++ b/test/sam_agreement.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python3 +"""Per-read agreement between two SAM files, primary records only. + + python3 test/sam_agreement.py star_Aligned.out.sam rustar_Aligned.out.sam + +Reports how many mates land at the same chromosome, position and CIGAR, and +how many carry the same NH. Ties broken differently by the two aligners show +up as position differences with equal NH, which is why the two figures are +printed apart: a drop in NH agreement is a different problem from a drop in +position agreement. +""" +import sys +def load(path): + out={} + for line in open(path): + if line.startswith("@"): continue + f=line.rstrip("\n").split("\t") + flag=int(f[1]) + if flag & 0x100 or flag & 0x800: continue + mate = 2 if flag & 0x80 else 1 + key=(f[0], mate) + nh=0 + for t in f[11:]: + if t.startswith("NH:i:"): nh=int(t[5:]) + out[key]=(f[2], f[3], f[5], nh, flag & 0x4) + return out +a,b=load(sys.argv[1]),load(sys.argv[2]) +keys=set(a)|set(b) +same=pos_same=nh_same=0 +only_a=only_b=0 +nh_hist_a={}; nh_hist_b={} +for k in keys: + x,y=a.get(k),b.get(k) + if x is None: only_b+=1; continue + if y is None: only_a+=1; continue + if x[4]==0: nh_hist_a[x[3]]=nh_hist_a.get(x[3],0)+1 + if y[4]==0: nh_hist_b[y[3]]=nh_hist_b.get(y[3],0)+1 + if x[:3]==y[:3]: pos_same+=1 + if x[:4]==y[:4]: same+=1 + if x[3]==y[3]: nh_same+=1 +n=len(keys) +print(f"records compared: {n} only in A: {only_a} only in B: {only_b}") +print(f"same chr/pos/CIGAR : {pos_same} ({pos_same/n:.4%})") +print(f"same incl. NH : {same} ({same/n:.4%})") +print(f"same NH : {nh_same} ({nh_same/n:.4%})") +print("max NH A:", max(nh_hist_a or {0:0}), " B:", max(nh_hist_b or {0:0})) diff --git a/test/yeast_tier.sh b/test/yeast_tier.sh new file mode 100755 index 00000000..628065cb --- /dev/null +++ b/test/yeast_tier.sh @@ -0,0 +1,69 @@ +#!/usr/bin/env bash +# Per-read agreement against STAR on the project's yeast tier. +# +# The full benchmark in CONTRIBUTING.md wants the whole ERR12389696 run; this +# fetches the first ~120 MB of each mate instead, which is 10 000 pairs and +# enough to see a faithfulness regression in the numbers that matter (position +# agreement, unique/multi counts, NH depth). It runs in about a minute. +# +# test/yeast_tier.sh # build, run, compare +# DATA=~/yeast-tier test/yeast_tier.sh # keep the downloads somewhere +# +# Requires STAR 2.7.11b on PATH. +set -euo pipefail + +DATA="${DATA:-/tmp/rustar-yeast-tier}" +RUSTAR="${RUSTAR:-./target/release/rustar-aligner}" +PAIRS="${PAIRS:-10000}" +THREADS="${THREADS:-4}" + +mkdir -p "$DATA" +cd "$DATA" + +if [ ! -f genome.fa ]; then + echo "fetching the yeast reference" + curl -sfL "https://ftp.ensembl.org/pub/release-110/fasta/saccharomyces_cerevisiae/dna/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz" -o genome.fa.gz + gunzip -kf genome.fa.gz +fi + +lines=$((PAIRS * 4)) +for m in 1 2; do + if [ ! -f "r${m}.fq" ]; then + echo "fetching mate ${m} (partial)" + # A byte range rather than the whole run: the reads are ordered, so the + # first N records of each mate still pair up. + curl -sfL -r 0-120000000 \ + "ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR123/096/ERR12389696/ERR12389696_${m}.fastq.gz" \ + -o "r${m}_part.fastq.gz" + gunzip -c "r${m}_part.fastq.gz" 2>/dev/null | head -"${lines}" > "r${m}.fq" || true + fi +done + +for pair in "rustar:$OLDPWD/$RUSTAR:ridx" "star:STAR:sidx"; do + name="${pair%%:*}"; rest="${pair#*:}"; exe="${rest%%:*}"; idx="${rest##*:}" + if [ ! -f "${idx}/SA" ]; then + echo "building the ${name} index" + mkdir -p "$idx" + "$exe" --runMode genomeGenerate --genomeDir "$idx" \ + --genomeFastaFiles genome.fa --genomeSAindexNbases 11 \ + --outFileNamePrefix "${name}_idx_" > /dev/null + fi +done + +echo "aligning with rustar-aligner" +"$OLDPWD/$RUSTAR" --runMode alignReads --genomeDir ridx --readFilesIn r1.fq r2.fq \ + --outSAMtype SAM --runThreadN "$THREADS" --outFileNamePrefix rustar_ > /dev/null + +echo "aligning with STAR" +STAR --genomeDir sidx --readFilesIn r1.fq r2.fq \ + --outSAMtype SAM --runThreadN "$THREADS" --outFileNamePrefix star_ > /dev/null + +echo +for f in star rustar; do + printf '%-7s ' "$f" + grep -E "Uniquely mapped reads number|mapped to multiple loci \|" "${f}_Log.final.out" \ + | sed 's/^ *//' | tr '\n' ' ' + echo +done +echo +python3 "$OLDPWD/test/sam_agreement.py" star_Aligned.out.sam rustar_Aligned.out.sam