parity_suite.sh: run on macOS, and flag runs that read no reads

CHANGELOG.md

@@ -13,6 +13,11 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Other changes
 
+- `test/parity_suite.sh` runs on macOS. It called GNU `time -f`, compared a
+  padded `wc -l` count, and passed `zcat`, which on macOS reads only `.Z`
+  files: both tools read 0 reads, exited 0, and every comparison came out
+  "identical". It now times with Python, uses `gzcat` where it exists, and
+  `parity_compare.py` warns when either tool read 0 reads.
 - Third-party notices completed. `LICENSES/` gains the licences of code
   ported from DegNorm (MIT) and David Elworthy's `sgt.h` via STAR (BSD
   2-clause), and of libc++ (Apache 2.0 with LLVM Exceptions), whose random

test/parity_compare.py

@@ -111,6 +111,10 @@ def align(a, b, name):
     diff = [f"{x.split('|')[0].strip()}: {x.split('|')[1].strip()} vs {y.split('|')[1].strip()}"
             for x, y in zip(la, lb) if x != y and "|" in x and "|" in y]
     print("   Log.final.out: " + ("identical" if la == lb else "; ".join(diff) or "differs"))
+    for tool, log in (("STAR", la), ("rustar-aligner", lb)):
+        n = next((x.split("|")[1].strip() for x in log if "Number of input reads" in x), "")
+        if n == "0":  # a --readFilesCommand that read nothing still exits 0
+            print(f"   WARNING: {tool} read 0 reads; check --readFilesCommand ({name})")
     sa = os.path.join(a, "Solo.out")
     if os.path.isdir(sa):
         sb = os.path.join(b, "Solo.out")

test/parity_suite.sh

@@ -5,7 +5,9 @@
 #   DATA=path/to/testdata STAR=path/to/STAR RUSTAR=target/release/rustar-aligner \
 #     test/parity_suite.sh fetch|index|run|compare|all
 #
-# Needs curl, seqtk, samtools and python3. THREADS (default 8) sets --runThreadN. Files already
+# Needs curl, seqtk, samtools and python3; runs on Linux and macOS. On macOS, Homebrew's STAR
+# stops on geneInfo.tab with --quantMode or --soloType, so use a STAR built for Linux there.
+# THREADS (default 8) sets --runThreadN. Files already
 # in $DATA are reused, so a slow download (the two ERR12389696 files are 1.6 GB each) can be
 # replaced by copies checked against ENA's MD5s.
 # Data (all public, pinned):
@@ -59,7 +61,7 @@ fetch() {
     local h=$P/chr21/ERR188021_2M_$m.fastq.gz
     if [ ! -s "$h" ]; then
       head_fq "https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR188/ERR188021/ERR188021_$m.fastq.gz" 8000000 > "$h.part"
-      [ "$(gunzip -c "$h.part" | wc -l)" = 8000000 ] && mv "$h.part" "$h" \
+      [ "$(gunzip -c "$h.part" | wc -l | tr -d " ")" = 8000000 ] && mv "$h.part" "$h" \
         || { echo "ERR188021_$m: download incomplete" >&2; exit 1; }
     fi
     get "$NFC/$MODULES/$g/test_starfusion_rnaseq_$m.fastq.gz" "$P/fusion/test_starfusion_rnaseq_$m.fastq.gz"
@@ -74,7 +76,18 @@ fetch() {
 
 TOOLS=(star rustar)
 bin() { [ "$1" = star ] && echo "$STAR" || echo "$RUSTAR"; }
-timed() { /usr/bin/time -f "%e %M" -o "$1" "${@:2}"; }
+# timed FILE CMD...: run CMD and write "wall-seconds peak-RSS-KB" to FILE, as GNU
+# `time -f "%e %M"` does; macOS's /usr/bin/time has no -f.
+TIMED='import resource, subprocess, sys, time
+t = time.time(); rc = subprocess.call(sys.argv[2:])
+rss = resource.getrusage(resource.RUSAGE_CHILDREN).ru_maxrss  # bytes on macOS, KB on Linux
+kb = rss // 1024 if sys.platform == "darwin" else rss
+open(sys.argv[1], "w").write("%.2f %d\n" % (time.time() - t, kb))
+sys.exit(rc)'
+timed() { python3 -c "$TIMED" "$@"; }
+# macOS's zcat only reads .Z files; gzcat is its gzip one. One word, as --readFilesCommand
+# takes one on main.
+ZCAT=$(command -v gzcat || echo zcat)
 
 index() {  # genome fasta gtf sjdbOverhang SAindexNbases [extra args]
   for t in "${TOOLS[@]}"; do
@@ -99,7 +112,7 @@ align() {  # name genome args...
   for t in "${TOOLS[@]}"; do
     local o=$P/out/${t}_$name/; mkdir -p "$o"
     timed "$o/time.txt" "$(bin $t)" --runThreadN "$THREADS" --genomeDir "$P/idx/${t}_$g" \
-      --readFilesCommand zcat --outFileNamePrefix "$o" "$@" > "$o/stdout.txt" 2>&1 \
+      --readFilesCommand "$ZCAT" --outFileNamePrefix "$o" "$@" > "$o/stdout.txt" 2>&1 \
       || echo "FAIL $t $name" >&2
   done
 }