genomicframe_core/formats/fastq/

reader.rs

//! FASTQ format support (sequence data with quality scores)
//!
//! FASTQ is a text-based format for representing nucleotide sequences
//! with per-base quality scores. Each record consists of 4 lines:
//! 1. Header: @ID description
//! 2. Sequence: DNA/RNA bases
//! 3. Separator: + (optionally followed by ID)
//! 4. Quality: ASCII-encoded Phred scores
//!
//! # Design Philosophy
//!
//! - **Streaming by default**: Records are never buffered unless explicitly requested
//! - **Minimal allocations**: Reuses buffers where possible
//! - **Automatic compression detection**: Handles `.fastq`, `.fq`, `.fastq.gz`, `.fq.gz`
//! - **Lazy evaluation**: Process millions of reads with O(1) memory
//! - **Quality score parsing**: Phred+33 (Sanger) and Phred+64 (Illumina 1.3+) support
//!
//! # Examples
//!
//! ```no_run
//! use genomicframe_core::formats::fastq::FastqReader;
//! use genomicframe_core::core::GenomicRecordIterator;
//!
//! // Streaming: O(1) memory, processes one read at a time
//! let mut reader = FastqReader::from_path("reads.fastq.gz")?;
//!
//! while let Some(record) = reader.next_record()? {
//!     let mean_quality = record.mean_quality();
//!     if mean_quality >= 30.0 {
//!         println!("High quality read: {}", record.id);
//!     }
//! }
//! # Ok::<(), genomicframe_core::error::Error>(())
//! ```

use crate::core::{GenomicReader, GenomicRecordIterator};
use crate::error::{Error, Result};
use crate::io::Compression;
use flate2::read::MultiGzDecoder;
use std::fs::File;
use std::io::{BufRead, BufReader};
use std::path::Path;

/// Quality score encoding scheme
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum QualityEncoding {
    /// Phred+33 (Sanger, Illumina 1.8+) - standard modern format
    Phred33,
    /// Phred+64 (Illumina 1.3-1.7) - legacy format
    Phred64,
}

impl QualityEncoding {
    /// Get the ASCII offset for this encoding
    pub fn offset(&self) -> u8 {
        match self {
            QualityEncoding::Phred33 => 33,
            QualityEncoding::Phred64 => 64,
        }
    }

    /// Detect encoding from quality string (heuristic)
    pub fn detect(quality: &str) -> Self {
        // Check for characters that only appear in Phred+64
        // Phred+33: ! to ~  (33-126)
        // Phred+64: @ to ~  (64-126)
        // If we see characters < 64, it must be Phred+33
        for byte in quality.bytes() {
            if byte < 64 {
                return QualityEncoding::Phred33;
            }
        }
        // Default to Phred+33 (modern standard)
        QualityEncoding::Phred33
    }
}

/// A single FASTQ record (DNA/RNA sequence with quality scores)
#[derive(Debug, Clone)]
pub struct FastqRecord {
    /// Sequence identifier (header line without '@')
    pub id: String,
    /// Optional description (text after first whitespace in header)
    pub description: Option<String>,
    /// Sequence data (DNA/RNA bases)
    pub sequence: String,
    /// Quality scores (ASCII-encoded Phred scores)
    pub quality: String,
}

impl FastqRecord {
    /// Get sequence length (number of bases)
    pub fn len(&self) -> usize {
        self.sequence.len()
    }

    /// Check if sequence is empty
    pub fn is_empty(&self) -> bool {
        self.sequence.is_empty()
    }

    /// Convert quality string to Phred scores
    pub fn phred_scores(&self, encoding: QualityEncoding) -> Vec<u8> {
        let offset = encoding.offset();
        self.quality
            .bytes()
            .map(|b| b.saturating_sub(offset))
            .collect()
    }

    /// Get mean quality score
    pub fn mean_quality(&self) -> f64 {
        self.mean_quality_with_encoding(QualityEncoding::Phred33)
    }

    /// Get mean quality score with specific encoding
    pub fn mean_quality_with_encoding(&self, encoding: QualityEncoding) -> f64 {
        let scores = self.phred_scores(encoding);
        if scores.is_empty() {
            return 0.0;
        }
        let sum: u32 = scores.iter().map(|&s| s as u32).sum();
        sum as f64 / scores.len() as f64
    }

    /// Check if sequence and quality have matching lengths
    pub fn is_valid(&self) -> bool {
        self.sequence.len() == self.quality.len()
    }

    /// Get GC content as a fraction (0.0 to 1.0)
    pub fn gc_content(&self) -> f64 {
        let gc_count = self
            .sequence
            .bytes()
            .filter(|&b| b == b'G' || b == b'C' || b == b'g' || b == b'c')
            .count();
        if self.sequence.is_empty() {
            0.0
        } else {
            gc_count as f64 / self.sequence.len() as f64
        }
    }

    /// Get N content (ambiguous bases) as a fraction (0.0 to 1.0)
    pub fn n_content(&self) -> f64 {
        let n_count = self
            .sequence
            .bytes()
            .filter(|&b| b == b'N' || b == b'n')
            .count();
        if self.sequence.is_empty() {
            0.0
        } else {
            n_count as f64 / self.sequence.len() as f64
        }
    }
}

/// FASTQ header metadata (minimal - FASTQ has no standard header)
#[derive(Debug, Clone, Default)]
pub struct FastqHeader {
    /// Detected quality encoding (set after reading first record)
    pub encoding: Option<QualityEncoding>,
}

/// FASTQ file reader (streaming, memory-efficient)
///
/// Create using `FastqReader::from_path()` or `FastqReader::from_buffer()`.
///
/// FASTQ format:
/// - Line 1: @ID description
/// - Line 2: Sequence
/// - Line 3: + (optional ID)
/// - Line 4: Quality
pub struct FastqReader {
    reader: Box<dyn BufRead>,
    header: FastqHeader,
    line_number: usize,
    // Reusable buffers to avoid allocations
    id_buffer: String,
    seq_buffer: String,
    sep_buffer: String,
    qual_buffer: String,
}

impl FastqReader {
    /// Create a FASTQ reader from a buffered reader
    ///
    /// This is the low-level constructor that accepts any buffered reader.
    /// Use this when you already have a BufRead source (e.g., from seeking
    /// to a specific position in a file for parallel processing).
    ///
    /// # Example
    ///
    /// ```no_run
    /// use genomicframe_core::formats::fastq::FastqReader;
    /// use std::io::BufReader;
    /// use std::fs::File;
    ///
    /// let file = File::open("reads.fastq")?;
    /// let buf_reader = BufReader::new(file);
    /// let mut reader = FastqReader::from_buffer(buf_reader);
    /// # Ok::<(), genomicframe_core::error::Error>(())
    /// ```
    pub fn from_buffer<R: BufRead + 'static>(reader: R) -> Self {
        Self {
            reader: Box::new(reader),
            header: FastqHeader::default(),
            line_number: 0,
            // Pre-allocate buffers with reasonable capacity (typical read ~150bp)
            id_buffer: String::with_capacity(256),
            seq_buffer: String::with_capacity(256),
            sep_buffer: String::with_capacity(4),
            qual_buffer: String::with_capacity(256),
        }
    }

    /// Open a FASTQ file - handles .fastq, .fq, .fastq.gz, .fq.gz automatically
    ///
    /// # Memory Behavior
    ///
    /// - Opens file handle (~8KB buffer)
    /// - **Does NOT load reads into memory**
    /// - Records are parsed on-demand as you iterate
    ///
    /// # Example
    ///
    /// ```no_run
    /// use genomicframe_core::formats::fastq::FastqReader;
    /// use genomicframe_core::core::GenomicRecordIterator;
    ///
    /// // Works with multi-GB files using only ~10KB RAM
    /// let mut reader = FastqReader::from_path("reads.fastq.gz")?;
    ///
    /// while let Some(record) = reader.next_record()? {
    ///     if record.mean_quality() >= 30.0 {
    ///         println!("High quality: {}", record.id);
    ///     }
    /// }
    /// # Ok::<(), genomicframe_core::error::Error>(())
    /// ```
    pub fn from_path<P: AsRef<Path>>(path: P) -> Result<Self> {
        let path = path.as_ref();
        let file = File::open(path)?;
        let compression = Compression::from_path(path);

        let reader: Box<dyn BufRead> = match compression {
            Compression::Gzip | Compression::Bgzip => {
                Box::new(BufReader::new(MultiGzDecoder::new(file)))
            }
            _ => Box::new(BufReader::new(file)),
        };

        Ok(Self {
            reader,
            header: FastqHeader::default(),
            line_number: 0,
            // Pre-allocate buffers with reasonable capacity (typical read ~150bp)
            id_buffer: String::with_capacity(256),
            seq_buffer: String::with_capacity(256),
            sep_buffer: String::with_capacity(4),
            qual_buffer: String::with_capacity(256),
        })
    }

    /// Get the FASTQ header metadata
    pub fn header(&self) -> &FastqHeader {
        &self.header
    }


    /// Parse a FASTQ record (4 lines)
    fn parse_record(&mut self) -> Result<Option<FastqRecord>> {
        // Line 1: Header (@ID description)
        self.id_buffer.clear();
        let bytes = self.reader.read_line(&mut self.id_buffer)?;
        if bytes > 0 {
            self.line_number += 1;
        } else {
            return Ok(None); // End of file
        }

        // Trim in-place by finding actual content length (strip trailing whitespace)
        let id_line = self.id_buffer.trim();
        if !id_line.starts_with('@') {
            return Err(Error::Parse(format!(
                "Line {}: Expected '@' at start of FASTQ record, got: {}",
                self.line_number,
                id_line.chars().take(50).collect::<String>()
            )));
        }

        // Parse ID and description
        let header_content = &id_line[1..]; // Skip '@'
        let (id, description) = if let Some((id, desc)) = header_content.split_once(char::is_whitespace) {
            (id.to_string(), Some(desc.trim().to_string()))
        } else {
            (header_content.to_string(), None)
        };

        // Line 2: Sequence
        self.seq_buffer.clear();
        let bytes = self.reader.read_line(&mut self.seq_buffer)?;
        if bytes > 0 {
            self.line_number += 1;
        } else {
            return Err(Error::Parse(format!(
                "Line {}: Unexpected end of file while reading sequence",
                self.line_number
            )));
        }
        let sequence = self.seq_buffer.trim().to_string();

        // Line 3: Separator (+ optional ID)
        self.sep_buffer.clear();
        let bytes = self.reader.read_line(&mut self.sep_buffer)?;
        if bytes > 0 {
            self.line_number += 1;
        } else {
            return Err(Error::Parse(format!(
                "Line {}: Unexpected end of file while reading separator",
                self.line_number
            )));
        }
        let separator = self.sep_buffer.trim();
        if !separator.starts_with('+') {
            return Err(Error::Parse(format!(
                "Line {}: Expected '+' separator, got: {}",
                self.line_number, separator
            )));
        }

        // Line 4: Quality
        self.qual_buffer.clear();
        let bytes = self.reader.read_line(&mut self.qual_buffer)?;
        if bytes > 0 {
            self.line_number += 1;
        } else {
            return Err(Error::Parse(format!(
                "Line {}: Unexpected end of file while reading quality",
                self.line_number
            )));
        }
        let quality = self.qual_buffer.trim().to_string();

        // Validate sequence and quality lengths match
        if sequence.len() != quality.len() {
            return Err(Error::Parse(format!(
                "Line {}: Sequence length ({}) does not match quality length ({})",
                self.line_number,
                sequence.len(),
                quality.len()
            )));
        }

        // Auto-detect quality encoding from first record
        if self.header.encoding.is_none() && !quality.is_empty() {
            self.header.encoding = Some(QualityEncoding::detect(&quality));
        }

        Ok(Some(FastqRecord {
            id,
            description,
            sequence,
            quality,
        }))
    }
}

impl GenomicRecordIterator for FastqReader {
    type Record = FastqRecord;

    fn next_raw(&mut self) -> Result<Option<Vec<u8>>> { 
        // TODO: Implement
        Ok(None)
    }

    fn next_record(&mut self) -> Result<Option<Self::Record>> {
        self.parse_record()
    }
}

impl GenomicReader for FastqReader {
    type Metadata = FastqHeader;

    fn metadata(&self) -> &Self::Metadata {
        &self.header
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::io::Write;
    use tempfile::NamedTempFile;

    #[test]
    fn test_fastq_basic_read() -> Result<()> {
        let fastq_data = "@READ1 description here\n\
                          ACGTACGTACGT\n\
                          +\n\
                          IIIIIIIIIIII\n\
                          @READ2\n\
                          GGGGCCCCAAAA\n\
                          +READ2\n\
                          HHHHHHHHHHHH\n";

        let mut temp_file = NamedTempFile::new()?;
        temp_file.write_all(fastq_data.as_bytes())?;
        temp_file.flush()?;

        let mut reader = FastqReader::from_path(temp_file.path())?;

        // First record
        let record1 = reader.next_record()?.expect("Should have first record");
        assert_eq!(record1.id, "READ1");
        assert_eq!(record1.description, Some("description here".to_string()));
        assert_eq!(record1.sequence, "ACGTACGTACGT");
        assert_eq!(record1.quality, "IIIIIIIIIIII");
        assert_eq!(record1.len(), 12);
        assert!(record1.is_valid());

        // Second record
        let record2 = reader.next_record()?.expect("Should have second record");
        assert_eq!(record2.id, "READ2");
        assert_eq!(record2.description, None);
        assert_eq!(record2.sequence, "GGGGCCCCAAAA");
        assert_eq!(record2.quality, "HHHHHHHHHHHH");

        // End of file
        assert!(reader.next_record()?.is_none());

        Ok(())
    }

    #[test]
    fn test_phred_scores() {
        let record = FastqRecord {
            id: "test".to_string(),
            description: None,
            sequence: "ACGT".to_string(),
            quality: "!5IA".to_string(), // Phred+33: 0, 20, 40, 32
        };

        let scores = record.phred_scores(QualityEncoding::Phred33);
        assert_eq!(scores, vec![0, 20, 40, 32]);

        let mean = record.mean_quality();
        assert!((mean - 23.0).abs() < 0.1);
    }

    #[test]
    fn test_gc_content() {
        let record = FastqRecord {
            id: "test".to_string(),
            description: None,
            sequence: "ACGTACGT".to_string(), // 4 GC out of 8 = 50%
            quality: "IIIIIIII".to_string(),
        };

        assert_eq!(record.gc_content(), 0.5);
    }

    #[test]
    fn test_invalid_record() {
        let fastq_data = "@READ1\n\
                          ACGT\n\
                          +\n\
                          III\n"; // Quality too short!

        let mut temp_file = NamedTempFile::new().unwrap();
        temp_file.write_all(fastq_data.as_bytes()).unwrap();
        temp_file.flush().unwrap();

        let mut reader = FastqReader::from_path(temp_file.path()).unwrap();
        let result = reader.next_record();
        assert!(result.is_err());
    }

    #[test]
    fn test_quality_encoding_detection() {
        // Phred+33 has characters < 64
        assert_eq!(
            QualityEncoding::detect("!5IA"),
            QualityEncoding::Phred33
        );

        // Phred+64 has no characters < 64
        assert_eq!(
            QualityEncoding::detect("@@@@"),
            QualityEncoding::Phred33 // Defaults to Phred33
        );
    }

    #[test]
    fn test_missing_separator() {
        let fastq_data = "@READ1\n\
                          ACGT\n\
                          -\n\
                          IIII\n"; // Wrong separator!

        let mut temp_file = NamedTempFile::new().unwrap();
        temp_file.write_all(fastq_data.as_bytes()).unwrap();
        temp_file.flush().unwrap();

        let mut reader = FastqReader::from_path(temp_file.path()).unwrap();
        let result = reader.next_record();
        assert!(result.is_err());
    }
}