rasusa 5.1.0

use crate::format::OutputEncoding;
use crate::source::RecordSource;
use crate::subsampler::seeded_rng;
use needletail::errors::ParseErrorKind::EmptyFile;
use rand::prelude::*;
use std::fs::File;
use std::io::{BufWriter, Write};
use std::path::{Path, PathBuf};
use thiserror::Error;

/// A collection of custom errors relating to the working with files for this package.
#[derive(Error, Debug)]
pub enum FastxError {
    /// Indicates that the specified input file could not be opened/read.
    #[error("Read error")]
    ReadError {
        source: needletail::errors::ParseError,
    },

    /// Indicates that a sequence record could not be parsed.
    #[error("Failed to parse record")]
    ParseError {
        source: needletail::errors::ParseError,
    },

    /// Indicates that the specified output file could not be created.
    #[error("Output file could not be created")]
    CreateError { source: std::io::Error },

    /// Indicates and error trying to create the compressor
    #[error(transparent)]
    CompressOutputError(#[from] niffler::Error),

    /// Indicates that some indices we expected to find in the input file weren't found.
    #[error("Some expected indices were not in the input file")]
    IndicesNotFound,

    /// Indicates that writing to the output file failed.
    #[error("Could not write to output file")]
    WriteError { source: anyhow::Error },

    /// Indicates an error when reading an unaligned alignment file (SAM/BAM/CRAM).
    #[error("Alignment read error: {source}")]
    AlignmentReadError { source: std::io::Error },

    /// Indicates that a paired one-pass run's two input files disagreed on read count, detected
    /// as soon as one reader ran out of records before the other.
    #[error(
        "Paired FASTA/Q inputs have different numbers of reads: after {matched} matched read \
         pair(s), the {longer} input had more. Paired Illumina files are assumed to have the \
         same number of reads."
    )]
    PairedCountMismatch {
        matched: usize,
        longer: &'static str,
    },

    /// Indicates that a mapped read was detected in the input alignment file.
    #[error("Error: Mapped read detected, please use `rasusa aln` for aligned data")]
    MappedReadDetected,

    /// Indicates that one-pass SAM/BAM/CRAM subsampling's name-grouped guard found evidence
    /// that the input is not grouped by read name, so streaming template-grouping (which only
    /// compares each record to the one immediately before it) cannot be trusted not to split
    /// mates apart.
    #[error(
        "Input does not appear to be grouped by read name: a paired/segmented read's records \
         must be adjacent for one-pass subsampling. Please collate the input first (e.g. \
         `samtools collate`), or mark it as name-grouped/name-sorted in the header (GO:query or \
         SO:queryname) if it already is."
    )]
    UngroupedAlignmentInput,
}

/// A `Struct` used for seamlessly dealing with either compressed or uncompressed fasta/fastq files.
#[derive(Debug, PartialEq)]
pub struct Fastx {
    /// The path for the file.
    path: PathBuf,
}

impl Fastx {
    /// Create a `Fastx` object from a `std::path::Path`.
    ///
    /// # Example
    ///
    /// ```rust
    /// use rasusa::fastx::Fastx;
    ///
    /// let path = std::path::Path::new("input.fa.gz");
    /// let fastx = Fastx::from_path(path);
    /// ```
    pub fn from_path(path: &Path) -> Self {
        Fastx {
            path: path.to_path_buf(),
        }
    }
}

impl Fastx {
    /// Opens the underlying (possibly compressed) file as a FASTA/Q record reader.
    ///
    /// Returns `Ok(None)` for a too-short/empty file, matching the pre-existing
    /// "no records" handling shared by [`RecordSource::read_lengths`] and
    /// [`RecordSource::count`].
    fn open_reader(&self) -> Result<Option<Box<dyn needletail::parser::FastxReader>>, FastxError> {
        let reader = match niffler::send::from_path(&self.path) {
            Ok((rdr, _)) => rdr,
            Err(niffler::error::Error::FileTooShort) => return Ok(None),
            Err(source) => return Err(FastxError::CompressOutputError(source)),
        };
        match needletail::parse_fastx_reader(reader) {
            Ok(rdr) => Ok(Some(rdr)),
            Err(e) if e.kind == EmptyFile => Ok(None),
            Err(source) => Err(FastxError::ReadError { source }),
        }
    }

    /// Reads the file once, writing each record to `write_to` independently with probability
    /// `fraction`, and returns how many records were seen and kept.
    ///
    /// Unlike [`RecordSource::filter_reads_into`] (which needs a pre-computed keep/discard
    /// decision per read, gathered by a prior pass), this makes each read's decision on the fly
    /// as it's parsed, so memory use stays constant in the size of the input: nothing scales with
    /// the number of reads. Output is written in input order, since each read is decided and
    /// emitted before the next one is parsed.
    pub fn subsample_one_pass(
        &self,
        fraction: f32,
        seed: Option<u64>,
        write_to: &mut dyn Write,
        fasta: bool,
    ) -> Result<OnePassStats, FastxError> {
        let mut rng = seeded_rng(seed);

        let mut reader = match self.open_reader()? {
            Some(rdr) => rdr,
            None => return Ok(OnePassStats::default()),
        };

        let mut stats = OnePassStats::default();

        while let Some(record) = reader.next() {
            let rec = record.map_err(|source| FastxError::ParseError { source })?;
            stats.reads_seen += 1;
            if rng.random_bool(fraction as f64) {
                crate::record::write_fastx_record(
                    write_to,
                    rec.id(),
                    &rec.seq(),
                    rec.qual(),
                    fasta,
                    &rec.line_ending().to_bytes(),
                )?;
                stats.reads_kept += 1;
            }
        }

        Ok(stats)
    }

    /// Reads two paired FASTA/Q files once, advancing both readers in lockstep: a read and its
    /// mate (a "template") are decided by a single shared random draw, so they are always kept
    /// or dropped together, never split.
    ///
    /// Neither file is buffered - both are read one record at a time - so peak memory stays
    /// constant in input size, the same as [`Fastx::subsample_one_pass`]. If the two files
    /// disagree on read count, that's detected the moment one reader runs out before the other,
    /// i.e. mid-stream, before either output is completely written - see
    /// [`FastxError::PairedCountMismatch`].
    #[allow(clippy::too_many_arguments)]
    pub fn subsample_one_pass_paired(
        &self,
        mate: &Fastx,
        fraction: f32,
        seed: Option<u64>,
        write_to: &mut dyn Write,
        mate_write_to: &mut dyn Write,
        fasta: bool,
        mate_fasta: bool,
    ) -> Result<OnePassStats, FastxError> {
        let mut rng = seeded_rng(seed);

        let mut reader = self.open_reader()?;
        let mut mate_reader = mate.open_reader()?;

        let mut stats = OnePassStats::default();

        loop {
            let record = reader.as_mut().and_then(|r| r.next());
            let mate_record = mate_reader.as_mut().and_then(|r| r.next());

            let (record, mate_record) = match (record, mate_record) {
                (None, None) => break,
                (Some(record), Some(mate_record)) => (record, mate_record),
                (Some(_), None) => {
                    return Err(FastxError::PairedCountMismatch {
                        matched: stats.reads_seen,
                        longer: "first",
                    })
                }
                (None, Some(_)) => {
                    return Err(FastxError::PairedCountMismatch {
                        matched: stats.reads_seen,
                        longer: "second",
                    })
                }
            };

            let rec = record.map_err(|source| FastxError::ParseError { source })?;
            let mate_rec = mate_record.map_err(|source| FastxError::ParseError { source })?;
            stats.reads_seen += 1;

            if rng.random_bool(fraction as f64) {
                crate::record::write_fastx_record(
                    write_to,
                    rec.id(),
                    &rec.seq(),
                    rec.qual(),
                    fasta,
                    &rec.line_ending().to_bytes(),
                )?;
                crate::record::write_fastx_record(
                    mate_write_to,
                    mate_rec.id(),
                    &mate_rec.seq(),
                    mate_rec.qual(),
                    mate_fasta,
                    &mate_rec.line_ending().to_bytes(),
                )?;
                stats.reads_kept += 1;
            }
        }

        Ok(stats)
    }
}

/// The outcome of a [`Fastx::subsample_one_pass`] run: how many reads were seen in the single
/// pass over the input, and how many of those were kept.
#[derive(Debug, Default, Clone, Copy, PartialEq, Eq)]
pub struct OnePassStats {
    pub reads_seen: usize,
    pub reads_kept: usize,
}

impl OnePassStats {
    /// The realised fraction of reads kept - `0.0` for an empty input, rather than `NaN`.
    pub fn realised_fraction(&self) -> f64 {
        if self.reads_seen == 0 {
            0.0
        } else {
            self.reads_kept as f64 / self.reads_seen as f64
        }
    }
}

impl RecordSource for Fastx {
    /// Returns a vector containing the lengths of all the reads in the file.
    ///
    /// # Errors
    /// If the file cannot be opened or there is an issue parsing any records then an
    /// `Err` containing a variant of [`FastxError`](#fastxerror) is returned.
    ///
    /// # Example
    ///
    /// ```rust
    /// use rasusa::fastx::Fastx;
    /// use rasusa::source::RecordSource;
    /// use std::io::Write;
    /// let text = "@read1\nACGT\n+\n!!!!\n@read2\nG\n+\n!";
    /// let mut file = tempfile::Builder::new().suffix(".fq").tempfile().unwrap();
    /// file.write_all(text.as_bytes()).unwrap();
    /// let fastx = Fastx::from_path(file.path());
    /// let actual = fastx.read_lengths().unwrap();
    /// let expected: Vec<u32> = vec![4, 1];
    /// assert_eq!(actual, expected)
    /// ```
    fn read_lengths(&self) -> Result<Vec<u32>, FastxError> {
        let mut read_lengths: Vec<u32> = vec![];

        let mut reader = match self.open_reader()? {
            Some(rdr) => rdr,
            None => return Ok(read_lengths),
        };

        while let Some(record) = reader.next() {
            match record {
                Ok(rec) => read_lengths.push(rec.num_bases() as u32),
                Err(err) => return Err(FastxError::ParseError { source: err }),
            }
        }
        Ok(read_lengths)
    }

    /// Returns the number of records in the file, without materializing their lengths.
    ///
    /// # Example
    ///
    /// ```rust
    /// use rasusa::fastx::Fastx;
    /// use rasusa::source::RecordSource;
    /// use std::io::Write;
    /// let text = "@read1\nACGT\n+\n!!!!\n@read2\nG\n+\n!";
    /// let mut file = tempfile::Builder::new().suffix(".fq").tempfile().unwrap();
    /// file.write_all(text.as_bytes()).unwrap();
    /// let fastx = Fastx::from_path(file.path());
    /// assert_eq!(fastx.count().unwrap(), 2)
    /// ```
    fn count(&self) -> Result<usize, FastxError> {
        let mut count: usize = 0;

        let mut reader = match self.open_reader()? {
            Some(rdr) => rdr,
            None => return Ok(count),
        };

        while let Some(record) = reader.next() {
            match record {
                Ok(_) => count += 1,
                Err(err) => return Err(FastxError::ParseError { source: err }),
            }
        }
        Ok(count)
    }

    /// Writes reads, with indices contained within `reads_to_keep`, to the specified handle
    /// `write_to`.
    ///
    /// # Errors
    /// This function could raise an `Err` instance of [`FastxError`](#fastxerror) in the following
    /// circumstances:
    /// -   If the file (of `self`) cannot be opened.
    /// -   If writing to `write_to` fails.
    /// -   If, after iterating through all reads in the file, there is still elements left in
    ///     `reads_to_keep`. *Note: in this case, this function still writes all reads where indices
    ///     were found in the file.*
    fn filter_reads_into(
        &self,
        reads_to_keep: &[bool],
        nb_reads_keep: usize,
        write_to: &mut dyn Write,
        encoding: OutputEncoding,
    ) -> Result<usize, FastxError> {
        // Callers never ask a FASTA/Q source for alignment-format output (`Reads::run` rejects
        // that combination before any `RecordSource` is touched), so this is the only branch
        // this impl needs to handle.
        let is_fasta = match encoding {
            OutputEncoding::Fastx { fasta } => fasta,
            OutputEncoding::Alignment(_) => {
                unreachable!("fastx sources never receive alignment output encoding")
            }
        };

        let mut total_len = 0;

        let (reader, _) = niffler::send::from_path(&self.path)?;
        let mut reader = needletail::parse_fastx_reader(reader)
            .map_err(|source| FastxError::ReadError { source })?;
        let mut read_idx: usize = 0;
        let mut nb_reads_written = 0;

        while let Some(record) = reader.next() {
            match record {
                Err(source) => return Err(FastxError::ParseError { source }),
                Ok(rec) if read_idx < reads_to_keep.len() && reads_to_keep[read_idx] => {
                    total_len += rec.num_bases();
                    crate::record::write_fastx_record(
                        write_to,
                        rec.id(),
                        &rec.seq(),
                        rec.qual(),
                        is_fasta,
                        &rec.line_ending().to_bytes(),
                    )?;
                    nb_reads_written += 1;
                    if nb_reads_keep == nb_reads_written {
                        break;
                    }
                }
                Ok(_) => (),
            }

            read_idx += 1;
        }

        if nb_reads_written == nb_reads_keep {
            Ok(total_len)
        } else {
            Err(FastxError::IndicesNotFound)
        }
    }
}

/// Create a file for writing.
///
/// # Errors
/// If the file cannot be created then an `Err` containing a variant of [`FastxError`](#fastxerror) is
/// returned.
pub fn create_output_writer(
    path: &Path,
    compression_lvl: Option<niffler::compression::Level>,
    compression_fmt: Option<niffler::compression::Format>,
) -> Result<Box<dyn Write>, FastxError> {
    let file = File::create(path).map_err(|source| FastxError::CreateError { source })?;
    let file_handle = Box::new(BufWriter::new(file));
    let fmt = compression_fmt.unwrap_or_else(|| crate::format::infer_compression_format(path));
    let compression_lvl =
        compression_lvl.unwrap_or_else(|| crate::format::default_compression_level(fmt));
    niffler::get_writer(file_handle, fmt, compression_lvl).map_err(FastxError::CompressOutputError)
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::any::Any;
    use std::io::{Read, Write};
    use std::path::Path;
    use tempfile::{Builder, NamedTempFile};

    /// Writes `text` to a temp `.fastq` file and wraps it in a `Fastx`. The returned
    /// `NamedTempFile` must be kept alive (bound to a variable) for as long as the `Fastx` is
    /// used, since dropping it deletes the underlying file.
    fn temp_fastx(text: &str) -> (NamedTempFile, Fastx) {
        let mut input = Builder::new().suffix(".fastq").tempfile().unwrap();
        input.write_all(text.as_bytes()).unwrap();
        let fastx = Fastx::from_path(input.path());
        (input, fastx)
    }

    #[test]
    fn fastx_from_fasta() {
        let path = Path::new("data/my.fa");

        let actual = Fastx::from_path(path);
        let expected = Fastx {
            path: path.to_path_buf(),
        };

        assert_eq!(actual, expected)
    }

    #[test]
    fn create_invalid_output_file_raises_error() {
        let path = Path::new("invalid/out/path.fq");

        let actual = create_output_writer(path, Some(niffler::Level::Eight), None)
            .err()
            .unwrap();
        let expected = FastxError::CreateError {
            source: std::io::Error::other(String::from("No such file or directory (os error 2)")),
        };

        assert_eq!(actual.type_id(), expected.type_id())
    }

    #[test]
    fn create_valid_output_file_and_can_write_to_it() {
        let file = Builder::new().suffix(".fastq").tempfile().unwrap();
        let mut writer =
            create_output_writer(file.path(), Some(niffler::Level::Eight), None).unwrap();

        let actual = writer.write(b"foo\nbar");

        assert!(actual.is_ok())
    }

    #[test]
    fn create_valid_compressed_output_file_and_can_write_to_it() {
        let file = Builder::new().suffix(".fastq.gz").tempfile().unwrap();
        let mut writer =
            create_output_writer(file.path(), Some(niffler::Level::Four), None).unwrap();

        let actual = writer.write(b"foo\nbar");

        assert!(actual.is_ok())
    }

    #[test]
    fn get_read_lengths_for_empty_fasta_returns_empty_vector() {
        let text = "";
        let mut file = Builder::new().suffix(".fa").tempfile().unwrap();
        file.write_all(text.as_bytes()).unwrap();
        let fastx = Fastx::from_path(file.path());

        let actual = fastx.read_lengths().unwrap();
        let expected: Vec<u32> = Vec::new();

        assert_eq!(actual, expected)
    }

    #[test]
    fn get_read_lengths_for_fasta() {
        let text = ">read1\nACGT\n>read2\nG";
        let mut file = Builder::new().suffix(".fa").tempfile().unwrap();
        file.write_all(text.as_bytes()).unwrap();
        let fastx = Fastx::from_path(file.path());

        let actual = fastx.read_lengths().unwrap();
        let expected: Vec<u32> = vec![4, 1];

        assert_eq!(actual, expected)
    }

    #[test]
    fn get_read_lengths_for_fastq() {
        let text = "@read1\nACGT\n+\n!!!!\n@read2\nG\n+\n!";
        let mut file = Builder::new().suffix(".fq").tempfile().unwrap();
        file.write_all(text.as_bytes()).unwrap();
        let fastx = Fastx::from_path(file.path());

        let actual = fastx.read_lengths().unwrap();
        let expected: Vec<u32> = vec![4, 1];

        assert_eq!(actual, expected)
    }

    #[test]
    fn filter_reads_empty_indices_no_output() {
        let text = "@read1\nACGT\n+\n!!!!";
        let mut input = Builder::new().suffix(".fastq").tempfile().unwrap();
        input.write_all(text.as_bytes()).unwrap();
        let fastx = Fastx::from_path(input.path());
        let reads_to_keep: Vec<bool> = vec![false];
        let output = Builder::new().suffix(".fastq").tempfile().unwrap();
        let mut out_fh = create_output_writer(output.path(), None, None).unwrap();
        let filter_result = fastx.filter_reads_into(
            &reads_to_keep,
            0,
            &mut out_fh,
            OutputEncoding::Fastx { fasta: false },
        );

        assert!(filter_result.is_ok());

        let mut actual = String::new();
        output.into_file().read_to_string(&mut actual).unwrap();
        let expected = String::new();

        assert_eq!(actual, expected)
    }

    #[test]
    fn filter_fastq_reads_one_index_matches_only_read() {
        let text = "@read1\nACGT\n+\n!!!!\n";
        let mut input = Builder::new().suffix(".fastq").tempfile().unwrap();
        input.write_all(text.as_bytes()).unwrap();
        let fastx = Fastx::from_path(input.path());
        let reads_to_keep: Vec<bool> = vec![true];
        let output = Builder::new().suffix(".fastq").tempfile().unwrap();
        {
            let mut out_fh = create_output_writer(output.path(), None, None).unwrap();
            let filter_result = fastx.filter_reads_into(
                &reads_to_keep,
                1,
                &mut out_fh,
                OutputEncoding::Fastx { fasta: false },
            );
            assert!(filter_result.is_ok());
        }

        let actual = std::fs::read_to_string(output).unwrap();
        let expected = text;

        assert_eq!(actual, expected)
    }

    #[test]
    fn filter_fasta_reads_one_index_matches_only_read() {
        let text = ">read1\nACGT\n";
        let mut input = Builder::new().suffix(".fa").tempfile().unwrap();
        input.write_all(text.as_bytes()).unwrap();
        let fastx = Fastx::from_path(input.path());
        let reads_to_keep: Vec<bool> = vec![true];
        let output = Builder::new().suffix(".fa").tempfile().unwrap();
        {
            let mut out_fh = create_output_writer(output.path(), None, None).unwrap();
            let filter_result = fastx.filter_reads_into(
                &reads_to_keep,
                1,
                &mut out_fh,
                OutputEncoding::Fastx { fasta: true },
            );
            assert!(filter_result.is_ok());
        }

        let actual = std::fs::read_to_string(output).unwrap();
        let expected = text;

        assert_eq!(actual, expected)
    }

    #[test]
    fn filter_fastq_reads_one_index_matches_one_of_two_reads() {
        let text = "@read1\nACGT\n+\n!!!!\n@read2\nCCCC\n+\n$$$$\n";
        let mut input = Builder::new().suffix(".fastq").tempfile().unwrap();
        input.write_all(text.as_bytes()).unwrap();
        let fastx = Fastx::from_path(input.path());
        let reads_to_keep: Vec<bool> = vec![false, true];
        let output = Builder::new().suffix(".fastq").tempfile().unwrap();
        {
            let mut out_fh = create_output_writer(output.path(), None, None).unwrap();
            let filter_result = fastx.filter_reads_into(
                &reads_to_keep,
                1,
                &mut out_fh,
                OutputEncoding::Fastx { fasta: false },
            );
            assert!(filter_result.is_ok());
        }

        let actual = std::fs::read_to_string(output).unwrap();
        let expected = "@read2\nCCCC\n+\n$$$$\n";

        assert_eq!(actual, expected)
    }

    #[test]
    fn filter_fastq_reads_two_indices_matches_first_and_last_reads() {
        let text = "@read1\nACGT\n+\n!!!!\n@read2\nCCCC\n+\n$$$$\n@read3\nA\n+\n$\n";
        let mut input = Builder::new().suffix(".fastq").tempfile().unwrap();
        input.write_all(text.as_bytes()).unwrap();
        let fastx = Fastx::from_path(input.path());
        let reads_to_keep: Vec<bool> = vec![true, false, true];
        let output = Builder::new().suffix(".fastq").tempfile().unwrap();
        {
            let mut out_fh = create_output_writer(output.path(), None, None).unwrap();
            let filter_result = fastx.filter_reads_into(
                &reads_to_keep,
                2,
                &mut out_fh,
                OutputEncoding::Fastx { fasta: false },
            );
            assert!(filter_result.is_ok());
        }

        let actual = std::fs::read_to_string(output).unwrap();
        let expected = "@read1\nACGT\n+\n!!!!\n@read3\nA\n+\n$\n";

        assert_eq!(actual, expected)
    }

    #[test]
    fn filter_fasta_reads_one_index_out_of_range() {
        let text = ">read1 length=4\nACGT\n>read2\nCCCC\n";
        let mut input = Builder::new().suffix(".fa").tempfile().unwrap();
        input.write_all(text.as_bytes()).unwrap();
        let fastx = Fastx::from_path(input.path());
        let reads_to_keep: Vec<bool> = vec![true, false, true];
        let output = Builder::new().suffix(".fa").tempfile().unwrap();
        {
            let mut out_fh =
                create_output_writer(output.path(), Some(niffler::Level::Four), None).unwrap();
            let filter_result = fastx.filter_reads_into(
                &reads_to_keep,
                2,
                &mut out_fh,
                OutputEncoding::Fastx { fasta: true },
            );
            assert!(filter_result.is_err());
        }

        let actual = std::fs::read_to_string(output).unwrap();
        let expected = ">read1 length=4\nACGT\n";

        assert_eq!(actual, expected)
    }

    #[test]
    fn filter_fastq_reads_one_index_out_of_range() {
        let text = "@read1 length=4\nACGT\n+\n!!!!\n@read2\nC\n+\n^\n";
        let mut input = Builder::new().suffix(".fq").tempfile().unwrap();
        input.write_all(text.as_bytes()).unwrap();
        let fastx = Fastx::from_path(input.path());
        let reads_to_keep: Vec<bool> = vec![true, false, true];
        let output = Builder::new().suffix(".fq").tempfile().unwrap();
        {
            let mut out_fh =
                create_output_writer(output.path(), Some(niffler::Level::Four), None).unwrap();
            let filter_result = fastx.filter_reads_into(
                &reads_to_keep,
                2,
                &mut out_fh,
                OutputEncoding::Fastx { fasta: false },
            );
            assert!(filter_result.is_err());
        }

        let actual = std::fs::read_to_string(output).unwrap();
        let expected = "@read1 length=4\nACGT\n+\n!!!!\n";

        assert_eq!(actual, expected)
    }

    #[test]
    fn one_pass_fraction_one_keeps_every_read() {
        let text = "@read1\nACGT\n+\n!!!!\n@read2\nCCCC\n+\n$$$$\n";
        let (_input, fastx) = temp_fastx(text);

        let mut out: Vec<u8> = Vec::new();
        let stats = fastx
            .subsample_one_pass(1.0, Some(1), &mut out, false)
            .unwrap();

        assert_eq!(stats.reads_seen, 2);
        assert_eq!(stats.reads_kept, 2);
        assert_eq!(String::from_utf8(out).unwrap(), text);
    }

    #[test]
    fn one_pass_fraction_zero_keeps_no_reads() {
        let text = "@read1\nACGT\n+\n!!!!\n@read2\nCCCC\n+\n$$$$\n";
        let (_input, fastx) = temp_fastx(text);

        let mut out: Vec<u8> = Vec::new();
        let stats = fastx
            .subsample_one_pass(0.0, Some(1), &mut out, false)
            .unwrap();

        assert_eq!(stats.reads_seen, 2);
        assert_eq!(stats.reads_kept, 0);
        assert!(out.is_empty());
    }

    #[test]
    fn one_pass_empty_file_returns_zeroed_stats() {
        let (_input, fastx) = temp_fastx("");

        let mut out: Vec<u8> = Vec::new();
        let stats = fastx
            .subsample_one_pass(0.5, Some(1), &mut out, false)
            .unwrap();

        assert_eq!(stats, OnePassStats::default());
        assert!(out.is_empty());
    }

    #[test]
    fn one_pass_same_seed_gives_same_result() {
        let text = "@read1\nACGT\n+\n!!!!\n@read2\nCCCC\n+\n$$$$\n@read3\nGGGG\n+\n####\n@read4\nTTTT\n+\n^^^^\n";
        let (_input, fastx) = temp_fastx(text);

        let mut out1: Vec<u8> = Vec::new();
        let stats1 = fastx
            .subsample_one_pass(0.5, Some(42), &mut out1, false)
            .unwrap();

        let mut out2: Vec<u8> = Vec::new();
        let stats2 = fastx
            .subsample_one_pass(0.5, Some(42), &mut out2, false)
            .unwrap();

        assert_eq!(stats1, stats2);
        assert_eq!(out1, out2);
    }

    #[test]
    fn one_pass_preserves_input_order() {
        let text = "@read1\nACGT\n+\n!!!!\n@read2\nCCCC\n+\n$$$$\n@read3\nGGGG\n+\n####\n@read4\nTTTT\n+\n^^^^\n@read5\nAAAA\n+\n%%%%\n";
        let (_input, fastx) = temp_fastx(text);

        let mut out: Vec<u8> = Vec::new();
        fastx
            .subsample_one_pass(0.6, Some(7), &mut out, false)
            .unwrap();

        let output = String::from_utf8(out).unwrap();
        let kept_ids: Vec<&str> = output.lines().filter(|l| l.starts_with('@')).collect();
        let mut sorted_ids = kept_ids.clone();
        sorted_ids.sort();

        assert_eq!(kept_ids, sorted_ids);
        assert!(!kept_ids.is_empty());
    }

    #[test]
    fn one_pass_writes_fasta_when_requested() {
        let (_input, fastx) = temp_fastx("@read1\nACGT\n+\n!!!!\n");

        let mut out: Vec<u8> = Vec::new();
        let stats = fastx
            .subsample_one_pass(1.0, Some(1), &mut out, true)
            .unwrap();

        assert_eq!(stats.reads_kept, 1);
        assert_eq!(String::from_utf8(out).unwrap(), ">read1\nACGT\n");
    }

    #[test]
    fn one_pass_stats_realised_fraction() {
        let stats = OnePassStats {
            reads_seen: 4,
            reads_kept: 1,
        };
        assert_eq!(stats.realised_fraction(), 0.25);
        assert_eq!(OnePassStats::default().realised_fraction(), 0.0);
    }

    #[test]
    fn one_pass_paired_fraction_one_keeps_every_template() {
        let r1 = "@read1\nACGT\n+\n!!!!\n@read2\nCCCC\n+\n$$$$\n";
        let r2 = "@read1\nTTTT\n+\n!!!!\n@read2\nGGGG\n+\n$$$$\n";
        let (_in1, fastx1) = temp_fastx(r1);
        let (_in2, fastx2) = temp_fastx(r2);

        let mut out1: Vec<u8> = Vec::new();
        let mut out2: Vec<u8> = Vec::new();
        let stats = fastx1
            .subsample_one_pass_paired(&fastx2, 1.0, Some(1), &mut out1, &mut out2, false, false)
            .unwrap();

        assert_eq!(stats.reads_seen, 2);
        assert_eq!(stats.reads_kept, 2);
        assert_eq!(String::from_utf8(out1).unwrap(), r1);
        assert_eq!(String::from_utf8(out2).unwrap(), r2);
    }

    #[test]
    fn one_pass_paired_fraction_zero_keeps_no_templates() {
        let r1 = "@read1\nACGT\n+\n!!!!\n@read2\nCCCC\n+\n$$$$\n";
        let r2 = "@read1\nTTTT\n+\n!!!!\n@read2\nGGGG\n+\n$$$$\n";
        let (_in1, fastx1) = temp_fastx(r1);
        let (_in2, fastx2) = temp_fastx(r2);

        let mut out1: Vec<u8> = Vec::new();
        let mut out2: Vec<u8> = Vec::new();
        let stats = fastx1
            .subsample_one_pass_paired(&fastx2, 0.0, Some(1), &mut out1, &mut out2, false, false)
            .unwrap();

        assert_eq!(stats.reads_seen, 2);
        assert_eq!(stats.reads_kept, 0);
        assert!(out1.is_empty());
        assert!(out2.is_empty());
    }

    #[test]
    fn one_pass_paired_keeps_mates_together() {
        let r1 = "@read1\nAAAA\n+\n!!!!\n@read2\nCCCC\n+\n$$$$\n@read3\nGGGG\n+\n####\n@read4\nTTTT\n+\n^^^^\n@read5\nAAAA\n+\n%%%%\n";
        let r2 = "@read1\nAAAA\n+\n!!!!\n@read2\nCCCC\n+\n$$$$\n@read3\nGGGG\n+\n####\n@read4\nTTTT\n+\n^^^^\n@read5\nAAAA\n+\n%%%%\n";
        let (_in1, fastx1) = temp_fastx(r1);
        let (_in2, fastx2) = temp_fastx(r2);

        let mut out1: Vec<u8> = Vec::new();
        let mut out2: Vec<u8> = Vec::new();
        let stats = fastx1
            .subsample_one_pass_paired(&fastx2, 0.6, Some(7), &mut out1, &mut out2, false, false)
            .unwrap();

        let ids1: Vec<String> = String::from_utf8(out1)
            .unwrap()
            .lines()
            .filter(|l| l.starts_with('@'))
            .map(str::to_owned)
            .collect();
        let ids2: Vec<String> = String::from_utf8(out2)
            .unwrap()
            .lines()
            .filter(|l| l.starts_with('@'))
            .map(str::to_owned)
            .collect();

        assert_eq!(ids1, ids2);
        assert_eq!(stats.reads_kept, ids1.len());
        assert!(!ids1.is_empty());
    }

    #[test]
    fn one_pass_paired_same_seed_gives_same_pair_selection() {
        let text = "@read1\nAAAA\n+\n!!!!\n@read2\nCCCC\n+\n$$$$\n@read3\nGGGG\n+\n####\n@read4\nTTTT\n+\n^^^^\n";
        let (_in1, fastx1) = temp_fastx(text);
        let (_in2, fastx2) = temp_fastx(text);
        let (_in1b, fastx1b) = temp_fastx(text);
        let (_in2b, fastx2b) = temp_fastx(text);

        let mut out1a: Vec<u8> = Vec::new();
        let mut out2a: Vec<u8> = Vec::new();
        let stats_a = fastx1
            .subsample_one_pass_paired(&fastx2, 0.5, Some(42), &mut out1a, &mut out2a, false, false)
            .unwrap();

        let mut out1b: Vec<u8> = Vec::new();
        let mut out2b: Vec<u8> = Vec::new();
        let stats_b = fastx1b
            .subsample_one_pass_paired(
                &fastx2b,
                0.5,
                Some(42),
                &mut out1b,
                &mut out2b,
                false,
                false,
            )
            .unwrap();

        assert_eq!(stats_a, stats_b);
        assert_eq!(out1a, out1b);
        assert_eq!(out2a, out2b);
    }

    #[test]
    fn one_pass_paired_read_count_mismatch_is_detected_mid_stream() {
        let r1 = "@read1\nAAAA\n+\n!!!!\n@read2\nCCCC\n+\n$$$$\n@read3\nGGGG\n+\n####\n";
        let r2 = "@read1\nAAAA\n+\n!!!!\n";
        let (_in1, fastx1) = temp_fastx(r1);
        let (_in2, fastx2) = temp_fastx(r2);

        let mut out1: Vec<u8> = Vec::new();
        let mut out2: Vec<u8> = Vec::new();
        // fraction 0.0 so nothing gets written, isolating the mismatch detection itself
        let result = fastx1.subsample_one_pass_paired(
            &fastx2,
            0.0,
            Some(1),
            &mut out1,
            &mut out2,
            false,
            false,
        );

        let err = result.expect_err("mismatched paired read counts should be an error");
        assert!(matches!(
            err,
            FastxError::PairedCountMismatch {
                matched: 1,
                longer: "first"
            }
        ));
    }

    #[test]
    fn one_pass_paired_both_empty_returns_zeroed_stats() {
        let (_in1, fastx1) = temp_fastx("");
        let (_in2, fastx2) = temp_fastx("");

        let mut out1: Vec<u8> = Vec::new();
        let mut out2: Vec<u8> = Vec::new();
        let stats = fastx1
            .subsample_one_pass_paired(&fastx2, 0.5, Some(1), &mut out1, &mut out2, false, false)
            .unwrap();

        assert_eq!(stats, OnePassStats::default());
        assert!(out1.is_empty());
        assert!(out2.is_empty());
    }
}