SEQ and Unmapped.out.mate keep the read's letters (lowercase, IUPAC)

CHANGELOG.md

@@ -232,6 +232,15 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Bug fixes
 
+- SEQ keeps the read's own letters, as STAR does: lowercase bases and IUPAC
+  codes other than N (`R`, `y`, ...) were written as uppercase `ACGTN`. On
+  the reverse strand they are reverse complemented with their case, as
+  STAR's `revComplementNucleotides`. BAM keeps the IUPAC codes (its 4-bit
+  code has no case). `Unmapped.out.mate*` has the read's letters too, and
+  for pairs now the whole mate, as STAR writes it: under `--clip3pNbases` or
+  `--clip5pNbases` it had the clipped mate. Reads in plain `ACGTN`, the usual
+  case, take no extra work.
+
 - **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

src/io/fastq.rs

@@ -22,19 +22,15 @@ impl UnmappedFastqWriter {
         })
     }
 
-    /// Write one FASTQ record. `seq` is in genome encoding (0=A,1=C,2=G,3=T,4=N).
-    /// `qual` is raw FASTQ quality bytes.
-    pub fn write_record(&mut self, name: &str, seq: &[u8], qual: &[u8]) -> Result<(), Error> {
+    /// Write one FASTQ record. `letters` are the read's bases as text
+    /// ([`EncodedRead::letters`]), `qual` raw FASTQ quality bytes.
+    pub fn write_record(&mut self, name: &str, letters: &[u8], qual: &[u8]) -> Result<(), Error> {
         self.writer.write_all(b"@").map_err(Error::from)?;
         self.writer
             .write_all(name.as_bytes())
             .map_err(Error::from)?;
         self.writer.write_all(b"\n").map_err(Error::from)?;
-        for &b in seq {
-            self.writer
-                .write_all(&[decode_base(b)])
-                .map_err(Error::from)?;
-        }
+        self.writer.write_all(letters).map_err(Error::from)?;
         self.writer.write_all(b"\n+\n").map_err(Error::from)?;
         self.writer.write_all(qual).map_err(Error::from)?;
         self.writer.write_all(b"\n").map_err(Error::from)
@@ -58,6 +54,20 @@ pub struct EncodedRead {
     /// Index of the `--readFilesIn` file this read came from (STAR
     /// `readFilesIndex`); selects the read group.
     pub file_index: u32,
+    /// The bases as written in the FASTQ, kept only when they are not all
+    /// `ACGTN` (lowercase, other IUPAC codes): STAR writes them as they are.
+    pub raw_letters: Option<Vec<u8>>,
+}
+
+impl EncodedRead {
+    /// The bases as text, as STAR writes them in SEQ and `Unmapped.out.mate*`
+    /// (its `Read0`).
+    pub fn letters(&self) -> Vec<u8> {
+        match &self.raw_letters {
+            Some(raw) => raw.clone(),
+            None => self.sequence.iter().map(|&b| decode_base(b)).collect(),
+        }
+    }
 }
 
 /// A paired-end read from two FASTQ files
@@ -209,6 +219,11 @@ impl FastqReader {
                     .to_string();
 
                 let sequence = record.sequence().iter().map(|&b| encode_base(b)).collect();
+                let raw_letters = record
+                    .sequence()
+                    .iter()
+                    .any(|b| !b"ACGTN".contains(b))
+                    .then(|| record.sequence().to_vec());
 
                 let name = match self
                     .name_separators
@@ -235,6 +250,7 @@ impl FastqReader {
                     sequence,
                     quality,
                     file_index: self.file_index,
+                    raw_letters,
                 }))
             }
             Some(Err(e)) => Err(Error::from(e)),
@@ -436,6 +452,38 @@ pub fn decode_base(encoded: u8) -> u8 {
     }
 }
 
+/// Reverse complement of read text as STAR's `revComplementNucleotides`
+/// (`SequenceFuns.cpp:16-54`): case and IUPAC codes kept, other bytes unchanged.
+pub fn revcomp_letters(letters: &[u8]) -> Vec<u8> {
+    letters
+        .iter()
+        .rev()
+        .map(|&b| {
+            let c = match b.to_ascii_uppercase() {
+                b'A' => b'T',
+                b'C' => b'G',
+                b'G' => b'C',
+                b'T' => b'A',
+                b'R' => b'Y',
+                b'Y' => b'R',
+                b'K' => b'M',
+                b'M' => b'K',
+                b'B' => b'V',
+                b'D' => b'H',
+                b'V' => b'B',
+                b'H' => b'D',
+                b'N' | b'S' | b'W' => b.to_ascii_uppercase(),
+                _ => return b,
+            };
+            if b.is_ascii_lowercase() {
+                c.to_ascii_lowercase()
+            } else {
+                c
+            }
+        })
+        .collect()
+}
+
 /// Complement an encoded base (A=0↔T=3, C=1↔G=2, N=4→N=4).
 pub fn complement_base(encoded: u8) -> u8 {
     match encoded {
@@ -496,6 +544,15 @@ mod tests {
         assert_eq!(encode_base(b'S'), 4);
     }
 
+    /// STAR's `revComplementNucleotides` table (`SequenceFuns.cpp:16-54`).
+    #[test]
+    fn revcomp_letters_as_star() {
+        assert_eq!(
+            revcomp_letters(b"ACGTNRYKMSWBDVHacgtnrykmswbdvh.X"),
+            b"X.dbhvwskmrynacgtDBHVWSKMRYNACGT".to_vec()
+        );
+    }
+
     #[test]
     fn test_decode_base() {
         assert_eq!(decode_base(0), b'A');

src/io/sam.rs

@@ -1399,6 +1399,47 @@ pub fn apply_solo_read_info(
 /// Apply `--outSAMflagOR` / `--outSAMflagAND` to a mapped record's FLAG:
 /// `(FLAG & flagAND) | flagOR`. Matches STAR/STAR-rs, which apply this only to
 /// mapped-mate records; unmapped and transcriptome-BAM records are untouched.
+/// Put a read's own letters (lowercase, IUPAC codes) back in SEQ, which the
+/// builders write from the ACGTN encoding. STAR writes its `Read0` text, reverse
+/// complemented by `revComplementNucleotides` on the reverse strand
+/// (`ReadAlign_outputTranscriptSAM.cpp:194-204`, `ReadAlign_alignBAM.cpp:489-506`).
+/// `mate1`/`mate2` are [`crate::io::fastq::EncodedRead::raw_letters`]; a record
+/// whose SEQ is not that read's bases is left as it is.
+pub fn restore_read_letters(records: &mut [RecordBuf], mate1: Option<&[u8]>, mate2: Option<&[u8]>) {
+    use crate::io::fastq::{encode_base, revcomp_letters};
+    for rec in records {
+        let flags = rec.flags();
+        let Some(letters) = (if flags.is_last_segment() {
+            mate2
+        } else {
+            mate1
+        }) else {
+            continue;
+        };
+        let full = if flags.is_reverse_complemented() && !flags.is_unmapped() {
+            revcomp_letters(letters)
+        } else {
+            letters.to_vec()
+        };
+        let hard5 = rec
+            .cigar()
+            .as_ref()
+            .first()
+            .filter(|op| op.kind() == cigar::op::Kind::HardClip)
+            .map_or(0, |op| op.len());
+        let Some(seq) = full.get(hard5..hard5 + rec.sequence().len()) else {
+            continue;
+        };
+        let same_bases = seq
+            .iter()
+            .zip(rec.sequence().as_ref())
+            .all(|(&a, &b)| decode_base(encode_base(a)) == b);
+        if same_bases {
+            *rec.sequence_mut() = Sequence::from(seq.to_vec());
+        }
+    }
+}
+
 pub(crate) fn apply_sam_flag_or_and(record: &mut RecordBuf, params: &Parameters) {
     let or_bits = params.out_sam_flag_or as u16;
     let and_bits = params.out_sam_flag_and.min(u16::MAX as u32) as u16;
@@ -2072,6 +2113,33 @@ mod tests {
     use noodles::sam::alignment::record::cigar;
     use tempfile::NamedTempFile;
 
+    #[test]
+    fn restore_read_letters_after_hard_clip_and_not_on_other_bases() {
+        use cigar::op::{Kind, Op};
+        use noodles::sam::alignment::record::Flags;
+        let rec = |flags: Flags, ops: Vec<Op>, seq: &[u8]| {
+            let mut r = RecordBuf::default();
+            *r.flags_mut() = flags;
+            *r.cigar_mut() = ops.into_iter().collect();
+            *r.sequence_mut() = Sequence::from(seq.to_vec());
+            r
+        };
+        let mut records = vec![
+            // Reverse strand, 3 bases hard-clipped: SEQ is that slice of the
+            // reverse complement `ttRacgt`.
+            rec(
+                Flags::REVERSE_COMPLEMENTED,
+                vec![Op::new(Kind::HardClip, 3), Op::new(Kind::Match, 4)],
+                b"ACGT",
+            ),
+            // Not this read's bases: left alone.
+            rec(Flags::empty(), vec![Op::new(Kind::Match, 7)], b"AAAAAAA"),
+        ];
+        restore_read_letters(&mut records, Some(b"acgtYaa"), None);
+        assert_eq!(records[0].sequence().as_ref(), b"acgt");
+        assert_eq!(records[1].sequence().as_ref(), b"AAAAAAA");
+    }
+
     fn make_test_genome() -> Genome {
         Genome {
             transform_blocks: None,

src/lib.rs

@@ -833,6 +833,23 @@ fn apply_rg_per_file(
     }
 }
 
+/// [`crate::io::sam::restore_read_letters`] over a batch's genomic and
+/// transcriptome records; `letters(i)` are the `i`-th read's mates' raw letters.
+fn restore_letters_batch<'a>(
+    out: &mut BatchOut<AlignmentBatchResults>,
+    letters: impl Fn(usize) -> (Option<&'a [u8]>, Option<&'a [u8]>),
+) {
+    for (i, res) in out.iter_mut().enumerate() {
+        let (m1, m2) = letters(i);
+        if let Ok(o) = res
+            && (m1.is_some() || m2.is_some())
+        {
+            crate::io::sam::restore_read_letters(&mut o.sam_records.records, m1, m2);
+            crate::io::sam::restore_read_letters(&mut o.transcriptome_records.records, m1, m2);
+        }
+    }
+}
+
 /// [`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(
@@ -2015,11 +2032,7 @@ fn align_reads_single_end<W: AlignmentWriter + ?Sized>(
                             buffer.push(record);
                         }
                         let unmapped_m1 = if write_unmapped_fastq {
-                            vec![(
-                                out_read_name.clone(),
-                                read.sequence.clone(),
-                                read.quality.clone(),
-                            )]
+                            vec![(out_read_name.clone(), read.letters(), read.quality.clone())]
                         } else {
                             Vec::new()
                         };
@@ -2231,11 +2244,7 @@ fn align_reads_single_end<W: AlignmentWriter + ?Sized>(
                         TrRecords::default()
                     };
                     let unmapped_m1 = if write_unmapped_fastq && is_unmapped_se {
-                        vec![(
-                            out_read_name.clone(),
-                            read.sequence.clone(),
-                            read.quality.clone(),
-                        )]
+                        vec![(out_read_name.clone(), read.letters(), read.quality.clone())]
                     } else {
                         Vec::new()
                     };
@@ -2252,6 +2261,7 @@ fn align_reads_single_end<W: AlignmentWriter + ?Sized>(
                     })
                 })
                 .collect();
+            restore_letters_batch(&mut out, |i| (batch[i].raw_letters.as_deref(), None));
             apply_rg_per_file_batch(params, &mut out, |i| batch[i].file_index);
             out
         };
@@ -2876,12 +2886,15 @@ 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 {
+                            crate::io::sam::restore_read_letters(
+                                &mut o.sam_records.records,
+                                batch[i].cdna.raw_letters.as_deref(),
+                                None,
+                            );
+                            let fi = batch[i].cdna.file_index;
+                            apply_rg_per_file(&mut o.sam_records.records, &rg_ids, fi);
                         }
                     }
                     out
@@ -3401,12 +3414,15 @@ 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 {
+                            crate::io::sam::restore_read_letters(
+                                &mut o.sam_records.records,
+                                batch[i].mate1.raw_letters.as_deref(),
+                                batch[i].mate2.raw_letters.as_deref(),
+                            );
+                            let fi = batch[i].mate1.file_index;
+                            apply_rg_per_file(&mut o.sam_records.records, &rg_ids, fi);
                         }
                     }
                     out
@@ -3646,12 +3662,12 @@ fn align_reads_paired_end<W: AlignmentWriter + ?Sized>(
                             (
                                 vec![(
                                     fastx_name1.clone(),
-                                    paired_read.mate1.sequence.clone(),
+                                    paired_read.mate1.letters(),
                                     paired_read.mate1.quality.clone(),
                                 )],
                                 vec![(
                                     fastx_name2.clone(),
-                                    paired_read.mate2.sequence.clone(),
+                                    paired_read.mate2.letters(),
                                     paired_read.mate2.quality.clone(),
                                 )],
                             )
@@ -3972,8 +3988,16 @@ fn align_reads_paired_end<W: AlignmentWriter + ?Sized>(
                         let pair_unmapped = results.is_empty() || has_half_mapped;
                         if pair_unmapped {
                             (
-                                vec![(fastx_name1.clone(), m1_seq.clone(), m1_qual.clone())],
-                                vec![(fastx_name2.clone(), m2_seq.clone(), m2_qual.clone())],
+                                vec![(
+                                    fastx_name1.clone(),
+                                    paired_read.mate1.letters(),
+                                    paired_read.mate1.quality.clone(),
+                                )],
+                                vec![(
+                                    fastx_name2.clone(),
+                                    paired_read.mate2.letters(),
+                                    paired_read.mate2.quality.clone(),
+                                )],
                             )
                         } else {
                             (Vec::new(), Vec::new())
@@ -3994,6 +4018,12 @@ fn align_reads_paired_end<W: AlignmentWriter + ?Sized>(
                     })
                 })
                 .collect();
+            restore_letters_batch(&mut out, |i| {
+                (
+                    batch[i].mate1.raw_letters.as_deref(),
+                    batch[i].mate2.raw_letters.as_deref(),
+                )
+            });
             apply_rg_per_file_batch(params, &mut out, |i| batch[i].mate1.file_index);
             out
         };