FASTA reads, blank lines and mate names as STAR; spill folder; --outBAMcompression -1

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

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.

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(&params.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,

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<Option<EncodedRead>, 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<Option<u8>, 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<String, Error> {
+        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<EncodedRead, Error> {
+        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");