[BenjaminDEMAILLE]

Part of #354 (lowercase and IUPAC letters in SEQ and Unmapped.out.mate*).

The other three items of the same audit (NH under --outSAMmultNmax, flag 0x200, --outFilterMismatchNoverReadLmax) are #375. This one is separate because it touches every place SEQ is written. The two PRs are independent; both add a test at the end of tests/alignment_features.rs and an entry under the same CHANGELOG heading, so the second to merge has a trivial conflict there.

Problem

rustar-aligner encodes the read to 0..4 (ACGTN) when it reads the FASTQ and builds SEQ from that encoding. So:

• lowercase bases came out uppercase (acgtACGT → ACGTACGT);
• IUPAC codes other than N came out as N (R → N), in SAM and in BAM;
• Unmapped.out.mate* had the same uppercase text.

While there: for pairs, Unmapped.out.mate* was written from the clipped mates, so --clip3pNbases 10 10 gave 90-base records where STAR writes all 100. Single-end already wrote the whole read.

What STAR does (2.7.11b source)

• ReadAlign_outputTranscriptSAM.cpp:194-204,232: SEQ is Read0[Mate], the read line as loaded, or on the reverse strand revComplementNucleotides(Read0[Mate], ...). Unmapped records write Read0[imate] as it is (:43).
• SequenceFuns.cpp:16-54: revComplementNucleotides complements A C G T N R Y K M S W B D V H in both cases, keeps the case, and leaves any other byte unchanged.
• ReadAlign_alignBAM.cpp:489-514, SequenceFuns.cpp:99-129: BAM packs the same text with nuclToNumBAM, which is case-blind and keeps the 16 IUPAC codes (anything else is N). noodles' BAM encoder does the same mapping, so putting the text in the record gives STAR's BAM.
• SequenceFuns.cpp:131-146: convertNucleotidesToNumbers makes the 0..4 copy used for mapping, lowercase as uppercase. rustar's encode_base already does this; mapping does not change.
• ReadAlign_outputAlignments.cpp:267,270: Unmapped.out.mate* gets Read0[im] and Qual0[im], the whole mate.

Change

• src/io/fastq.rs: EncodedRead.raw_letters keeps the read line, only when it has a byte outside ACGTN, so plain reads cost nothing. EncodedRead::letters() gives the text either way. revcomp_letters is STAR's table.
• src/io/sam.rs: restore_read_letters puts those letters in SEQ after the records are built: the mate is chosen by the LAST_SEGMENT flag, the text is reverse complemented for reverse-strand mapped records, and a leading hard clip is skipped. It replaces SEQ only where the record's bases are that text's bases in ACGTN form, so records from anything else are untouched.
• src/lib.rs: called next to the per-file read-group pass on the genomic and transcriptome records, in the single-end, paired-end and both STARsolo paths. UnmappedFastqWriter::write_record now takes the text; the four call sites pass letters() and, for pairs, the whole mate's quality.

Measured against STAR 2.7.11b

Synthetic reads on a synthetic genome, macOS, each tool with its own index. SEQ columns of Aligned.out.sam, samtools view of Aligned.out.bam, and Unmapped.out.mate1:

Read	STAR	rustar before	rustar after
SE forward, lowercase head and R	acgcag...TCTTRAGAAC...	ACGCAG...TCTTNAGAAC...	as STAR
SE reverse, lowercase, R and n	...GTGnGCTGACTGTYCGG...TGTgtttcg...	...GTGNGCTGACTGTNCGG...TGTGTTTCG...	as STAR
SE unmapped, all lowercase	lowercase (SAM and Unmapped.out.mate1)	uppercase	as STAR
PE mate 1 all lowercase / mate 2 with n	lowercase / n	uppercase / N	as STAR
BAM, reverse read with R	...CTGTYCGG...	...CTGTNCGG...	as STAR
PE unmapped, --clip3pNbases 10 10 --clip5pNbases 3 3, Unmapped.out.mate1 length	100, 100	87, 87	100, 100

Over the two synthetic sets (8 single-end reads, 5 pairs, plus the test's 3 reads and 2 pairs) every SEQ in SAM and BAM and every sequence line of Unmapped.out.mate* is identical to STAR's after the change. The Unmapped.out.mate* name lines still differ (0:Y: and the mate-mapped field); that is #364.

Tests

• tests/alignment_features.rs::test_seq_keeps_read_letters_like_star: single-end forward, reverse and unmapped reads and a mapped and an unmapped pair, with lowercase and IUPAC letters and --clip3pNbases, checking FLAG and SEQ of every record and Unmapped.out.mate1. I ran STAR on the same input; its output is what the test expects. It fails without the fix (ACGCAG...N... instead of acgcag...R...).
• src/io/fastq.rs::revcomp_letters_as_star: the whole revComplementNucleotides table, both cases, and two bytes outside it.
• src/io/sam.rs::restore_read_letters_after_hard_clip_and_not_on_other_bases.
• cargo test --release, cargo clippy --all-targets -- -D warnings and cargo fmt --check are clean on macOS.

Not checked

• Hard-clipped chimeric records (--chimOutType WithinBAM HardClip) against STAR: covered by the unit test only.
• Aligned.toTranscriptome.out.bam and the STARsolo paths: the code runs there, not compared with STAR.
• FASTA input and --readFilesType SAM: not run.
• The quality line of Unmapped.out.mate* under --outQSconversionAdd (STAR writes Qual0; not compared).
• No real-data harness; reads with lowercase or IUPAC letters are rare in real FASTQ, and for plain ones nothing changes (raw_letters stays None).

🤖 Generated with Claude Code