From 0862c003665aa00e96394505d9e558f03d3ab0ec Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 10:59:02 +0100 Subject: [PATCH] SEQ and Unmapped.out.mate keep the read's letters Lowercase bases and IUPAC codes are written as in the FASTQ, reverse complemented with case on the reverse strand (STAR's Read0 and revComplementNucleotides). Paired Unmapped.out.mate records have the whole mate, not the clipped one. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 9 +++ src/io/fastq.rs | 73 +++++++++++++++++++--- src/io/sam.rs | 68 ++++++++++++++++++++ src/lib.rs | 82 ++++++++++++++++-------- src/solo/mod.rs | 8 +++ tests/alignment_features.rs | 121 ++++++++++++++++++++++++++++++++++++ 6 files changed, 327 insertions(+), 34 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ab8f0e..06581e1 100644 --- a/CHANGELOG.md +++ b/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`, `/Features.stats` and `/UMIperCellSorted.txt` were not written, SJ and Velocyto had no `Summary.csv`, and the Gene/GeneFull diff --git a/src/io/fastq.rs b/src/io/fastq.rs index 48f9673..342eeb0 100644 --- a/src/io/fastq.rs +++ b/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>, +} + +impl EncodedRead { + /// The bases as text, as STAR writes them in SEQ and `Unmapped.out.mate*` + /// (its `Read0`). + pub fn letters(&self) -> Vec { + 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 { + 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'); diff --git a/src/io/sam.rs b/src/io/sam.rs index 40a03e5..a3e8110 100644 --- a/src/io/sam.rs +++ b/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, 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, diff --git a/src/lib.rs b/src/lib.rs index 1d0c993..05a44f2 100644 --- a/src/lib.rs +++ b/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, + 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( 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( 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( }) }) .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( }) }) .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( }) }) .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( ( 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( 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( }) }) .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 }; diff --git a/src/solo/mod.rs b/src/solo/mod.rs index b91d1fb..cc36bbd 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -2249,6 +2249,7 @@ mod tests { sequence: seq.bytes().map(encode_base).collect(), quality: qual.bytes().collect(), file_index: 0, + raw_letters: None, } } @@ -2750,6 +2751,7 @@ mod tests { sequence: seq.bytes().map(encode_base).collect(), quality: vec![b'I'; seq.len()], file_index: 0, + raw_letters: None, }; // Default: exactly CB+UMI = 26 bases. let checked = BarcodeReadLength { @@ -2890,6 +2892,7 @@ mod tests { sequence: vec![2, 2, 2, 2, 0, 1, 2, 3, 3, 3, 3, 3], quality: vec![b'I'; 12], file_index: 0, + raw_letters: None, }; let bc = layout.extract(&read).expect("adapter present"); assert_eq!(bc.cb_seq, vec![2, 2, 2, 2]); @@ -2901,6 +2904,7 @@ mod tests { sequence: vec![1, 1, 2, 2, 2, 2, 0, 1, 2, 3, 3, 3, 3, 3], quality: vec![b'I'; 14], file_index: 0, + raw_letters: None, }; let bc = layout.extract(&read).expect("adapter present, shifted"); assert_eq!(bc.cb_seq, vec![2, 2, 2, 2]); @@ -2912,6 +2916,7 @@ mod tests { sequence: vec![2, 2, 2, 2, 0, 1, 2, 0, 3, 3, 3, 3], quality: vec![b'I'; 12], file_index: 0, + raw_letters: None, }; assert!(layout.extract(&read).is_ok()); @@ -2922,6 +2927,7 @@ mod tests { sequence: vec![2; 12], quality: vec![b'I'; 12], file_index: 0, + raw_letters: None, }; assert_eq!(layout.extract(&read), Err(BarcodeMiss::NoAdapter)); } @@ -2999,6 +3005,7 @@ mod tests { sequence: vec![2; 12], quality: vec![b'I'; 12], file_index: 0, + raw_letters: None, }; assert_eq!(layout.extract(&read), Err(BarcodeMiss::NoAdapter)); } @@ -3029,6 +3036,7 @@ mod tests { sequence: vec![0, 1, 2, 3, 3, 3, 3], quality: vec![b'I'; 7], file_index: 0, + raw_letters: None, }; assert_eq!(layout.extract(&read), Err(BarcodeMiss::NoCb)); } diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 7e71c04..7c9e0d5 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -3673,3 +3673,124 @@ fn solo_out_files(output_dir: &Path) -> Vec { out.sort(); out } + +/// Reverse complement keeping case and the IUPAC codes used below. +fn rc_letters(seq: &[u8]) -> Vec { + seq.iter() + .rev() + .map(|&b| { + let c = match b.to_ascii_uppercase() { + b'A' => b'T', + b'T' => b'A', + b'C' => b'G', + b'G' => b'C', + b'R' => b'Y', + b'Y' => b'R', + b'N' => b'N', + _ => panic!("no complement for {b}"), + }; + if b.is_ascii_lowercase() { + c.to_ascii_lowercase() + } else { + c + } + }) + .collect() +} + +/// SEQ and `Unmapped.out.mate*` keep the read's letters, lowercase and IUPAC +/// codes, as STAR 2.7.11b writes them (its `Read0`, reverse complemented with +/// case kept on the reverse strand; `ReadAlign_outputTranscriptSAM.cpp:194-204`, +/// `ReadAlign_outputAlignments.cpp:267`). `Unmapped.out.mate*` has the whole mate, +/// before clipping. Checked against STAR's output on the same input. +#[test] +fn test_seq_keeps_read_letters_like_star() { + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", None); + + // Forward read: lowercase head, an R. Reverse read: the reverse complement of + // a forward text with a lowercase tail and a y, so its SEQ is that text. + let mut fwd = genome[1000..1100].to_vec(); + fwd[..20].make_ascii_lowercase(); + fwd[50] = b'R'; + let mut rev_seq = genome[3000..3100].to_vec(); + rev_seq[70..].make_ascii_lowercase(); + rev_seq[39] = b'y'; + let noise = lcg_seq(777, 100).to_ascii_lowercase(); + // Pair: lowercase mate 1, mate 2 with an n. + let m1 = genome[5000..5100].to_ascii_lowercase(); + let mut m2_seq = genome[5200..5300].to_vec(); + m2_seq[10] = b'n'; + let m2 = rc_letters(&m2_seq); + let noise2 = lcg_seq(999, 100).to_ascii_lowercase(); + + let write = |name: &str, recs: &[(&str, &[u8])]| { + let path = tmpdir.path().join(name); + let mut text = String::new(); + for (n, s) in recs { + let s = String::from_utf8_lossy(s); + text.push_str(&format!("@{n}\n{s}\n+\n{}\n", "I".repeat(s.len()))); + } + fs::write(&path, text).unwrap(); + path.to_str().unwrap().to_string() + }; + let se = write( + "se.fq", + &[ + ("fwd", &fwd), + ("rev", &rc_letters(&rev_seq)), + ("noise", &noise), + ], + ); + let pe1 = write("pe1.fq", &[("pair", &m1), ("noise", &noise)]); + let pe2 = write("pe2.fq", &[("pair", &m2), ("noise", &noise2)]); + + let run = |out: &str, reads: &[&str]| -> (Vec<(String, String)>, String) { + let prefix = format!("{}/{out}/", tmpdir.path().display()); + cargo_bin_cmd!("rustar-aligner") + .args(["--genomeDir", genome_dir.to_str().unwrap(), "--readFilesIn"]) + .args(reads) + .args(["--outFileNamePrefix", &prefix]) + .args(["--outSAMunmapped", "Within", "--outReadsUnmapped", "Fastx"]) + .args(["--clip3pNbases"]) + .args(reads.iter().map(|_| "10")) + .assert() + .success(); + let sam = fs::read_to_string(format!("{prefix}Aligned.out.sam")).unwrap(); + let mut recs: Vec<(String, String)> = sam + .lines() + .filter(|l| !l.starts_with('@')) + .map(|l| { + let f: Vec<&str> = l.split('\t').collect(); + (format!("{} {}", f[0], f[1]), f[9].to_string()) + }) + .collect(); + recs.sort(); + let mate1 = fs::read_to_string(format!("{prefix}Unmapped.out.mate1")).unwrap(); + (recs, mate1) + }; + let s = |b: &[u8]| String::from_utf8(b.to_vec()).unwrap(); + + let (recs, mate1) = run("se", &[&se]); + let want = vec![ + ("fwd 0".to_string(), s(&fwd)), + ("noise 4".to_string(), s(&noise)), + ("rev 16".to_string(), s(&rev_seq)), + ]; + assert_eq!(recs, want); + assert_eq!(mate1.lines().nth(1), Some(s(&noise).as_str())); + + let (recs, mate1) = run("pe", &[&pe1, &pe2]); + let want = vec![ + ("noise 141".to_string(), s(&noise2)), + ("noise 77".to_string(), s(&noise)), + ("pair 147".to_string(), s(&m2_seq)), + ("pair 99".to_string(), s(&m1)), + ]; + assert_eq!(recs, want); + // The whole mate, not the clipped 90 bases. + assert_eq!(mate1.lines().nth(1), Some(s(&noise).as_str())); +}