test: yeast-tier differential in a minute, and a SAM agreement tool

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.

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}))

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