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())); +}