Port STAR's WithinBAM output (supersedes #282); PE nM is the pair's

[Psy-Fer]

--chimOutType WithinBAM ported from STAR; paired-end nM fix (supersedes #282)

WithinBAM output didn't match STAR for a single chimeric read: 0/20 SE and 0/12 PE on yeast 10k, compared record by record. The problems were structural:

• Duplicate output: the read's normal alignment was still written, plus two chimeric records.
• Incomplete records: the chimeric records had no qualities, no NH/HI/nM, wrong AS and no mate fields.
• Representative: the donor was always the representative.
• PE bracketing chimeras got a supplementary record.

#282 patched three of these on the old per-mate builder. This PR replaces that builder with a port of STAR's own code instead.

What changed

New src/chimeric/bam.rs. Ports ChimericAlign::chimericBAMoutput and the mapped branch of ReadAlign::alignBAM. It works from the same STAR-shaped window transcripts the detector (#284) uses, so a paired segment covering both mates becomes two records, as in STAR.

• Representative (chimRepresent / chimType):
  • a segment covering both mates is the representative;
  • mates on either side of the junction are both written as normal records, with no supplementary and no SA;
  • in SE, the higher-scoring segment is the representative.
• The other segment is supplementary (0x800): hard-clipped on the junction side (-11/-12, with SEQ/QUAL trimmed), or soft-clipped with --chimOutType ... SoftClip (-13).
• Each record gets STAR's flags (paired / proper / strand / mate-strand bits, 0x100 for non-best chimeras), MAPQ from the number of chimeras, mate fields per chimericBAMoutput, TLEN for two-mate records, and NH/HI/AS/nM/NM/XS/jM/jI/MD. AS and nM come from STAR's re-scoring of the placed transcripts.
• SA tags point the representative and supplementary at each other, using the other record's final CIGAR, MAPQ and NM.

The read's normal output is replaced, not supplemented. STAR returns before outputAlignments when it writes a chimera into the BAM (ReadAlign_oneRead.cpp:99), and so do we now:

• no normal records for the read;
• not counted as uniquely or multi-mapped in Log.final.out, but still counted as input. Unmapped reasons still count, since mappedFilter runs first;
• no SJ.out.tab, gene counts, transcriptome BAM or unmapped FASTQ for the read.

The four writer-thread loops that appended chimeric records after each batch are removed. Records are now emitted in read order from the alignment workers.

--chimOutType (ParametersChimeric_initialize.cpp):

• WithinBAM requires --outSAMtype BAM; with SAM output it is now an error, as in STAR.
• WithinBAM adds NM to the output attributes, as STAR does. The SA tag's last field reads it.
• Up to three values are accepted, e.g. Junctions WithinBAM HardClip. The last of HardClip/SoftClip wins, and hard clipping is the default.
• Unknown values are rejected. SeparateSAMold is refused rather than accepted and silently ignored; it is not implemented.

Paired-end nM (normal output too). STAR writes the pair's mismatch count (trOut.nMM) on both mates, the same way it writes AS. We wrote each mate's own count. That differed from STAR on 3,444 of 8,280 yeast PE pairs whose alignments otherwise match, and went unnoticed because compare_pe.py checks positions and CIGARs, not tags. After the fix, all 8,280 match STAR including every tag. The two unit tests that asserted the per-mate value were updated.

Results vs STAR 2.7.11b

Yeast 10k ($DATA), --chimSegmentMin 12 --outSAMtype BAM .... Every record of every chimeric read was compared on fields 2–11 plus all tags.

	before	after
SE WithinBAM (HardClip / SoftClip / sorted)	0 / 20	19 / 20
SE Junctions WithinBAM --chimMultimapNmax 10	—	40 / 42
PE WithinBAM (HardClip / SoftClip / sorted)	0 / 12	10 / 12
PE multimap	—	14 / 16
PE, all 84 chimeric reads incl. mates on either side of the junction	—	82 / 84

The remaining differences are known and outside this PR:

• SE: the one read whose 21-base partner our seeding doesn't produce.
• PE: two copy-choice ties (same records at the other copy of a repeat).

Log.final.out: SE now matches STAR's WithinBAM run exactly, including uniquely mapped (8249 vs 8265 without WithinBAM) and average mapped length. The only difference is the chimeric count, 19 vs 20, from the same SE miss. PE differs only in "splices: Total / Non-canonical" (273/22 vs 271/20), which is pre-existing: the same difference appears with no chimeric options at all.

Other modes:

• --outStd BAM_Unsorted gives the same records.
• 1 thread and 16 threads give identical records in the same order.
• WithinBAM with SAM output is refused.
• Without chimeric options, SE output is identical to main, and PE differs only in nM.

Tests

• bam.rs:
  • an SE chimera: representative by score, hard-clip side, trimmed SEQ, SA tags with the final CIGARs, and NM present;
  • mates on either side of the junction: two normal records, cross-pointing mate fields, no SA;
  • SoftClip mode.
• params: chim_out_type_follows_star, covering the BAM requirement, NM added, three values, last clip mode winning, and unknown / SeparateSAMold rejected.
• cargo fmt --check, cargo clippy --all-targets -- -D warnings (clippy 1.99) and cargo test (638 tests) all pass.

Docs

guides/chimeric.md still described the five detection tiers #284 removed, and two of its column descriptions were wrong (columns 11 and 13 are genomic start positions). It now describes STAR's two detectors, WithinBAM's real behaviour, and the p CIGAR. cli-parameters.md and output-files.md are updated to match.

R="$DATA/reads/ERR12389696_sub_1_10k.fastq.gz $DATA/reads/ERR12389696_sub_2_10k.fastq.gz"
STAR --genomeDir "$DATA/indices_star" --readFilesCommand zcat --readFilesIn $R \
  --chimSegmentMin 12 --chimOutType WithinBAM --outSAMtype BAM Unsorted --outFileNamePrefix "$DATA/wb_star_"
./target/release/rustar-aligner --runMode alignReads --genomeDir "$DATA/indices_rustar" \
  --readFilesCommand zcat --readFilesIn $R \
  --chimSegmentMin 12 --chimOutType WithinBAM --outSAMtype BAM Unsorted --outFileNamePrefix "$DATA/wb_rs_"
# compare per read: fields 2-11 plus sorted tags
for f in star rs; do samtools view "$DATA/wb_${f}_Aligned.out.bam" \
  | awk 'BEGIN{OFS="\t"}{t=""; for(i=12;i<=NF;i++) t=t" "$i; print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,t}' | sort > "/tmp/wb_$f.txt"; done
diff /tmp/wb_star.txt /tmp/wb_rs.txt | grep -c '^[<>]'

Supersedes #282.