NH over all loci, 0x200 from the Illumina filter flag, --outFilterMismatchNoverReadLmax

CHANGELOG.md

@@ -67,6 +67,14 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Features
 
+- `--outFilterMismatchNoverReadLmax` is accepted (default 1.0), as in STAR
+  2.7.11b. It caps the mismatches of a read at this fraction of its length
+  (both mates, after clipping), together with `--outFilterMismatchNmax`, in
+  the same places STAR does: while stitching and extending, in the mapped
+  filter and for `--quantMode TranscriptomeSAM`. On 4 000 reads and 4 000
+  pairs with 0 to 9 mismatches each, at 0.04, the SAM records are identical
+  to STAR's.
+
 - **`--clipAdapterType CellRanger4` now matches STAR exactly**, both halves of the
   clip (#198). The 5' TSO trim is an overlap alignment against the first 91 bases of the
   read, replicating `ClipCR4`'s Opal call (ACGTN alphabet; match +1, mismatch -2,
@@ -232,6 +240,14 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Bug fixes
 
+- `NH` counts every locus of a multimapper, also when `--outSAMmultNmax`
+  writes fewer. A read with 3 loci under `--outSAMmultNmax 2` had `NH:i:2`,
+  STAR writes `NH:i:3`. Single-end and paired-end.
+- A read whose FASTQ name has the Illumina filter flag (`@r1 1:Y:0:ACGT`) gets
+  SAM flag 0x200 (not passing quality control) on all its records, mapped or
+  unmapped, and on both mates (from mate 1's name), in SAM, BAM and
+  `Aligned.toTranscriptome.out.bam`, as STAR does. It was never set.
+
 - **STARsolo's stats files as STAR writes them.** `Solo.out/Barcodes.stats`,
   `<feature>/Features.stats` and `<feature>/UMIperCellSorted.txt` were not
   written, SJ and Velocyto had no `Summary.csv`, and the Gene/GeneFull

docs/src/content/docs/reference/cli-parameters.md

@@ -65,6 +65,7 @@ Run `rustar-aligner --help` for the full machine-generated listing.
 | `--outFilterMultimapScoreRange` | `1` | Score range for keeping multi-mappers within best score. |
 | `--outFilterMismatchNmax` | `10` | Max mismatches per pair. |
 | `--outFilterMismatchNoverLmax` | `0.3` | Max ratio of mismatches to mapped length. |
+| `--outFilterMismatchNoverReadLmax` | `1.0` | Max ratio of mismatches to read length (both mates). |
 | `--outFilterScoreMin` | `0` | Min absolute alignment score. |
 | `--outFilterScoreMinOverLread` | `0.66` | Min alignment score normalized to read length. |
 | `--outFilterMatchNmin` | `0` | Min absolute matched bases. |

src/align/read_align.rs

@@ -409,7 +409,8 @@ fn align_read_inner(
     // does NOT filter clusters by seed coverage. Removed to match STAR.
 
     // Step 3: Stitch seeds within each cluster
-    let scorer = AlignmentScorer::from_params(params);
+    let mut scorer = AlignmentScorer::from_params(params);
+    scorer.n_mm_max = params.mismatch_nmax_total(read_seq.len());
     let mut transcripts = Vec::new();
     // STAR's `trAll`: every window's transcripts, before any cross-window
     // dedup or filtering, which is what its chimeric detection reads.
@@ -577,7 +578,7 @@ fn align_read_inner(
                 || n_match < (params.out_filter_match_nmin_over_lread * lread_m1) as u32
             {
                 Some(UnmappedReason::TooShort)
-            } else if best.n_mismatch > params.out_filter_mismatch_nmax
+            } else if best.n_mismatch > scorer.n_mm_max
                 || f64::from(best.n_mismatch) / f64::from(r_length.max(1))
                     > params.out_filter_mismatch_nover_lmax
             {
@@ -709,7 +710,8 @@ pub fn align_paired_read(
     let len2 = mate2_seq.len();
 
     let debug_pe = !params.read_name_filter.is_empty() && read_name == params.read_name_filter;
-    let scorer = AlignmentScorer::from_params(params);
+    let mut scorer = AlignmentScorer::from_params(params);
+    scorer.n_mm_max = params.mismatch_nmax_total(len1 + len2);
     let junction_db = if index.junction_db.is_empty() {
         None
     } else {
@@ -1081,6 +1083,9 @@ pub fn align_paired_read(
             merge_params.out_filter_match_nmin = 0;
             merge_params.out_filter_match_nmin_over_lread = 0.0;
             merge_params.out_filter_score_min_over_lread = 0.0;
+            // STAR's copyRead hands the merged read the pair's mismatch cap.
+            merge_params.out_filter_mismatch_nmax = scorer.n_mm_max;
+            merge_params.out_filter_mismatch_nover_read_lmax = 1.0;
             let (merged_transcripts, _merged_chim, _n_mapq, _unmapped, _best) =
                 align_read_inner(&merge.merged, read_name, index, &merge_params, false)?;
             let mut converted: Vec<PairedAlignment> = Vec::new();
@@ -1559,7 +1564,7 @@ fn filter_paired_transcripts(
             return Some(UnmappedReason::TooShort);
         }
 
-        if combined_nm > params.out_filter_mismatch_nmax
+        if combined_nm > params.mismatch_nmax_total(mate1_seq.len() + mate2_seq.len())
             || f64::from(combined_nm) / f64::from(combined_r_length.max(1))
                 > params.out_filter_mismatch_nover_lmax
         {

src/io/fastq.rs

@@ -45,6 +45,20 @@ impl UnmappedFastqWriter {
     }
 }
 
+/// `Y` when the first word of a FASTQ description has `:Y:` from its second byte,
+/// STAR's Illumina filter flag (`ReadAlignChunk_processChunks.cpp:121-127`).
+fn filter_flag(description: &[u8]) -> u8 {
+    let field = description
+        .split(u8::is_ascii_whitespace)
+        .next()
+        .unwrap_or(&[]);
+    if field.get(1..4) == Some(b":Y:") {
+        b'Y'
+    } else {
+        b'N'
+    }
+}
+
 /// A read from a FASTQ file with encoded bases
 #[derive(Debug, Clone)]
 pub struct EncodedRead {
@@ -58,6 +72,8 @@ pub struct EncodedRead {
     /// Index of the `--readFilesIn` file this read came from (STAR
     /// `readFilesIndex`); selects the read group.
     pub file_index: u32,
+    /// `Y` when the second field of a FASTQ name line reads `n:Y:...`, else `N`.
+    pub filter_flag: u8,
 }
 
 /// A paired-end read from two FASTQ files
@@ -235,6 +251,7 @@ impl FastqReader {
                     sequence,
                     quality,
                     file_index: self.file_index,
+                    filter_flag: filter_flag(record.description()),
                 }))
             }
             Some(Err(e)) => Err(Error::from(e)),
@@ -496,6 +513,26 @@ mod tests {
         assert_eq!(encode_base(b'S'), 4);
     }
 
+    /// Expected flags are STAR's on the same name lines.
+    #[test]
+    fn filter_flag_needs_y_as_the_second_colon_field_of_a_one_byte_first_field() {
+        for (description, flag) in [
+            (&b"1:Y:0:ACGT"[..], b'Y'),
+            (b"2:Y:18:GG", b'Y'),
+            (b"1:Y:0 2:N:0", b'Y'),
+            (b"1:N:0:ACGT", b'N'),
+            (b"Y", b'N'),
+            (b":Y:", b'N'),
+            (b"1:Y", b'N'),
+            (b"ab:Y:", b'N'),
+            (b"1:y:0", b'N'),
+            (b"x\ty:Y:z", b'N'),
+            (b"", b'N'),
+        ] {
+            assert_eq!(filter_flag(description), flag, "{description:?}");
+        }
+    }
+
     #[test]
     fn test_decode_base() {
         assert_eq!(decode_base(0), b'A');

src/io/sam.rs

@@ -213,7 +213,7 @@ impl SamWriter {
                 read_qual,
                 genome,
                 mapq,
-                max_output,    // NH = number of reported alignments
+                n_alignments, // NH counts every locus, not only the written ones (STAR nTr)
                 hit_index + 1, // 1-based
                 params.out_sam_attr_ih_start,
                 attrs,
@@ -363,7 +363,7 @@ impl SamWriter {
                 aligned_qual,
                 genome,
                 mapq,
-                max_output,    // NH = number of reported alignments
+                n_alignments, // NH counts every locus, not only the written ones (STAR nTr)
                 hit_index + 1, // 1-based
                 params.out_sam_attr_ih_start,
                 attrs,
@@ -482,7 +482,7 @@ impl SamWriter {
                     paired_aln.insert_size,
                     params,
                 ),
-                max_output, // NH = number of reported alignments
+                n_alignments, // NH counts every locus, not only the written ones (STAR nTr)
                 hit_index,
                 params.out_sam_attr_ih_start,
                 combined_score,
@@ -513,7 +513,7 @@ impl SamWriter {
                     -paired_aln.insert_size, // Negative for mate2 (mode 1 default)
                     params,
                 ),
-                max_output, // NH = number of reported alignments
+                n_alignments, // NH counts every locus, not only the written ones (STAR nTr)
                 hit_index,
                 params.out_sam_attr_ih_start,
                 combined_score,
@@ -3861,13 +3861,13 @@ mod tests {
         // Only 3 records output despite 5 transcripts
         assert_eq!(records.len(), 3);
 
-        // NH should be 3 (number of reported alignments)
+        // NH counts all 5 loci, as STAR's does (`outputTranscriptSAM(.., nTrOutSAM, ..)`)
         for rec in &records {
             let data = rec.data();
             assert_eq!(
                 data.get(&Tag::ALIGNMENT_HIT_COUNT),
-                Some(&Value::from(3_i32)),
-                "NH should be 3 (capped by outSAMmultNmax)"
+                Some(&Value::from(5_i32)),
+                "NH is not capped by outSAMmultNmax"
             );
         }
 

src/lib.rs

@@ -833,6 +833,51 @@ fn apply_rg_per_file(
     }
 }
 
+/// Flag every record of a read whose FASTQ name has the Illumina filter flag
+/// `n:Y:` as not passing quality control (0x200), mapped or not, as STAR does
+/// (`ReadAlign_outputTranscriptSAM.cpp:27,91`, `ReadAlign_alignBAM.cpp:146,198`).
+/// STAR applies `--outSAMflagAND/OR` after setting the bit, on mapped records
+/// only, so they are applied again here.
+fn apply_qc_fail(
+    records: &mut [noodles::sam::alignment::record_buf::RecordBuf],
+    filter_flag: u8,
+    params: &Parameters,
+    and_or: bool,
+) {
+    if filter_flag != b'Y' {
+        return;
+    }
+    for rec in records {
+        rec.flags_mut()
+            .insert(noodles::sam::alignment::record::Flags::QC_FAIL);
+        if and_or && !rec.flags().is_unmapped() {
+            crate::io::sam::apply_sam_flag_or_and(rec, params);
+        }
+    }
+}
+
+/// [`apply_rg_per_file`] and [`apply_qc_fail`] over a batch's genomic and
+/// transcriptome records; `read(i)` is the batch's `i`-th read (mate 1 for pairs).
+fn apply_per_read_batch(
+    params: &Parameters,
+    out: &mut BatchOut<AlignmentBatchResults>,
+    read: impl Fn(usize) -> (u32, u8),
+) {
+    for (i, res) in out.iter_mut().enumerate() {
+        if let Ok(o) = res {
+            let filter_flag = read(i).1;
+            apply_qc_fail(&mut o.sam_records.records, filter_flag, params, true);
+            apply_qc_fail(
+                &mut o.transcriptome_records.records,
+                filter_flag,
+                params,
+                false,
+            );
+        }
+    }
+    apply_rg_per_file_batch(params, out, |i| read(i).0);
+}
+
 /// [`apply_rg_per_file`] over a batch's genomic and transcriptome records;
 /// `file_index(i)` is the file of the batch's `i`-th read.
 fn apply_rg_per_file_batch(
@@ -2252,7 +2297,9 @@ fn align_reads_single_end<W: AlignmentWriter + ?Sized>(
                     })
                 })
                 .collect();
-            apply_rg_per_file_batch(params, &mut out, |i| batch[i].file_index);
+            apply_per_read_batch(params, &mut out, |i| {
+                (batch[i].file_index, batch[i].filter_flag)
+            });
             out
         };
         Arc::new(process)
@@ -2876,12 +2923,11 @@ fn align_reads_solo<W: AlignmentWriter + ?Sized>(
                             })
                         })
                         .collect();
-                    if !rg_ids.is_empty() {
-                        for (i, res) in out.iter_mut().enumerate() {
-                            if let Ok(o) = res {
-                                let fi = batch[i].cdna.file_index;
-                                apply_rg_per_file(&mut o.sam_records.records, &rg_ids, fi);
-                            }
+                    for (i, res) in out.iter_mut().enumerate() {
+                        if let Ok(o) = res {
+                            let r = &batch[i].cdna;
+                            apply_qc_fail(&mut o.sam_records.records, r.filter_flag, params, true);
+                            apply_rg_per_file(&mut o.sam_records.records, &rg_ids, r.file_index);
                         }
                     }
                     out
@@ -3401,12 +3447,11 @@ fn align_reads_solo_pe<W: AlignmentWriter + ?Sized>(
                             })
                         })
                         .collect();
-                    if !rg_ids.is_empty() {
-                        for (i, res) in out.iter_mut().enumerate() {
-                            if let Ok(o) = res {
-                                let fi = batch[i].mate1.file_index;
-                                apply_rg_per_file(&mut o.sam_records.records, &rg_ids, fi);
-                            }
+                    for (i, res) in out.iter_mut().enumerate() {
+                        if let Ok(o) = res {
+                            let r = &batch[i].mate1;
+                            apply_qc_fail(&mut o.sam_records.records, r.filter_flag, params, true);
+                            apply_rg_per_file(&mut o.sam_records.records, &rg_ids, r.file_index);
                         }
                     }
                     out
@@ -3994,7 +4039,9 @@ fn align_reads_paired_end<W: AlignmentWriter + ?Sized>(
                     })
                 })
                 .collect();
-            apply_rg_per_file_batch(params, &mut out, |i| batch[i].mate1.file_index);
+            apply_per_read_batch(params, &mut out, |i| {
+                (batch[i].mate1.file_index, batch[i].mate1.filter_flag)
+            });
             out
         };
         Arc::new(process)

src/params/mod.rs

@@ -861,6 +861,10 @@ pub struct Parameters {
     #[arg(long = "outFilterMismatchNoverLmax", default_value_t = 0.3)]
     pub out_filter_mismatch_nover_lmax: f64,
 
+    /// Max ratio of mismatches to read length (both mates)
+    #[arg(long = "outFilterMismatchNoverReadLmax", default_value_t = 1.0)]
+    pub out_filter_mismatch_nover_read_lmax: f64,
+
     /// Min alignment score (absolute)
     #[arg(long = "outFilterScoreMin", default_value_t = 0)]
     pub out_filter_score_min: i32,
@@ -1670,6 +1674,14 @@ impl Parameters {
             .collect()
     }
 
+    /// STAR's per-read `outFilterMismatchNmaxTotal` (`ReadAlign_oneRead.cpp:78`):
+    /// `outFilterMismatchNmax` capped by `outFilterMismatchNoverReadLmax` times the
+    /// read length, both mates after clipping.
+    pub fn mismatch_nmax_total(&self, read_len: usize) -> u32 {
+        let by_len = (self.out_filter_mismatch_nover_read_lmax * read_len as f64) as u32;
+        self.out_filter_mismatch_nmax.min(by_len)
+    }
+
     /// Read group ID emitted on SAM records (the first RG line's `ID:` value).
     /// Returns `None` when no RG line is configured.
     pub fn primary_rg_id(&self) -> Result<Option<String>, crate::error::Error> {

src/quant/transcriptome.rs

@@ -1592,9 +1592,7 @@ fn prepare_for_projection(
 /// STAR's `(nMM + nMM1) > min(outFilterMismatchNmaxTotal,
 /// outFilterMismatchNoverLmax*(Lread-1))` test; `lread` is STAR's `Lread`.
 fn within_mismatch_budget(n_mm: u32, read_len_sum: u32, lread: u32, params: &Parameters) -> bool {
-    // `outFilterMismatchNmaxTotal = min(outFilterMismatchNmax,
-    // outFilterMismatchNoverReadLmax*(readLength[0]+readLength[1]))`, ratio 1.
-    let nmax_total = params.out_filter_mismatch_nmax.min(read_len_sum);
+    let nmax_total = params.mismatch_nmax_total(read_len_sum as usize);
     let rel =
         (params.out_filter_mismatch_nover_lmax * f64::from(lread.saturating_sub(1))).floor() as u32;
     n_mm <= nmax_total.min(rel)