[BenjaminDEMAILLE]

Part of #354 (NH under --outSAMmultNmax, SAM flag 0x200, --outFilterMismatchNoverReadLmax).

Three record and parameter differences from STAR, each small. Lowercase bases in SEQ, the fourth item of the same audit, is #376, separate because it touches every place SEQ is written.

Problem

1. NH was the number of records written, not the number of loci. A single-end read with 3 loci under --outSAMmultNmax 2 had NH:i:2; STAR writes NH:i:3. Paired-end had the same bug.
2. SAM flag 0x200 (not passing quality control) was never set. STAR sets it on every record of a read whose FASTQ name has the Illumina filter flag, e.g. @r1 1:Y:0:ACGT: 512 when mapped, 516 when unmapped.
3. --outFilterMismatchNoverReadLmax was refused by the CLI (tests/parameter_surface.rs listed it as not yet accepted). STAR 2.7.11b accepts it.

What STAR does (2.7.11b source)

NH:

• ReadAlign_outputAlignments.cpp:168: nTrOutWrite=min(P.outSAMmultNmax,nTrOutSAM) is only the number of records written.
• ReadAlign_outputAlignments.cpp:177,184: each record gets nTrOutSAM, all loci, as nTrOut, which is written as NH (ReadAlign_outputTranscriptSAM.cpp:281, ReadAlign_alignBAM.cpp:293). Same code for single-end and paired-end.

0x200:

• ReadAlignChunk_processChunks.cpp:121-127: readFilter is Y when the second word of the read-name line has :Y: at bytes 1 to 3. Only mate 1's name line is read; both mates get the same value (:133-140). readLoad.cpp:27 reads it back.
• ReadAlign_outputTranscriptSAM.cpp:27-28 (unmapped) and :91 (mapped), ReadAlign_alignBAM.cpp:146,198: 0x200 when readFilter=='Y'. The transcriptome BAM goes through alignBAM too.
• ReadAlign_outputTranscriptSAM.cpp:219, ReadAlign_alignBAM.cpp:543: --outSAMflagAND/OR are applied after 0x200, on mapped records only.

--outFilterMismatchNoverReadLmax:

• parametersDefault:473: default 1.0. Parameters.cpp:152: a double.
• ReadAlign_oneRead.cpp:78: outFilterMismatchNmaxTotal=min(outFilterMismatchNmax, (uint)(outFilterMismatchNoverReadLmax*(readLength[0]+readLength[1]))), per read, after clipping.
• That total replaces outFilterMismatchNmax everywhere the aligner checks mismatches: stitching (stitchAlignToTranscript.cpp:310,314,374,391), extension (ReadAlign_stitchWindowSeeds.cpp:58,83,146,165, stitchWindowAligns.cpp:45,66), the mapped filter (ReadAlign_mappedFilter.cpp:12) and the transcriptome projection (ReadAlign_quantTranscriptome.cpp:57).
• ReadAlign_waspMap.cpp:120 (copyRead): the merged read of --peOverlapNbasesMin keeps the pair's total.

Change

• src/io/sam.rs: transcript_to_record gets n_alignments (all loci) instead of max_output (records written), in write_alignment, build_alignment_records and build_paired_records.
• src/io/fastq.rs: EncodedRead.filter_flag, parsed from the FASTQ description by filter_flag(). src/lib.rs: apply_qc_fail sets 0x200 on the read's genomic and transcriptome records, from mate 1 for pairs, and applies --outSAMflagAND/OR again on mapped records so the order is STAR's. It runs next to the per-file read-group pass, in the single-end, paired-end and both STARsolo paths.
• src/params/mod.rs: the new parameter and Parameters::mismatch_nmax_total(read_len). src/align/read_align.rs sets the per-read scorer's n_mm_max from it (stitching and extension already read n_mm_max) and uses it in both mapped filters; the merged-mates read gets the pair's total. src/quant/transcriptome.rs used the formula with the ratio fixed at 1; it now calls the same function.

Overlap with #364. #364 (fix/unmapped-mate-names) adds the same filter_flag() function and EncodedRead.filter_flag field, for the Unmapped.out.mate* names. This PR copies them as they are there, including the unit test, so whichever merges second has a conflict in src/io/fastq.rs and the src/solo/mod.rs test fixtures that resolves by keeping one copy. This PR does not touch the Unmapped.out.mate* names.

Measured against STAR 2.7.11b

Synthetic genome (3 x 20 kb, one 400 bp segment on each chromosome), macOS, each tool with its own index.

Records of the same reads, --outSAMmultNmax 2 --outSAMunmapped Within:

Read	STAR	rustar before	rustar after
SE, 3 loci	FLAG 0 / 256, NH:i:3	0 / 256, NH:i:2	0 / 256, NH:i:3
PE, 3 loci	99 / 147 / 355 / 403, NH:i:3	same flags, NH:i:2	same flags, NH:i:3
SE 1:Y: mapped	512	0	512
SE 1:Y: unmapped	516	4	516
PE 1:Y: mapped	611 / 659	99 / 147	611 / 659
PE 1:Y: unmapped	589 / 653	77 / 141	589 / 653

The same in BAM (samtools view): identical to STAR after, apart from the IUPAC letter, which is #376. Aligned.toTranscriptome.out.bam of the mapped 1:Y: read: FLAG 512 after (rustar only; STAR's quant mode does not run on this macOS build).

--outFilterMismatchNoverReadLmax, 4 000 single-end reads and 4 000 pairs with 0 to 9 random mismatches per mate. Before, rustar refused the option (exit 2). After:

Run	STAR records / unique / too many MM	rustar records / unique / too many MM	records differing
SE, 0.04	3 323 / 3 002 / 0	3 323 / 3 002 / 0	0
PE, 0.04	8 422 / 3 754 / 0	8 422 / 3 754 / 0	0
SE, 1.0	4 153 / 3 911 / 0	4 153 / 3 911 / 0	0
PE, 1.0	8 244 / 3 898 / 0	8 244 / 3 898 / 0	0
SE, 0.04, EndToEnd	2 092 / 1 982 / 1 980	2 095 / 1 987 / 1 969	7
PE, 0.04, EndToEnd	5 442 / 2 678 / 1 304	5 998 / 2 936 / 995	2 158

The EndToEnd gaps are not from this option. In both tools, ratio 0.04 gives exactly the same records as --outFilterMismatchNmax 4 (SE) and --outFilterMismatchNmax 8 (PE), and at ratio 1.0 EndToEnd PE already differs (6 810 vs 7 202 records). That is a separate, existing difference under --alignEndsType EndToEnd for pairs.

Tests

• tests/alignment_features.rs::test_nh_qc_fail_and_mismatch_over_read_length_match_star: maps single-end and paired-end reads on a genome with a repeat planted three times and checks FLAG and NH of every record against STAR 2.7.11b's output on the same input (I ran STAR on it). Without the fix it fails: first because the CLI refuses the option, and with only the NH and 0x200 parts reverted, on NH:i:2 and the missing 0x200.
• src/io/sam.rs::test_out_sam_mult_nmax asserted the old NH (3 of 5 loci); it now expects 5.
• src/io/fastq.rs::filter_flag_needs_y_as_the_second_colon_field_of_a_one_byte_first_field (the same as in Unmapped.out.mate* names and sequences as STAR writes them #364).
• tests/parameter_surface.rs: the option leaves NOT_YET_ACCEPTED.
• cargo test --release, cargo clippy --all-targets -- -D warnings and cargo fmt --check are clean on macOS.

Not checked

• STARsolo paths and --outSAMflagAND/OR together with a 1:Y: read: the code path is there, not compared with STAR.
• The transcriptome BAM against STAR (STAR's quant mode does not run on the macOS build here).
• A half-mapped pair still writes NH:i:1 for its mapped mate whatever the number of loci. That is a separate path (build_half_mapped_records), not changed here.
• Unmapped.out.mate* names (0:Y:) are Unmapped.out.mate* names and sequences as STAR writes them #364's.
• No real-data harness (parity_suite.sh); the cases above are synthetic.

🤖 Generated with Claude Code