fastx/

reader.rs

//! Streaming FASTA/FASTQ reader.

use std::fmt;
use std::fs::File;
use std::io::{self, BufReader, Read};
use std::path::Path;

use crate::borrowed::SequenceRef;
use crate::error::{Error, ParseError, Result};
use crate::format::{Compression, Format};
use crate::qual::{self, QualityEncoding};
use crate::record::Sequence;

/// Default read buffer, large enough to hold a full 150 bp read set line and to
/// keep syscall overhead negligible.
pub const DEFAULT_BUFFER_SIZE: usize = 128 * 1024;

/// Smallest buffer we will honour; smaller values just make parsing slower.
const MIN_BUFFER_SIZE: usize = 4 * 1024;

/// A streaming reader for FASTA and FASTQ.
///
/// The reader owns a single growable buffer and never holds more than one record
/// plus the buffer in memory, so a 300 GB FASTQ costs the same as a 300 byte one.
/// Records with lines longer than the buffer are handled by growing the buffer on
/// demand, which means a chromosome-on-one-line FASTA works too.
///
/// # Examples
///
/// ```
/// use fastx::FastxReader;
///
/// let data = b">read1 first\nACGT\nACGT\n>read2\nTTTT\n";
/// let mut reader = FastxReader::new(&data[..]);
///
/// let first = reader.next().unwrap()?;
/// assert_eq!(first.id, "read1");
/// assert_eq!(first.description.as_deref(), Some("first"));
/// assert_eq!(first.seq, b"ACGTACGT"); // multi-line sequences are joined
///
/// assert_eq!(reader.count(), 1); // one record left
/// # Ok::<(), fastx::Error>(())
/// ```
///
/// Reuse one record to parse without allocating:
///
/// ```
/// use fastx::{FastxReader, Sequence};
///
/// let data = b"@r1\nACGT\n+\nIIII\n@r2\nTTTT\n+\n!!!!\n";
/// let mut reader = FastxReader::new(&data[..]);
/// let mut record = Sequence::default();
/// let mut bases = 0;
/// while reader.read_into(&mut record)? {
///     bases += record.len();
/// }
/// assert_eq!(bases, 8);
/// # Ok::<(), fastx::Error>(())
/// ```
pub struct FastxReader<R: Read> {
    inner: R,
    buf: Vec<u8>,
    /// Cursor of the next unconsumed byte in `buf`.
    pos: usize,
    /// Number of valid bytes in `buf`.
    end: usize,
    eof: bool,
    format: Option<Format>,
    line: u64,
    quality_encoding: QualityEncoding,
    max_line_length: Option<usize>,
    max_record_length: Option<usize>,
    /// The previous record's description buffer, kept so that parsing a header
    /// does not allocate once per record. See `Sequence::set_header_reusing`.
    spare_description: String,
    /// Where `read_ref` joins a multi-line sequence, which cannot be borrowed
    /// from the buffer because the newlines are in the way. Untouched for
    /// single-line records, which is the common case.
    joined_seq: Vec<u8>,
    /// The same for a multi-line quality string.
    joined_quality: Vec<u8>,
    /// Ranges of the last scanned record's parts, in `buf`. Reused between
    /// records so that `read_ref` allocates nothing per record either.
    scan_header: (usize, usize),
    scan_seq: Vec<(usize, usize)>,
    scan_quality: Vec<(usize, usize)>,
    /// True when the record spans several lines and so had to be joined.
    joined: bool,
}

/// A cursor that walks lines in the buffered region without consuming them.
///
/// The scanners need to look ahead across a whole record and then either commit
/// or ask for more input, so they cannot advance the reader as they go.
struct Lines {
    pos: usize,
    end: usize,
    eof: bool,
    /// How many lines the cursor has walked, to advance the line counter on
    /// commit.
    consumed: u64,
}

impl Lines {
    fn new(pos: usize, end: usize, eof: bool) -> Lines {
        Lines {
            pos,
            end,
            eof,
            consumed: 0,
        }
    }

    /// True when the cursor has reached the end of a stream that has no more
    /// data coming.
    fn exhausted(&self) -> bool {
        self.pos == self.end && self.eof
    }

    /// The first byte of the next line, without consuming it.
    fn peek_first_byte(&self, buf: &[u8]) -> Option<u8> {
        buf.get(self.pos).copied().filter(|_| self.pos < self.end)
    }

    /// Consume one line, returning its bounds with any line terminator trimmed.
    ///
    /// `None` means the buffer holds no complete line — either more input is
    /// needed, or the stream has ended.
    fn next(&mut self, buf: &[u8]) -> Option<(usize, usize)> {
        match memchr::memchr(b'\n', &buf[self.pos..self.end]) {
            Some(offset) => {
                let newline = self.pos + offset;
                let start = self.pos;
                let mut stop = newline;
                if stop > start && buf[stop - 1] == b'\r' {
                    stop -= 1;
                }
                self.pos = newline + 1;
                self.consumed += 1;
                Some((start, stop))
            }
            // A final line without a trailing newline still counts, but only once
            // the reader knows nothing more is coming.
            None if self.eof && self.pos < self.end => {
                let start = self.pos;
                let mut stop = self.end;
                if stop > start && buf[stop - 1] == b'\r' {
                    stop -= 1;
                }
                self.pos = self.end;
                self.consumed += 1;
                Some((start, stop))
            }
            None => None,
        }
    }
}

/// End offset of `region[start..stop]` with a trailing `\r` removed.
fn trim_cr(region: &[u8], start: usize, stop: usize) -> usize {
    if stop > start && region[stop - 1] == b'\r' {
        stop - 1
    } else {
        stop
    }
}

/// Concatenate byte ranges of `buf` into `out`, which is cleared first.
fn join_ranges(ranges: &[(usize, usize)], buf: &[u8], out: &mut Vec<u8>) {
    out.clear();
    for &(start, stop) in ranges {
        out.extend_from_slice(&buf[start..stop]);
    }
}

impl<R: Read> FastxReader<R> {
    /// A reader that determines the format from the first record header.
    pub fn new(inner: R) -> FastxReader<R> {
        FastxReader::with_capacity(inner, DEFAULT_BUFFER_SIZE)
    }

    /// A reader with an explicit format, skipping auto-detection.
    pub fn with_format(inner: R, format: Format) -> FastxReader<R> {
        let mut reader = FastxReader::new(inner);
        reader.format = Some(format);
        reader
    }

    /// A reader with a custom buffer size (clamped to at least 4 KiB).
    pub fn with_capacity(inner: R, capacity: usize) -> FastxReader<R> {
        FastxReader {
            inner,
            buf: vec![0; capacity.max(MIN_BUFFER_SIZE)],
            pos: 0,
            end: 0,
            eof: false,
            format: None,
            line: 0,
            quality_encoding: QualityEncoding::Phred33,
            max_line_length: None,
            max_record_length: None,
            spare_description: String::new(),
            joined_seq: Vec::new(),
            joined_quality: Vec::new(),
            scan_header: (0, 0),
            scan_seq: Vec::new(),
            scan_quality: Vec::new(),
            joined: false,
        }
    }

    /// The format, once known. `None` before the first record has been read on
    /// an auto-detecting reader.
    pub fn format(&self) -> Option<Format> {
        self.format
    }

    /// The quality encoding of the *input*.
    ///
    /// Records themselves always come out as Phred+33: anything else is
    /// normalised while parsing, so downstream code never has to ask.
    pub fn quality_encoding(&self) -> QualityEncoding {
        self.quality_encoding
    }

    /// The 1-based number of the last line consumed; useful for error messages.
    pub fn line_number(&self) -> u64 {
        self.line
    }

    /// Unwrap the underlying reader, discarding any buffered bytes.
    pub fn into_inner(self) -> R {
        self.inner
    }

    /// Parse the next record into `record`, reusing its allocations.
    ///
    /// Returns `Ok(false)` at end of input. This is the allocation-free core of
    /// the reader; [`Iterator::next`] is a thin wrapper over it.
    pub fn read_into(&mut self, record: &mut Sequence) -> Result<bool> {
        let format = match self.format {
            Some(format) => {
                if !self.skip_blank_lines()? {
                    return Ok(false);
                }
                format
            }
            None => match self.detect_format()? {
                Some(format) => {
                    self.format = Some(format);
                    format
                }
                None => return Ok(false),
            },
        };
        // Take the previous description's buffer back before `clear` drops it, so
        // that one allocation serves the whole file instead of one per record.
        if let Some(mut previous) = record.description.take() {
            if previous.capacity() > self.spare_description.capacity() {
                previous.clear();
                self.spare_description = previous;
            }
        }
        record.clear();
        match format {
            Format::Fasta => self.read_fasta_into(record)?,
            Format::Fastq => self.read_fastq_into(record)?,
        }
        Ok(true)
    }

    /// Parse the next record into a fresh [`Sequence`].
    pub fn read_record(&mut self) -> Result<Option<Sequence>> {
        let mut record = Sequence::default();
        if self.read_into(&mut record)? {
            Ok(Some(record))
        } else {
            Ok(None)
        }
    }

    /// Borrowing iterator over the remaining records.
    ///
    /// Prefer this over consuming the reader when you still need it afterwards.
    pub fn records(&mut self) -> Records<'_, R> {
        Records { reader: self }
    }

    /// Run `f` on every remaining record, reusing a single buffer.
    ///
    /// This is the fastest way to consume a file and the one to reach for in
    /// pipelines: no per-record allocation, no `Result` per record to unwrap.
    ///
    /// ```
    /// use fastx::FastxReader;
    ///
    /// let data = b">a\nACGT\n>b\nGGCC\n";
    /// let mut total = 0;
    /// FastxReader::new(&data[..]).for_each_record(|r| { total += r.len(); Ok(()) })?;
    /// assert_eq!(total, 8);
    /// # Ok::<(), fastx::Error>(())
    /// ```
    pub fn for_each_record<F>(&mut self, mut f: F) -> Result<()>
    where
        F: FnMut(&Sequence) -> Result<()>,
    {
        let mut record = Sequence::default();
        while self.read_into(&mut record)? {
            f(&record)?;
        }
        Ok(())
    }

    /// Parse the next record without copying it, borrowing from the buffer.
    ///
    /// The returned record is valid until the next one is read, which is what
    /// lets this skip the copy that [`FastxReader::read_into`] makes. Nothing is
    /// copied for FASTQ or for single-line FASTA; a sequence spread over several
    /// lines has to be joined, so those records are assembled in a scratch buffer
    /// the reader owns and borrowed from there instead.
    ///
    /// Because each record borrows the reader, this cannot be an [`Iterator`] —
    /// use it in a `while let` loop, or reach for
    /// [`FastxReader::for_each_ref`].
    ///
    /// ```
    /// use fastx::FastxReader;
    ///
    /// let data = b"@r1\nACGT\n+\nIIII\n@r2\nTT\n+\n!!\n";
    /// let mut reader = FastxReader::new(&data[..]);
    /// let mut bases = 0;
    /// while let Some(record) = reader.read_ref()? {
    ///     bases += record.len();
    /// }
    /// assert_eq!(bases, 6);
    /// # Ok::<(), fastx::Error>(())
    /// ```
    pub fn read_ref(&mut self) -> Result<Option<SequenceRef<'_>>> {
        let format = match self.format {
            Some(format) => {
                if !self.skip_blank_lines()? {
                    return Ok(None);
                }
                format
            }
            None => match self.detect_format()? {
                Some(format) => {
                    self.format = Some(format);
                    format
                }
                None => return Ok(None),
            },
        };
        match format {
            Format::Fasta => self.read_fasta_ref(),
            Format::Fastq => self.read_fastq_ref(),
        }
        .map(Some)
    }

    /// Run `f` on every remaining record without copying any of them.
    ///
    /// The borrow-free way to use [`FastxReader::read_ref`]: the closure gets
    /// each record in turn, and none of them outlive the call.
    ///
    /// ```
    /// use fastx::FastxReader;
    ///
    /// let data = b">a\nACGT\n>b\nGGCC\n";
    /// let mut gc = 0;
    /// FastxReader::new(&data[..]).for_each_ref(|record| {
    ///     gc += record.base_counts().g + record.base_counts().c;
    ///     Ok(())
    /// })?;
    /// assert_eq!(gc, 6);
    /// # Ok::<(), fastx::Error>(())
    /// ```
    pub fn for_each_ref<F>(&mut self, mut f: F) -> Result<()>
    where
        F: FnMut(SequenceRef<'_>) -> Result<()>,
    {
        while let Some(record) = self.read_ref()? {
            f(record)?;
        }
        Ok(())
    }

    /// Count the remaining records without keeping them.
    pub fn count_records(&mut self) -> Result<u64> {
        let mut n = 0;
        let mut record = Sequence::default();
        while self.read_into(&mut record)? {
            n += 1;
        }
        Ok(n)
    }

    // ----- parsing ---------------------------------------------------------

    fn read_fasta_into(&mut self, record: &mut Sequence) -> Result<()> {
        let (start, end) = match self.read_line()? {
            Some(range) => range,
            None => {
                return Err(Error::parse(
                    self.line,
                    ParseError::UnexpectedEof {
                        expected: "a FASTA header",
                    },
                ))
            }
        };
        if self.buf[start] != b'>' && self.buf[start] != b';' {
            return Err(Error::parse(
                self.line,
                ParseError::ExpectedHeader {
                    found: self.buf[start],
                },
            ));
        }
        record.set_header_reusing(&self.buf[start + 1..end], &mut self.spare_description);
        if record.id.is_empty() {
            return Err(Error::parse(self.line, ParseError::EmptyId));
        }
        loop {
            match self.peek_byte()? {
                None | Some(b'>') => break,
                _ => {
                    let (start, end) = self.read_line()?.expect("peeked byte is available");
                    record.seq.extend_from_slice(&self.buf[start..end]);
                    self.check_record_limit(record.seq.len(), "sequence")?;
                }
            }
        }
        Ok(())
    }

    fn read_fastq_into(&mut self, record: &mut Sequence) -> Result<()> {
        let (start, end) = match self.read_line()? {
            Some(range) => range,
            None => {
                return Err(Error::parse(
                    self.line,
                    ParseError::UnexpectedEof {
                        expected: "a FASTQ header",
                    },
                ))
            }
        };
        if self.buf[start] != b'@' {
            return Err(Error::parse(
                self.line,
                ParseError::ExpectedHeader {
                    found: self.buf[start],
                },
            ));
        }
        record.set_header_reusing(&self.buf[start + 1..end], &mut self.spare_description);
        if record.id.is_empty() {
            return Err(Error::parse(self.line, ParseError::EmptyId));
        }

        // Sequence lines, up to the '+' separator. Multi-line FASTQ is rare but
        // legal, and a '+' can never start a sequence line.
        loop {
            match self.peek_byte()? {
                None => {
                    return Err(Error::parse(
                        self.line,
                        ParseError::UnexpectedEof {
                            expected: "a FASTQ '+' separator",
                        },
                    ))
                }
                Some(b'+') => {
                    self.read_line()?;
                    break;
                }
                _ => {
                    let (start, end) = self.read_line()?.expect("peeked byte is available");
                    record.seq.extend_from_slice(&self.buf[start..end]);
                    self.check_record_limit(record.seq.len(), "sequence")?;
                }
            }
        }

        // Quality lines. Because a quality character may itself be '@', the only
        // safe terminator is having collected as many scores as bases.
        let quality = record.quality.get_or_insert_with(Vec::new);
        while quality.len() < record.seq.len() {
            match self.read_line()? {
                Some((start, end)) => quality.extend_from_slice(&self.buf[start..end]),
                None => {
                    return Err(Error::LengthMismatch {
                        id: record.id.clone(),
                        seq: record.seq.len(),
                        quality: quality.len(),
                    })
                }
            }
        }
        if quality.len() != record.seq.len() {
            return Err(Error::LengthMismatch {
                id: record.id.clone(),
                seq: record.seq.len(),
                quality: quality.len(),
            });
        }
        // Normalise to Phred+33 so that a `Sequence` has exactly one encoding,
        // whatever the file used.
        if self.quality_encoding != QualityEncoding::Phred33 {
            let from = self.quality_encoding.offset();
            for c in quality.iter_mut() {
                *c = qual::encode(qual::score(*c, from), qual::PHRED33);
            }
        }
        Ok(())
    }

    /// Fail rather than buffer without bound when a line exceeds its limit.
    fn check_line_limit(&self, length: usize) -> Result<()> {
        match self.max_line_length {
            Some(limit) if length > limit => Err(Error::TooLarge {
                line: self.line + 1,
                what: "line",
                limit,
            }),
            _ => Ok(()),
        }
    }

    /// Fail rather than grow without bound when a record exceeds its limit.
    fn check_record_limit(&self, length: usize, what: &'static str) -> Result<()> {
        match self.max_record_length {
            Some(limit) if length > limit => Err(Error::TooLarge {
                line: self.line,
                what,
                limit,
            }),
            _ => Ok(()),
        }
    }

    // ----- borrowed parsing --------------------------------------------------

    /// Parse a FASTQ record already known to be buffered, or borrow one.
    ///
    /// Every offset is computed after the final refill, because refilling
    /// compacts the buffer and moves the bytes the offsets refer to.
    fn read_fastq_ref(&mut self) -> Result<SequenceRef<'_>> {
        loop {
            if self.scan_fastq_quick()? || self.scan_fastq()? {
                break;
            }
            if self.eof {
                return Err(Error::parse(
                    self.line,
                    ParseError::UnexpectedEof {
                        expected: "a complete FASTQ record",
                    },
                ));
            }
            self.check_record_limit(self.end - self.pos, "record")?;
            self.refill()?;
        }

        if self.joined {
            join_ranges(&self.scan_seq, &self.buf, &mut self.joined_seq);
            join_ranges(&self.scan_quality, &self.buf, &mut self.joined_quality);
        }
        let (seq, quality) = if self.joined {
            (&self.joined_seq[..], &self.joined_quality[..])
        } else {
            let seq = self.scan_seq[0];
            let quality = self.scan_quality[0];
            (&self.buf[seq.0..seq.1], &self.buf[quality.0..quality.1])
        };
        let (id, description) =
            crate::record::split_header(&self.buf[self.scan_header.0..self.scan_header.1]);
        if id.is_empty() {
            return Err(Error::parse(self.line, ParseError::EmptyId));
        }
        Ok(SequenceRef::new(id, description, seq, Some(quality)))
    }

    /// The FASTA equivalent: a record runs until the next header or end of input.
    fn read_fasta_ref(&mut self) -> Result<SequenceRef<'_>> {
        loop {
            if self.scan_fasta()? {
                break;
            }
            self.check_record_limit(self.end - self.pos, "record")?;
            self.refill()?;
        }

        if self.joined {
            join_ranges(&self.scan_seq, &self.buf, &mut self.joined_seq);
        }
        let seq = if self.joined {
            &self.joined_seq[..]
        } else {
            match self.scan_seq.first() {
                Some(&(start, stop)) => &self.buf[start..stop],
                None => &[][..],
            }
        };
        let (id, description) =
            crate::record::split_header(&self.buf[self.scan_header.0..self.scan_header.1]);
        if id.is_empty() {
            return Err(Error::parse(self.line, ParseError::EmptyId));
        }
        Ok(SequenceRef::new(id, description, seq, None))
    }

    /// The shape almost every FASTQ record has: four lines, one each.
    ///
    /// Worth a path of its own because the general scanner pays for generality it
    /// does not need here — a `Vec` push and clear per part, a second pass over
    /// those `Vec`s to check limits and a third to sum lengths. This finds the
    /// four newlines with a single `memchr_iter`, so one SIMD setup covers the
    /// whole record instead of four, and keeps the offsets in locals.
    ///
    /// Returns `false` when the record is not this shape, or when more input is
    /// needed; either way the general scanner then runs and decides which it was.
    fn scan_fastq_quick(&mut self) -> Result<bool> {
        let region = &self.buf[self.pos..self.end];
        let mut newlines = memchr::memchr_iter(b'\n', region);
        let (Some(a), Some(b), Some(c), Some(d)) = (
            newlines.next(),
            newlines.next(),
            newlines.next(),
            newlines.next(),
        ) else {
            return Ok(false);
        };

        // Bounds within the buffer, with the terminators trimmed.
        let base = self.pos;
        let header = (base, trim_cr(region, 0, a) + base);
        let seq = (base + a + 1, trim_cr(region, a + 1, b) + base);
        let plus = (base + b + 1, trim_cr(region, b + 1, c) + base);
        let quality = (base + c + 1, trim_cr(region, c + 1, d) + base);

        // Anything unusual — a missing separator, or a record whose parts run
        // over several lines — goes to the general scanner rather than being
        // half-handled here.
        if self.buf[header.0] != b'@'
            || self.buf[plus.0] != b'+'
            || quality.1 - quality.0 != seq.1 - seq.0
        {
            return Ok(false);
        }

        self.check_line_limit(seq.1 - seq.0)?;
        self.check_line_limit(quality.1 - quality.0)?;
        self.check_record_limit(seq.1 - seq.0, "sequence")?;

        self.scan_header = (header.0 + 1, header.1);
        self.scan_seq.clear();
        self.scan_seq.push(seq);
        self.scan_quality.clear();
        self.scan_quality.push(quality);
        self.joined = false;
        self.pos = base + d + 1;
        self.line += 4;
        Ok(true)
    }

    /// Look for a complete FASTQ record in the buffered region without consuming
    /// it. Returns false when more input is needed.
    ///
    /// On success `pos` and `line` have advanced past the record and the scan
    /// fields describe where its parts are.
    fn scan_fastq(&mut self) -> Result<bool> {
        let mut lines = Lines::new(self.pos, self.end, self.eof);

        let header = match lines.next(&self.buf) {
            Some(range) => range,
            None => return Ok(false),
        };
        if self.buf[header.0] != b'@' {
            return Err(Error::parse(
                self.line + 1,
                ParseError::ExpectedHeader {
                    found: self.buf[header.0],
                },
            ));
        }

        // Sequence lines, up to the '+' separator: '+' can never start one.
        self.scan_seq.clear();
        let mut seq_len = 0;
        loop {
            match lines.peek_first_byte(&self.buf) {
                None => return Ok(false),
                Some(b'+') => {
                    lines.next(&self.buf);
                    break;
                }
                Some(_) => match lines.next(&self.buf) {
                    None => return Ok(false),
                    Some(range) => {
                        seq_len += range.1 - range.0;
                        self.scan_seq.push(range);
                    }
                },
            }
        }

        // Quality lines, until as many scores as bases: a quality character may
        // itself be '@', so length is the only safe terminator.
        self.scan_quality.clear();
        let mut quality_len = 0;
        while quality_len < seq_len {
            match lines.next(&self.buf) {
                None => return Ok(false),
                Some(range) => {
                    quality_len += range.1 - range.0;
                    self.scan_quality.push(range);
                }
            }
        }
        if quality_len != seq_len {
            let (id, _) = crate::record::split_header(&self.buf[header.0 + 1..header.1]);
            return Err(Error::LengthMismatch {
                id: String::from_utf8_lossy(id).into_owned(),
                seq: seq_len,
                quality: quality_len,
            });
        }

        self.finish_scan(header, &lines)
    }

    /// The FASTA equivalent: consume lines until the next header or end of input.
    fn scan_fasta(&mut self) -> Result<bool> {
        let mut lines = Lines::new(self.pos, self.end, self.eof);

        let header = match lines.next(&self.buf) {
            Some(range) => range,
            None => return Ok(false),
        };
        if self.buf[header.0] != b'>' && self.buf[header.0] != b';' {
            return Err(Error::parse(
                self.line + 1,
                ParseError::ExpectedHeader {
                    found: self.buf[header.0],
                },
            ));
        }

        self.scan_seq.clear();
        self.scan_quality.clear();
        loop {
            match lines.peek_first_byte(&self.buf) {
                // A header ends the record; end of input does too, but only once
                // the reader knows there is no more coming.
                Some(b'>') => break,
                None if lines.exhausted() => break,
                None => return Ok(false),
                Some(_) => match lines.next(&self.buf) {
                    None => return Ok(false),
                    Some(range) => self.scan_seq.push(range),
                },
            }
        }

        self.finish_scan(header, &lines)
    }

    /// Commit a successful scan: record where the parts are and consume the
    /// input the scanner walked over.
    fn finish_scan(&mut self, header: (usize, usize), lines: &Lines) -> Result<bool> {
        for &(start, stop) in self.scan_seq.iter().chain(self.scan_quality.iter()) {
            self.check_line_limit(stop - start)?;
        }
        let length: usize = self.scan_seq.iter().map(|&(a, b)| b - a).sum();
        self.check_record_limit(length, "sequence")?;

        self.scan_header = (header.0 + 1, header.1);
        // More than one line for either part means the bytes are not contiguous,
        // so they have to be joined into a scratch buffer rather than borrowed.
        self.joined = self.scan_seq.len() > 1 || self.scan_quality.len() > 1;
        self.pos = lines.pos;
        self.line += lines.consumed;
        Ok(true)
    }

    /// Advance past blank lines. Returns false at end of input.
    fn skip_blank_lines(&mut self) -> Result<bool> {
        loop {
            match self.peek_byte()? {
                None => return Ok(false),
                Some(b'\n') => {
                    self.pos += 1;
                    self.line += 1;
                }
                Some(b'\r') => self.pos += 1,
                Some(_) => return Ok(true),
            }
        }
    }

    /// Sniff the format from the first meaningful byte without consuming it.
    fn detect_format(&mut self) -> Result<Option<Format>> {
        if !self.skip_blank_lines()? {
            return Ok(None);
        }
        let byte = self.buf[self.pos];
        match Format::from_first_byte(byte) {
            Some(format) => Ok(Some(format)),
            None => Err(Error::parse(
                self.line + 1,
                ParseError::ExpectedHeader { found: byte },
            )),
        }
    }

    // ----- buffer management ------------------------------------------------

    /// Consume one line, returning its bounds in `self.buf` without the line
    /// terminator. Returns `None` only at end of input.
    fn read_line(&mut self) -> Result<Option<(usize, usize)>> {
        let mut search_from = self.pos;
        loop {
            if let Some(offset) = memchr::memchr(b'\n', &self.buf[search_from..self.end]) {
                let newline = search_from + offset;
                let start = self.pos;
                let mut stop = newline;
                if stop > start && self.buf[stop - 1] == b'\r' {
                    stop -= 1;
                }
                self.check_line_limit(stop - start)?;
                self.pos = newline + 1;
                self.line += 1;
                return Ok(Some((start, stop)));
            }
            if self.eof {
                if self.pos == self.end {
                    return Ok(None);
                }
                // Final line without a trailing newline.
                let start = self.pos;
                let mut stop = self.end;
                if stop > start && self.buf[stop - 1] == b'\r' {
                    stop -= 1;
                }
                self.check_line_limit(stop - start)?;
                self.pos = self.end;
                self.line += 1;
                return Ok(Some((start, stop)));
            }
            // No newline yet, so everything buffered belongs to the current line.
            // Checking here as well is what stops the buffer growing without
            // bound on input that never supplies a newline at all.
            self.check_line_limit(self.end - self.pos)?;
            let previous_end = self.end;
            let shift = self.refill()?;
            search_from = previous_end - shift;
        }
    }

    /// Ensure at least one byte is buffered and return it without consuming.
    fn peek_byte(&mut self) -> Result<Option<u8>> {
        while self.pos == self.end && !self.eof {
            self.refill()?;
        }
        if self.pos == self.end {
            Ok(None)
        } else {
            Ok(Some(self.buf[self.pos]))
        }
    }

    /// Move unconsumed bytes to the front, grow if the buffer is full, then read.
    /// Returns how far indices into `buf` shifted left.
    fn refill(&mut self) -> Result<usize> {
        let mut shift = 0;
        if self.pos > 0 {
            self.buf.copy_within(self.pos..self.end, 0);
            shift = self.pos;
            self.end -= self.pos;
            self.pos = 0;
        }
        if self.end == self.buf.len() {
            // A single line longer than the buffer: double it.
            let grown = self.buf.len().saturating_mul(2).max(MIN_BUFFER_SIZE);
            self.buf.resize(grown, 0);
        }
        loop {
            match self.inner.read(&mut self.buf[self.end..]) {
                Ok(0) => {
                    self.eof = true;
                    break;
                }
                Ok(n) => {
                    self.end += n;
                    break;
                }
                Err(e) if e.kind() == io::ErrorKind::Interrupted => continue,
                Err(e) => return Err(Error::Io(e)),
            }
        }
        Ok(shift)
    }
}

impl<R: Read> fmt::Debug for FastxReader<R> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.debug_struct("FastxReader")
            .field("format", &self.format)
            .field("buffer_size", &self.buf.len())
            .field("buffered", &(self.end - self.pos))
            .field("line", &self.line)
            .field("eof", &self.eof)
            .finish_non_exhaustive()
    }
}

impl<R: Read> Iterator for FastxReader<R> {
    type Item = Result<Sequence>;

    fn next(&mut self) -> Option<Self::Item> {
        match self.read_record() {
            Ok(Some(record)) => Some(Ok(record)),
            Ok(None) => None,
            Err(e) => Some(Err(e)),
        }
    }
}

/// Borrowing iterator returned by [`FastxReader::records`].
pub struct Records<'a, R: Read> {
    reader: &'a mut FastxReader<R>,
}

impl<R: Read> Iterator for Records<'_, R> {
    type Item = Result<Sequence>;

    fn next(&mut self) -> Option<Self::Item> {
        match self.reader.read_record() {
            Ok(Some(record)) => Some(Ok(record)),
            Ok(None) => None,
            Err(e) => Some(Err(e)),
        }
    }
}

/// Configuration for a [`FastxReader`].
///
/// ```
/// use fastx::{Format, ReaderBuilder};
///
/// let data = b">a\nACGT\n";
/// let mut reader = ReaderBuilder::new()
///     .format(Format::Fasta)
///     .buffer_size(64 * 1024)
///     .build(&data[..]);
/// assert_eq!(reader.next().unwrap()?.id, "a");
/// # Ok::<(), fastx::Error>(())
/// ```
#[derive(Debug, Clone)]
pub struct ReaderBuilder {
    format: Option<Format>,
    buffer_size: usize,
    quality_encoding: QualityEncoding,
    max_line_length: Option<usize>,
    max_record_length: Option<usize>,
}

impl Default for ReaderBuilder {
    fn default() -> Self {
        ReaderBuilder {
            format: None,
            buffer_size: DEFAULT_BUFFER_SIZE,
            quality_encoding: QualityEncoding::Phred33,
            max_line_length: None,
            max_record_length: None,
        }
    }
}

impl ReaderBuilder {
    /// A builder with default settings.
    pub fn new() -> ReaderBuilder {
        ReaderBuilder::default()
    }

    /// Force a format instead of detecting it.
    pub fn format(mut self, format: Format) -> Self {
        self.format = Some(format);
        self
    }

    /// Size of the internal read buffer, in bytes.
    pub fn buffer_size(mut self, bytes: usize) -> Self {
        self.buffer_size = bytes;
        self
    }

    /// The quality encoding of the input.
    ///
    /// Old Illumina 1.3–1.7 files are Phred+64. Set this and the reader will
    /// convert quality strings to Phred+33 as it parses, so every [`Sequence`]
    /// this crate produces uses one encoding.
    ///
    /// ```
    /// use fastx::{qual::QualityEncoding, ReaderBuilder};
    ///
    /// // 'h' is Q40 in Phred+64.
    /// let data = b"@old\nACGT\n+\nhhhh\n";
    /// let record = ReaderBuilder::new()
    ///     .quality_encoding(QualityEncoding::Phred64)
    ///     .build(&data[..])
    ///     .read_record()?
    ///     .unwrap();
    ///
    /// assert_eq!(record.quality.as_deref(), Some(&b"IIII"[..])); // now Phred+33
    /// assert_eq!(record.mean_quality(), Some(40.0));
    /// # Ok::<(), fastx::Error>(())
    /// ```
    pub fn quality_encoding(mut self, encoding: QualityEncoding) -> Self {
        self.quality_encoding = encoding;
        self
    }

    /// Refuse lines longer than `bytes` instead of growing the buffer.
    ///
    /// Unlimited by default, because a chromosome legitimately arrives on a
    /// single line. Set it when reading files you do not control: without a
    /// limit, one unterminated line can grow the buffer until the process is
    /// killed.
    pub fn max_line_length(mut self, bytes: usize) -> Self {
        self.max_line_length = Some(bytes);
        self
    }

    /// Refuse records whose sequence exceeds `bytes`.
    ///
    /// Unlimited by default. This bounds a record assembled from many short
    /// lines, which `max_line_length` alone does not catch.
    pub fn max_record_length(mut self, bytes: usize) -> Self {
        self.max_record_length = Some(bytes);
        self
    }

    /// Build a reader around any [`Read`].
    pub fn build<R: Read>(&self, inner: R) -> FastxReader<R> {
        let mut reader = FastxReader::with_capacity(inner, self.buffer_size);
        reader.format = self.format;
        reader.quality_encoding = self.quality_encoding;
        reader.max_line_length = self.max_line_length;
        reader.max_record_length = self.max_record_length;
        reader
    }

    /// Open a path, transparently decompressing gzip and inferring the format.
    pub fn open<P: AsRef<Path>>(&self, path: P) -> Result<FastxReader<Box<dyn Read + Send>>> {
        let path = path.as_ref();
        let mut builder = self.clone();
        if builder.format.is_none() {
            builder.format = Format::from_path(path);
        }
        Ok(builder.build(open_reader(path)?))
    }
}

/// The boxed reader type produced by [`open`] and [`from_stdin`].
pub type BoxedReader = FastxReader<Box<dyn Read + Send>>;

/// Open a FASTA/FASTQ file, transparently handling gzip.
///
/// The format is taken from the extension when recognisable and otherwise from
/// the first byte of the (decompressed) stream. gzip is detected from the file's
/// magic bytes, so a compressed file without a `.gz` suffix works as well.
///
/// Requires the `gzip` feature for compressed input.
pub fn open<P: AsRef<Path>>(path: P) -> Result<BoxedReader> {
    ReaderBuilder::default().open(path)
}

/// Read records from standard input (gzip is detected from the magic bytes).
pub fn from_stdin() -> Result<BoxedReader> {
    let stream = decompress(Box::new(io::stdin()))?;
    Ok(FastxReader::new(stream))
}

/// Wrap a file in a decompressing reader when needed.
fn open_reader(path: &Path) -> Result<Box<dyn Read + Send>> {
    let file = File::open(path)
        .map_err(|e| Error::Io(io::Error::new(e.kind(), format!("{}: {e}", path.display()))))?;
    decompress(Box::new(BufReader::with_capacity(64 * 1024, file)))
}

/// Peek at the magic bytes and wrap the stream in a decompressor if needed.
///
/// The bytes read for sniffing are put back in front of the stream, so this
/// works on a pipe as well as on a file — nothing is seeked.
fn decompress(mut stream: Box<dyn Read + Send>) -> Result<Box<dyn Read + Send>> {
    // Enough to cover a whole BGZF block header, so that BGZF can be told from
    // plain gzip here and decompressed across cores. zstd needs four bytes,
    // gzip two.
    let mut magic = [0u8; HEADER_LEN_PROBE];
    let mut filled = 0;
    while filled < magic.len() {
        match stream.read(&mut magic[filled..]) {
            Ok(0) => break,
            Ok(n) => filled += n,
            Err(e) if e.kind() == io::ErrorKind::Interrupted => continue,
            Err(e) => return Err(Error::Io(e)),
        }
    }
    let head = io::Cursor::new(magic[..filled].to_vec());
    let rejoined = head.chain(stream);
    match Compression::from_magic(&magic[..filled]) {
        Compression::None => Ok(Box::new(rejoined)),
        // BGZF is gzip, and reading it sequentially needs no special handling —
        // MultiGzDecoder walks the members. Only random access cares, and that
        // goes through `crate::bgzf`.
        Compression::Gzip | Compression::Bgzf => gunzip(rejoined, &magic[..filled]),
        Compression::Zstd => unzstd(rejoined),
    }
}

/// Bytes sniffed before deciding how to decompress: a BGZF header is 18 bytes
/// with the standard single `BC` extra subfield.
const HEADER_LEN_PROBE: usize = 18;

/// Decode gzip, using every core when the input is BGZF.
///
/// BGZF blocks are independent, so a whole-file pass over one can inflate across
/// cores. Plain gzip is a single deflate stream and cannot: there, and without
/// the `parallel` feature, `MultiGzDecoder` walks the members one at a time.
#[cfg(feature = "gzip")]
fn gunzip<R: Read + Send + 'static>(stream: R, head: &[u8]) -> Result<Box<dyn Read + Send>> {
    #[cfg(feature = "parallel")]
    if crate::bgzf::is_bgzf(head) {
        return Ok(Box::new(crate::bgzf::ParallelBgzfReader::new(stream)));
    }
    let _ = head;
    Ok(Box::new(flate2::read::MultiGzDecoder::new(stream)))
}

#[cfg(not(feature = "gzip"))]
fn gunzip<R: Read + Send + 'static>(_stream: R, _head: &[u8]) -> Result<Box<dyn Read + Send>> {
    Err(Error::FeatureDisabled("gzip"))
}

/// Decode a Zstandard stream, including one made of several frames.
#[cfg(feature = "zstd")]
fn unzstd<R: Read + Send + 'static>(stream: R) -> Result<Box<dyn Read + Send>> {
    Ok(Box::new(zstd::stream::read::Decoder::new(stream)?))
}

#[cfg(not(feature = "zstd"))]
fn unzstd<R: Read + Send + 'static>(_stream: R) -> Result<Box<dyn Read + Send>> {
    Err(Error::FeatureDisabled("zstd"))
}

#[cfg(test)]
mod tests {
    use super::*;

    fn ids(data: &[u8]) -> Vec<String> {
        FastxReader::new(data).map(|r| r.unwrap().id).collect()
    }

    #[test]
    fn reads_simple_fasta() {
        let data = b">a desc here\nACGT\n>b\nTTTT\nGGGG\n";
        let records: Vec<_> = FastxReader::new(&data[..])
            .collect::<Result<Vec<_>>>()
            .unwrap();
        assert_eq!(records.len(), 2);
        assert_eq!(records[0].id, "a");
        assert_eq!(records[0].description.as_deref(), Some("desc here"));
        assert_eq!(records[0].seq, b"ACGT");
        assert_eq!(records[1].seq, b"TTTTGGGG");
        assert!(records[1].quality.is_none());
    }

    #[test]
    fn reads_simple_fastq() {
        let data = b"@a\nACGT\n+\nIIII\n@b desc\nTT\n+b desc\n!!\n";
        let records: Vec<_> = FastxReader::new(&data[..])
            .collect::<Result<Vec<_>>>()
            .unwrap();
        assert_eq!(records.len(), 2);
        assert_eq!(records[0].quality.as_deref(), Some(&b"IIII"[..]));
        assert_eq!(records[1].id, "b");
        assert_eq!(records[1].description.as_deref(), Some("desc"));
        assert_eq!(records[1].quality.as_deref(), Some(&b"!!"[..]));
    }

    #[test]
    fn detects_format() {
        let mut reader = FastxReader::new(&b">a\nAC\n"[..]);
        assert_eq!(reader.format(), None);
        reader.next().unwrap().unwrap();
        assert_eq!(reader.format(), Some(Format::Fasta));

        let mut reader = FastxReader::new(&b"@a\nAC\n+\nII\n"[..]);
        reader.next().unwrap().unwrap();
        assert_eq!(reader.format(), Some(Format::Fastq));
    }

    #[test]
    fn handles_crlf_and_missing_final_newline() {
        let data = b">a\r\nACGT\r\nAC\r\n>b\r\nTT";
        let records: Vec<_> = FastxReader::new(&data[..])
            .collect::<Result<Vec<_>>>()
            .unwrap();
        assert_eq!(records[0].seq, b"ACGTAC");
        assert_eq!(records[1].seq, b"TT");
    }

    #[test]
    fn handles_blank_lines_between_records() {
        let data = b"\n\n>a\nACGT\n\n\n>b\nTT\n\n";
        assert_eq!(ids(&data[..]), ["a", "b"]);
        // The blank line inside record `a` must not become part of the sequence.
        let records: Vec<_> = FastxReader::new(&data[..])
            .collect::<Result<Vec<_>>>()
            .unwrap();
        assert_eq!(records[0].seq, b"ACGT");
    }

    #[test]
    fn handles_empty_input() {
        assert_eq!(FastxReader::new(&b""[..]).count(), 0);
        assert_eq!(FastxReader::new(&b"\n\n\n"[..]).count(), 0);
    }

    #[test]
    fn multi_line_fastq() {
        let data = b"@a\nACGT\nACGT\n+\nIIII\nJJJJ\n@b\nTT\n+\n!!\n";
        let records: Vec<_> = FastxReader::new(&data[..])
            .collect::<Result<Vec<_>>>()
            .unwrap();
        assert_eq!(records[0].seq, b"ACGTACGT");
        assert_eq!(records[0].quality.as_deref(), Some(&b"IIIIJJJJ"[..]));
        assert_eq!(records[1].id, "b");
    }

    #[test]
    fn quality_starting_with_at_sign() {
        // '@' is a legal quality character (Q31 in Phred+33).
        let data = b"@a\nACGT\n+\n@@@@\n@b\nTTTT\n+\nIIII\n";
        let records: Vec<_> = FastxReader::new(&data[..])
            .collect::<Result<Vec<_>>>()
            .unwrap();
        assert_eq!(records.len(), 2);
        assert_eq!(records[0].quality.as_deref(), Some(&b"@@@@"[..]));
        assert_eq!(records[1].id, "b");
    }

    #[test]
    fn tiny_buffer_still_parses() {
        // Force many refills and a line longer than the initial buffer.
        let long = "A".repeat(50_000);
        let data = format!(">a\n{long}\n>b\nACGT\n");
        let mut reader = FastxReader::with_capacity(data.as_bytes(), 1);
        let records: Vec<_> = reader.records().collect::<Result<Vec<_>>>().unwrap();
        assert_eq!(records.len(), 2);
        assert_eq!(records[0].seq.len(), 50_000);
        assert_eq!(records[1].seq, b"ACGT");
    }

    #[test]
    fn description_buffer_is_reused_without_leaking_between_records() {
        // The reader keeps the previous description's buffer to avoid allocating
        // per record. The case that breaks a naive version of that is a record
        // with a description followed by one without: the buffer has to be
        // reclaimed while the field still reads as `None`.
        let data = b"@a first one\nAC\n+\nII\n\
                     @b\nGT\n+\nII\n\
                     @c third\nTT\n+\nII\n\
                     @d\nCC\n+\nII\n\
                     @e a much longer description than any before it\nGG\n+\nII\n";
        let mut reader = FastxReader::new(&data[..]);
        let mut record = Sequence::default();
        let mut seen = Vec::new();
        while reader.read_into(&mut record).unwrap() {
            seen.push((record.id.clone(), record.description.clone()));
        }
        assert_eq!(
            seen,
            vec![
                ("a".to_string(), Some("first one".to_string())),
                ("b".to_string(), None),
                ("c".to_string(), Some("third".to_string())),
                ("d".to_string(), None),
                (
                    "e".to_string(),
                    Some("a much longer description than any before it".to_string())
                ),
            ]
        );

        // Reading into a record the caller already filled must not blend the two.
        let mut record = Sequence::fasta("old", b"AAAA").with_description("stale description");
        let data = b"@new\nAC\n+\nII\n";
        assert!(FastxReader::new(&data[..]).read_into(&mut record).unwrap());
        assert_eq!(record.id, "new");
        assert_eq!(record.description, None);
        assert_eq!(record.seq, b"AC");
    }

    #[test]
    fn read_into_reuses_allocations() {
        let data = b">a\nACGT\n>b\nTT\n";
        let mut reader = FastxReader::new(&data[..]);
        let mut record = Sequence::default();
        assert!(reader.read_into(&mut record).unwrap());
        assert_eq!(record.id, "a");
        assert!(reader.read_into(&mut record).unwrap());
        assert_eq!(record.id, "b");
        assert_eq!(record.seq, b"TT");
        assert!(!reader.read_into(&mut record).unwrap());
    }

    #[test]
    fn empty_fasta_record_is_allowed() {
        let data = b">a\n>b\nACGT\n";
        let records: Vec<_> = FastxReader::new(&data[..])
            .collect::<Result<Vec<_>>>()
            .unwrap();
        assert_eq!(records[0].seq, b"");
        assert_eq!(records[1].seq, b"ACGT");
    }

    #[test]
    fn rejects_garbage() {
        let err = FastxReader::new(&b"not a sequence file\n"[..])
            .next()
            .unwrap()
            .unwrap_err();
        assert!(matches!(
            err,
            Error::Parse {
                kind: ParseError::ExpectedHeader { found: b'n' },
                ..
            }
        ));
    }

    #[test]
    fn rejects_truncated_fastq() {
        let err = FastxReader::new(&b"@a\nACGT\n"[..])
            .next()
            .unwrap()
            .unwrap_err();
        assert!(matches!(
            err,
            Error::Parse {
                kind: ParseError::UnexpectedEof { .. },
                ..
            }
        ));

        let err = FastxReader::new(&b"@a\nACGT\n+\nII\n"[..])
            .next()
            .unwrap()
            .unwrap_err();
        assert!(matches!(
            err,
            Error::LengthMismatch {
                seq: 4,
                quality: 2,
                ..
            }
        ));
    }

    #[test]
    fn rejects_empty_id() {
        let err = FastxReader::new(&b">\nACGT\n"[..])
            .next()
            .unwrap()
            .unwrap_err();
        assert!(matches!(
            err,
            Error::Parse {
                kind: ParseError::EmptyId,
                ..
            }
        ));
    }

    #[test]
    fn reports_line_numbers() {
        let data = b">a\nACGT\n>b\nACGT\nnope";
        let mut reader = FastxReader::with_format(&data[..], Format::Fasta);
        reader.next().unwrap().unwrap();
        assert_eq!(reader.line_number(), 2);
    }

    #[test]
    fn phred64_input_is_normalised_to_phred33() {
        // 'h' is Q40 in Phred+64, 'B' is Q2.
        let data = b"@old\nACGT\n+\nhhhB\n";

        let record = ReaderBuilder::new()
            .quality_encoding(QualityEncoding::Phred64)
            .build(&data[..])
            .read_record()
            .unwrap()
            .unwrap();
        assert_eq!(record.quality.as_deref(), Some(&b"III#"[..]));
        assert_eq!(record.quality_scores().unwrap(), vec![40, 40, 40, 2]);

        // Without the setting the same bytes are read as Phred+33 verbatim.
        let record = FastxReader::new(&data[..]).read_record().unwrap().unwrap();
        assert_eq!(record.quality.as_deref(), Some(&b"hhhB"[..]));
    }

    #[test]
    fn line_length_limit_is_enforced() {
        let long = format!(">a\n{}\n", "A".repeat(10_000));
        let err = ReaderBuilder::new()
            .max_line_length(1_000)
            .build(long.as_bytes())
            .read_record()
            .unwrap_err();
        assert!(
            matches!(
                err,
                Error::TooLarge {
                    what: "line",
                    limit: 1_000,
                    ..
                }
            ),
            "{err}"
        );

        // Under the limit it parses normally.
        let record = ReaderBuilder::new()
            .max_line_length(1_000_000)
            .build(long.as_bytes())
            .read_record()
            .unwrap()
            .unwrap();
        assert_eq!(record.seq.len(), 10_000);
    }

    #[test]
    fn record_length_limit_catches_many_short_lines() {
        // 200 lines of 50 bases: every line is small, the record is not.
        let mut data = String::from(">a\n");
        for _ in 0..200 {
            data.push_str(&"A".repeat(50));
            data.push('\n');
        }
        let err = ReaderBuilder::new()
            .max_line_length(1_000)
            .max_record_length(5_000)
            .build(data.as_bytes())
            .read_record()
            .unwrap_err();
        assert!(
            matches!(
                err,
                Error::TooLarge {
                    what: "sequence",
                    limit: 5_000,
                    ..
                }
            ),
            "{err}"
        );
    }

    #[test]
    fn limits_are_unlimited_by_default() {
        // A single line far larger than the buffer must still be accepted.
        let long = format!(">chrom\n{}\n", "ACGT".repeat(50_000));
        let record = FastxReader::with_capacity(long.as_bytes(), 4096)
            .read_record()
            .unwrap()
            .unwrap();
        assert_eq!(record.seq.len(), 200_000);
    }

    #[test]
    fn forced_format_reads_fasta_as_written() {
        let data = b">a\nACGT\n";
        let mut reader = FastxReader::with_format(&data[..], Format::Fasta);
        assert_eq!(reader.format(), Some(Format::Fasta));
        assert_eq!(reader.next().unwrap().unwrap().seq, b"ACGT");
    }
}