diff --git a/CHANGELOG.md b/CHANGELOG.md index 1a57e84..857f467 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,15 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- FASTA reads are accepted, in one or several lines, mixed with FASTQ, as in STAR. Their records have `*` for quality, and they go to `Unmapped.out.mate*` as FASTA. +- A blank line at the end of a read file no longer stops the run. A blank line with reads after it is an error. See `DIVERGENCE.md` ยง4.8. +- Mates named `SRR1.5.1` / `SRR1.5.2` are accepted as a pair. `SRR1.1` / `SRR1.2` are still refused, since there the number is the read's. +- Both records of a pair carry mate 1's full name as QNAME, as STAR does. A `_1`, `.R1` or `.1` suffix was removed from it. +- Sort spill files stay in the output folder when `--outFileNamePrefix` ends in `/`. They went to the folder above it. +- `--outBAMcompression -1` compresses like level 1, as STAR does. It wrote uncompressed BAM. + +### Bug fixes + - Streamed FASTQ input (FIFOs, `<(...)`, `/dev/stdin`) fails instead of passing for complete. A compressed stream cut in a way the FASTQ parser took for a clean end of file (gzip trailer, member header, or decoded text ending on a diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 6d65f0a..bb7e84c 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -341,6 +341,18 @@ rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand` **Source.** `FastqReader::open_with_command` in `src/io/fastq.rs`; `join_read_files_command` in `src/params/mod.rs`. +### 4.8 A blank line before more reads is an error + +**What STAR does.** Its read loop stops at the first line starting with a space or a newline (observed: `@q1 ... \n\n@q2 ...` maps `q1`, ignores `q2`, exit 0). + +**What rustar-aligner does.** Blank lines at the end of a file are fine, as in STAR. Data after a blank line stops the run with an error naming the file. + +**Why.** Reads after the blank line would be lost with exit code 0. + +**Impact.** Only input with a blank line in the middle. Those runs now fail rather than align part of the file. + +**Source.** `src/io/fastq.rs` (`expect_only_blank_lines`), locked by `blank_lines_end_a_file_but_not_before_more_reads`. + ## 5. Known residual single-read differences These are **not** deliberate divergences. They are the differences `test/parity_suite.sh` finds, kept here so that each one is fixed or explained before v0.2.0 (README, *Comparison with STAR 2.7.11b*). Read names are from that suite's datasets. diff --git a/src/io/bam.rs b/src/io/bam.rs index be8b040..220a648 100644 --- a/src/io/bam.rs +++ b/src/io/bam.rs @@ -36,11 +36,13 @@ impl BufferedBamRecords { /// Convert STAR's `--outBAMcompression` integer (-1..9) to a noodles level. /// -/// STAR mapping: -1 or 0 = uncompressed, 1-9 = deflate levels, default 1. +/// 0 = uncompressed, 1-9 = deflate levels. STAR's -1 compresses like level 1 +/// (measured: 59,925 bytes against 59,924 at level 1 and 635,969 at level 0). fn bgzf_compression(level: i32) -> bgzf::io::writer::CompressionLevel { use bgzf::io::writer::CompressionLevel; match level { - n if n <= 0 => CompressionLevel::NONE, + 0 => CompressionLevel::NONE, + n if n < 0 => CompressionLevel::try_from(1).unwrap_or_default(), n if n >= 9 => CompressionLevel::BEST, n => CompressionLevel::try_from(n as u8).unwrap_or_default(), } @@ -166,6 +168,13 @@ fn sort_temp_dir(params: &Parameters) -> std::path::PathBuf { return std::path::PathBuf::from(¶ms.out_tmp_dir); } let prefix = params.output_path(""); + // `Path::parent` of `out/sub/` is `out`, one level above the output. + if params + .out_file_name_prefix + .ends_with(['/', std::path::MAIN_SEPARATOR]) + { + return prefix; + } match prefix.parent() { Some(parent) if !parent.as_os_str().is_empty() => parent.to_path_buf(), _ => std::path::PathBuf::from("."), @@ -714,6 +723,22 @@ mod tests { use noodles::sam::alignment::record::cigar; use tempfile::NamedTempFile; + #[test] + fn minus_one_compresses_like_level_one_and_zero_does_not() { + assert_eq!(bgzf_compression(-1), bgzf_compression(1)); + assert_ne!(bgzf_compression(0), bgzf_compression(1)); + } + + #[test] + fn spill_files_stay_inside_a_directory_prefix() { + let args = + |prefix: &'static str| ["x", "--readFilesIn", "r.fq", "--outFileNamePrefix", prefix]; + let dir = |prefix| sort_temp_dir(&Parameters::parse_from(args(prefix))); + assert_eq!(dir("/data/run1/"), std::path::Path::new("/data/run1/")); + assert_eq!(dir("/data/run1/s_"), std::path::Path::new("/data/run1")); + assert_eq!(dir("s_"), std::path::Path::new(".")); + } + fn create_test_genome() -> Genome { Genome { transform_blocks: None, diff --git a/src/io/fastq.rs b/src/io/fastq.rs index 6cae978..dfc37fd 100644 --- a/src/io/fastq.rs +++ b/src/io/fastq.rs @@ -25,7 +25,9 @@ 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> { - self.writer.write_all(b"@").map_err(Error::from)?; + // A FASTA input read (no quality) stays FASTA, as STAR writes it. + let marker: &[u8] = if qual.is_empty() { b">" } else { b"@" }; + self.writer.write_all(marker).map_err(Error::from)?; self.writer .write_all(name.as_bytes()) .map_err(Error::from)?; @@ -35,6 +37,9 @@ impl UnmappedFastqWriter { .write_all(&[decode_base(b)]) .map_err(Error::from)?; } + if qual.is_empty() { + return self.writer.write_all(b"\n").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) @@ -63,7 +68,7 @@ pub struct EncodedRead { /// A paired-end read from two FASTQ files #[derive(Debug, Clone)] pub struct PairedRead { - /// Base read name (without /1 or /2 suffix) + /// Mate 1's read name, which both mates' records carry pub name: String, /// First mate in pair pub mate1: EncodedRead, @@ -253,6 +258,15 @@ impl FastqReader { } fn next_in_current_file(&mut self) -> Result, Error> { + match self.peek_first_byte()? { + None => return Ok(None), + Some(b'>') => return self.next_fasta().map(Some), + Some(b'\n' | b'\r' | b' ') => { + self.expect_only_blank_lines()?; + return Ok(None); + } + Some(_) => {} + } match self.inner.records().next() { Some(Ok(record)) => { // STAR's check (`readLoad.cpp`). It is also what catches plain @@ -271,27 +285,10 @@ impl FastqReader { )); } - let name = std::str::from_utf8(record.name()) - .map_err(|e| { - Error::from(std::io::Error::new( - std::io::ErrorKind::InvalidData, - format!("invalid UTF-8 in read name: {e}"), - )) - })? - .to_string(); - + let name = self.read_name(record.name())?; + self.check_not_empty(record.sequence(), &name)?; let sequence = record.sequence().iter().map(|&b| encode_base(b)).collect(); - let name = match self - .name_separators - .iter() - .filter_map(|&sep| name.as_bytes().iter().position(|&b| b == sep)) - .min() - { - Some(cut) => name[..cut].to_string(), - None => name, - }; - let quality = if self.qual_shift == 0 { record.quality_scores().to_vec() } else { @@ -314,6 +311,108 @@ impl FastqReader { } } + /// The next byte of the current file, not consumed; `None` at its end. + fn peek_first_byte(&mut self) -> Result, Error> { + let buf = self + .inner + .get_mut() + .fill_buf() + .map_err(|e| Error::io(e, &self.path))?; + Ok(buf.first().copied()) + } + + /// A read name as the rest of the run sees it, cut at `--readNameSeparator`. + fn read_name(&self, raw: &[u8]) -> Result { + let name = std::str::from_utf8(raw).map_err(|e| { + Error::from(std::io::Error::new( + std::io::ErrorKind::InvalidData, + format!("invalid UTF-8 in read name: {e}"), + )) + })?; + let cut = self + .name_separators + .iter() + .filter_map(|&sep| name.as_bytes().iter().position(|&b| b == sep)) + .min() + .unwrap_or(name.len()); + Ok(name[..cut].to_string()) + } + + /// STAR stops on a read with no bases (`readLoad.cpp`, exit 104). + fn check_not_empty(&self, sequence: &[u8], name: &str) -> Result<(), Error> { + if !sequence.is_empty() { + return Ok(()); + } + Err(Error::io( + std::io::Error::new( + std::io::ErrorKind::InvalidData, + format!("read {name}: short read sequence line: 0"), + ), + &self.path, + )) + } + + /// One FASTA read: a `>name` line, then its sequence over any number of lines. + /// Quality stays empty, which the SAM and BAM writers output as `*`, as STAR does. + fn next_fasta(&mut self) -> Result { + let mut line = Vec::new(); + self.inner + .get_mut() + .read_until(b'\n', &mut line) + .map_err(|e| Error::io(e, &self.path))?; + let header = line[1..].trim_ascii_end(); + let raw_name = header.split(u8::is_ascii_whitespace).next().unwrap_or(&[]); + let name = self.read_name(raw_name)?; + + let mut sequence = Vec::new(); + while let Some(first) = self.peek_first_byte()? { + if matches!(first, b'>' | b'@' | b'\n' | b'\r' | b' ') { + break; + } + line.clear(); + self.inner + .get_mut() + .read_until(b'\n', &mut line) + .map_err(|e| Error::io(e, &self.path))?; + sequence.extend(line.trim_ascii_end().iter().map(|&b| encode_base(b))); + } + self.check_not_empty(&sequence, &name)?; + Ok(EncodedRead { + name, + sequence, + quality: Vec::new(), + file_index: self.file_index, + }) + } + + /// Blank lines are fine at the end of a file. STAR stops at the first one and + /// exits 0 even with reads after it; those would be lost, so that is an error + /// here (`DIVERGENCE.md` 4.8). + fn expect_only_blank_lines(&mut self) -> Result<(), Error> { + let mut line = Vec::new(); + loop { + line.clear(); + let n = self + .inner + .get_mut() + .read_until(b'\n', &mut line) + .map_err(|e| Error::io(e, &self.path))?; + if n == 0 { + return Ok(()); + } + if !line.trim_ascii().is_empty() { + return Err(Error::io( + std::io::Error::new( + std::io::ErrorKind::InvalidData, + "data after a blank line; remove the blank line, or reads after \ + it would be lost", + ), + &self.path, + )); + } + } + } + /// Read a batch of encoded reads for parallel processing /// /// # Arguments @@ -404,7 +503,8 @@ impl PairedFastqReader { } Ok(Some(PairedRead { - name: name1_base, + // STAR's QNAME is mate 1's name; the suffix only matters above. + name: read1.name.clone(), mate1: read1, mate2: read2, })) @@ -446,6 +546,7 @@ impl PairedFastqReader { /// - /1 or /2 (Illumina convention) /// - .R1 or .R2 (alternative convention) /// - _1 or _2 (another convention) +/// - .1 or .2 after a read number (`SRR1.5.1`) /// - space and everything after (e.g., "READ_NAME 1:N:0:0" -> "READ_NAME") /// /// # Arguments @@ -469,6 +570,11 @@ pub fn strip_mate_suffix(name: &str) -> String { name[..name.len() - 3].to_string() } else if name.ends_with("_1") || name.ends_with("_2") { name[..name.len() - 2].to_string() + } else if (name.ends_with(".1") || name.ends_with(".2")) && name[..name.len() - 2].contains('.') + { + // `SRR1.5.1`. Not `SRR1.1`: there `.1` is the spot number, and stripping it + // would hide mates that are out of step. + name[..name.len() - 2].to_string() } else { name.to_string() } @@ -541,7 +647,9 @@ pub fn clip_read(seq: &[u8], qual: &[u8], clip5p: usize, clip3p: usize) -> (Vec< let start = clip5p; let end = len - clip3p; - (seq[start..end].to_vec(), qual[start..end].to_vec()) + // FASTA reads have no quality. + let clipped_qual = qual.get(start..end).unwrap_or_default().to_vec(); + (seq[start..end].to_vec(), clipped_qual) } #[cfg(test)] @@ -797,6 +905,77 @@ mod tests { assert!(reader.next_encoded().unwrap().is_none()); } + fn reader_over(contents: &str) -> FastqReader { + let mut f = NamedTempFile::new().unwrap(); + f.write_all(contents.as_bytes()).unwrap(); + FastqReader::open(f.path(), None).unwrap() + } + + #[test] + fn fasta_reads_have_empty_quality_and_may_span_lines() { + let mut r = reader_over(">f1 desc\nAC\nGT\n>f2/1\nTT\n@q3\nAA\n+\nII\n>f4\nCC"); + let f1 = r.next_encoded().unwrap().unwrap(); + assert_eq!( + (f1.name.as_str(), f1.sequence, f1.quality), + ("f1", vec![0, 1, 2, 3], vec![]) + ); + assert_eq!(r.next_encoded().unwrap().unwrap().name, "f2"); + assert_eq!(r.next_encoded().unwrap().unwrap().quality, b"II"); + assert_eq!(r.next_encoded().unwrap().unwrap().sequence, vec![1, 1]); + assert!(r.next_encoded().unwrap().is_none()); + } + + /// STAR exits 104 on a read with no bases, FASTA or FASTQ. + #[test] + fn a_read_with_no_bases_is_an_error() { + let mut r = reader_over(">e0\nACGT\n>e1\n>e2\nAC\n"); + assert_eq!(r.next_encoded().unwrap().unwrap().name, "e0"); + let err = r.next_encoded().unwrap_err().to_string(); + assert!(err.contains("e1: short read sequence line: 0"), "{err}"); + let err = reader_over("@q\n\n+\n\n") + .next_encoded() + .unwrap_err() + .to_string(); + assert!(err.contains("short read sequence line: 0"), "{err}"); + } + + #[test] + fn blank_lines_end_a_file_but_not_before_more_reads() { + let mut r = reader_over("@a\nAC\n+\nII\n\n \n\r\n"); + assert!(r.next_encoded().unwrap().is_some()); + assert!(r.next_encoded().unwrap().is_none()); + + let mut r = reader_over("@a\nAC\n+\nII\n\n@b\nAC\n+\nII\n"); + assert!(r.next_encoded().unwrap().is_some()); + let err = r.next_encoded().unwrap_err().to_string(); + assert!(err.contains("blank line"), "{err}"); + } + + #[test] + fn clip_read_of_a_fasta_read_has_no_quality() { + let (seq, qual) = clip_read(&[0, 1, 2, 3], &[], 1, 1); + assert_eq!((seq, qual), (vec![1, 2], vec![])); + } + + #[test] + fn dotted_mate_suffix_needs_a_read_number_before_it() { + assert_eq!(strip_mate_suffix("SRR1.5.1"), "SRR1.5"); + assert_eq!(strip_mate_suffix("SRR1.5.2"), "SRR1.5"); + // Here `.1` and `.2` are spot numbers: mates out of step must still differ. + assert_eq!(strip_mate_suffix("SRR1.1"), "SRR1.1"); + assert_ne!(strip_mate_suffix("SRR1.1"), strip_mate_suffix("SRR1.2")); + } + + #[test] + fn pair_is_named_after_mate1_as_star_does() { + let mut f1 = NamedTempFile::new().unwrap(); + let mut f2 = NamedTempFile::new().unwrap(); + f1.write_all(b"@p_1\nAC\n+\nII\n").unwrap(); + f2.write_all(b"@p_2\nAC\n+\nII\n").unwrap(); + let mut r = PairedFastqReader::open(f1.path(), f2.path(), None).unwrap(); + assert_eq!(r.next_paired().unwrap().unwrap().name, "p_1"); + } + #[test] fn test_strip_mate_suffix_slash() { assert_eq!(strip_mate_suffix("read123/1"), "read123"); diff --git a/src/io/sam.rs b/src/io/sam.rs index 40a03e5..5a9c9c7 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -132,6 +132,15 @@ pub(crate) fn insert_mate_cigar_tags(rec1: &mut RecordBuf, rec2: &mut RecordBuf) .insert(tag, Value::String(BString::from(c1))); } +/// The aligned part of a read's quality. A FASTA read has none, whatever its clips. +fn aligned_quality(qual: &[u8], clip5p: usize, clip3p: usize) -> &[u8] { + if qual.is_empty() { + qual + } else { + &qual[clip5p..qual.len() - clip3p] + } +} + /// SAM file writer pub struct SamWriter { writer: sam::io::Writer>, @@ -330,7 +339,7 @@ impl SamWriter { // (CIGAR core, MD, NM) is built from that aligned slice, then apply_read_clips // restores the full read + soft-clips. clip5p/clip3p == 0 => aligned == full. let aligned_seq = &read_seq[clip5p..read_seq.len() - clip3p]; - let aligned_qual = &read_qual[clip5p..read_qual.len() - clip3p]; + let aligned_qual = aligned_quality(read_qual, clip5p, clip3p); let n_alignments = transcripts.len(); let max_output = if params.out_sam_mult_nmax < 0 { @@ -588,7 +597,7 @@ impl SamWriter { }; // Aligned slice of the mapped mate (core SEQ/CIGAR/MD built from this). let aligned_mapped_seq = &mapped_seq[mapped_clip5p..mapped_seq.len() - mapped_clip3p]; - let aligned_mapped_qual = &mapped_qual[mapped_clip5p..mapped_qual.len() - mapped_clip3p]; + let aligned_mapped_qual = aligned_quality(mapped_qual, mapped_clip5p, mapped_clip3p); // --- Build mapped mate record --- let mut mapped_rec = RecordBuf::default(); @@ -1977,7 +1986,7 @@ fn build_paired_mate_record( // Core SEQ/QUAL/MD are built from the aligned (clipped) slice; apply_read_clips // below restores the full mate read + soft-clips (clip==0 => aligned == mate_seq). let aligned_seq = &mate_seq[clip5p..mate_seq.len() - clip3p]; - let aligned_qual = &mate_qual[clip5p..mate_qual.len() - clip3p]; + let aligned_qual = aligned_quality(mate_qual, clip5p, clip3p); // Sequence and quality scores (reverse complement for reverse strand) if transcript.is_reverse {