genomicframe_core/formats/sam/

reader.rs

//! SAM (Sequence Alignment/Map) format reader
//!
//! SAM files are tab-delimited text files containing alignment information.
//! Each line represents one alignment record with 11 mandatory fields plus optional tags.
//!
//! # Design Philosophy
//!
//! - **Streaming by default**: Records are never buffered unless explicitly requested
//! - **O(1) memory**: Process multi-GB files with constant memory
//! - **Lazy evaluation**: Parse alignments on-demand as you iterate
//! - **Text-based**: Human readable (unlike BAM's binary format)
//!
//! # SAM Format Specification
//!
//! Mandatory fields (tab-separated):
//! 1. QNAME: Query template name
//! 2. FLAG: Bitwise flags
//! 3. RNAME: Reference sequence name
//! 4. POS: 1-based leftmost mapping position
//! 5. MAPQ: Mapping quality
//! 6. CIGAR: CIGAR string
//! 7. RNEXT: Reference name of mate/next read
//! 8. PNEXT: Position of mate/next read
//! 9. TLEN: Template length
//! 10. SEQ: Segment sequence
//! 11. QUAL: ASCII Phred quality scores
//! 12+: Optional tags (TAG:TYPE:VALUE)
//!
//! # Examples
//!
//! ```no_run
//! use genomicframe_core::formats::sam::SamReader;
//! use genomicframe_core::core::GenomicRecordIterator;
//!
//! // Streaming: O(1) memory, processes one alignment at a time
//! let mut reader = SamReader::from_path("alignments.sam")?;
//!
//! // Read the header first (required)
//! let header = reader.read_header()?;
//! println!("References: {}", header.references.len());
//!
//! // Iterate through alignments
//! while let Some(record) = reader.next_record()? {
//!     if !record.is_unmapped() && record.mapq >= 30 {
//!         println!("{}: {} at {}", record.qname, record.cigar_string(), record.pos);
//!     }
//! }
//! # Ok::<(), genomicframe_core::error::Error>(())
//! ```

use crate::core::{GenomicReader, GenomicRecordIterator};
use crate::error::{Error, Result};
use crate::formats::bam::{CigarOp, RefSeq, TagValue};
use crate::io::Compression;
use flate2::read::MultiGzDecoder;
use std::collections::HashMap;
use std::fs::File;
use std::io::{BufRead, BufReader};
use std::path::Path;


/// SAM file header
#[derive(Debug, Clone)]
pub struct SamHeader {
    /// Header lines (e.g., @HD, @SQ, @RG, @PG, @CO)
    pub lines: Vec<String>,
    /// Reference sequences (from @SQ lines)
    pub references: Vec<RefSeq>,
    /// Reference name to index mapping (for fast lookup)
    ref_index: HashMap<String, usize>,
}

impl SamHeader {
    /// Create a new empty SAM header
    pub fn new() -> Self {
        Self {
            lines: Vec::new(),
            references: Vec::new(),
            ref_index: HashMap::new(),
        }
    }

    /// Get reference index by name
    pub fn ref_index(&self, name: &str) -> Option<usize> {
        self.ref_index.get(name).copied()
    }

    /// Get reference name by index
    pub fn ref_name(&self, idx: usize) -> Option<&str> {
        self.references.get(idx).map(|r| r.name.as_str())
    }
}

impl Default for SamHeader {
    fn default() -> Self {
        Self::new()
    }
}

/// SAM alignment record
///
/// This struct is identical to BAM's BamRecord but parsed from text instead of binary.
#[derive(Debug, Clone)]
pub struct SamRecord {
    /// Query template name
    pub qname: String,
    /// Bitwise FLAG
    pub flag: u16,
    /// Reference sequence name
    pub rname: String,
    /// 1-based leftmost mapping position (0 for unmapped)
    pub pos: i32,
    /// Mapping quality (255 if unavailable)
    pub mapq: u8,
    /// CIGAR string
    pub cigar: Vec<CigarOp>,
    /// Reference name of mate/next read
    pub rnext: String,
    /// Position of mate/next read (0 if unavailable)
    pub pnext: i32,
    /// Template length
    pub tlen: i32,
    /// Segment sequence (empty if unavailable)
    pub seq: Vec<u8>,
    /// Phred-scaled base qualities (empty if unavailable)
    pub qual: Vec<u8>,
    /// Optional tags
    pub tags: HashMap<String, TagValue>,
}

impl SamRecord {
    /// Check if read is paired
    pub fn is_paired(&self) -> bool {
        (self.flag & 0x1) != 0
    }

    /// Check if read is mapped in proper pair
    pub fn is_proper_pair(&self) -> bool {
        (self.flag & 0x2) != 0
    }

    /// Check if read is unmapped
    pub fn is_unmapped(&self) -> bool {
        (self.flag & 0x4) != 0
    }

    /// Check if mate is unmapped
    pub fn is_mate_unmapped(&self) -> bool {
        (self.flag & 0x8) != 0
    }

    /// Check if read is reverse strand
    pub fn is_reverse(&self) -> bool {
        (self.flag & 0x10) != 0
    }

    /// Check if mate is reverse strand
    pub fn is_mate_reverse(&self) -> bool {
        (self.flag & 0x20) != 0
    }

    /// Check if read is first in pair
    pub fn is_first_in_pair(&self) -> bool {
        (self.flag & 0x40) != 0
    }

    /// Check if read is second in pair
    pub fn is_second_in_pair(&self) -> bool {
        (self.flag & 0x80) != 0
    }

    /// Check if alignment is secondary
    pub fn is_secondary(&self) -> bool {
        (self.flag & 0x100) != 0
    }

    /// Check if read fails quality checks
    pub fn is_qc_fail(&self) -> bool {
        (self.flag & 0x200) != 0
    }

    /// Check if read is a PCR or optical duplicate
    pub fn is_duplicate(&self) -> bool {
        (self.flag & 0x400) != 0
    }

    /// Check if alignment is supplementary
    pub fn is_supplementary(&self) -> bool {
        (self.flag & 0x800) != 0
    }

    /// Get CIGAR string representation
    pub fn cigar_string(&self) -> String {
        if self.cigar.is_empty() {
            "*".to_string()
        } else {
            self.cigar.iter().map(|op| op.to_string()).collect()
        }
    }

    /// Get sequence as string
    pub fn seq_string(&self) -> String {
        if self.seq.is_empty() {
            "*".to_string()
        } else {
            String::from_utf8_lossy(&self.seq).to_string()
        }
    }
}

/// SAM file reader (streaming, memory-efficient)
pub struct SamReader {
    reader: Box<dyn BufRead>,
    header: Option<SamHeader>,
    line_buffer: String,
    line_number: usize,
}

impl SamReader {
    /// Open a SAM file from a path - handles .sam and .sam.gz automatically
    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 => Box::new(BufReader::new(MultiGzDecoder::new(file))),
            _ => Box::new(BufReader::new(file)),
        };

        Ok(Self {
            reader,
            header: None,
            line_buffer: String::with_capacity(512),
            line_number: 0,
        })
    }

    /// Create a new SAM reader from any BufRead source
    pub fn new<R: BufRead + 'static>(reader: R) -> Self {
        Self {
            reader: Box::new(reader),
            header: None,
            line_buffer: String::with_capacity(512),
            line_number: 0,
        }
    }

    /// Read and parse the SAM header
    ///
    /// This must be called before reading records.
    /// Header lines start with '@', record lines don't.
    pub fn read_header(&mut self) -> Result<&SamHeader> {
        if self.header.is_some() {
            return Ok(self.header.as_ref().unwrap());
        }

        let mut header = SamHeader::new();

        loop {
            self.line_buffer.clear();
            let bytes = self.reader.read_line(&mut self.line_buffer)?;
            if bytes == 0 {
                // Empty file
                break;
            }

            self.line_number += 1;

            let line = self.line_buffer.trim();

            if line.is_empty() {
                continue;
            }

            if !line.starts_with('@') {
                // First non-header line - put it back by not consuming
                // We'll re-read it in next_record()
                self.line_number -= 1;
                break;
            }

            header.lines.push(line.to_string());

            // Parse @SQ (reference sequence) lines
            if line.starts_with("@SQ") {
                let mut name = None;
                let mut length = None;

                // Parse tab-separated fields
                for field in line.split('\t').skip(1) {
                    if let Some((key, value)) = field.split_once(':') {
                        match key {
                            "SN" => name = Some(value.to_string()),
                            "LN" => length = value.parse::<u32>().ok(),
                            _ => {}
                        }
                    }
                }

                if let (Some(name), Some(length)) = (name, length) {
                    let idx = header.references.len();
                    header.ref_index.insert(name.clone(), idx);
                    header.references.push(RefSeq { name, length });
                }
            }
        }

        self.header = Some(header);
        Ok(self.header.as_ref().unwrap())
    }

    /// Get the header (must call read_header first)
    pub fn header(&self) -> Option<&SamHeader> {
        self.header.as_ref()
    }

    /// Parse CIGAR string into operations (delegates to shared implementation)
    fn parse_cigar(cigar_str: &str) -> Result<Vec<CigarOp>> {
        CigarOp::parse_cigar_string(cigar_str)
    }

    /// Parse a SAM record line
    fn parse_record(&mut self, line: &str) -> Result<SamRecord> {
        let fields: Vec<&str> = line.split('\t').collect();

        if fields.len() < 11 {
            return Err(Error::Parse(format!(
                "Line {}: SAM record must have at least 11 fields, got {}",
                self.line_number,
                fields.len()
            )));
        }

        // Parse mandatory fields
        let qname = fields[0].to_string();
        let flag = fields[1].parse::<u16>().map_err(|_| {
            Error::Parse(format!("Line {}: Invalid FLAG: {}", self.line_number, fields[1]))
        })?;
        let rname = fields[2].to_string();
        let pos = fields[3].parse::<i32>().map_err(|_| {
            Error::Parse(format!("Line {}: Invalid POS: {}", self.line_number, fields[3]))
        })?;
        let mapq = fields[4].parse::<u8>().map_err(|_| {
            Error::Parse(format!("Line {}: Invalid MAPQ: {}", self.line_number, fields[4]))
        })?;
        let cigar = Self::parse_cigar(fields[5])?;
        let rnext = fields[6].to_string();
        let pnext = fields[7].parse::<i32>().map_err(|_| {
            Error::Parse(format!("Line {}: Invalid PNEXT: {}", self.line_number, fields[7]))
        })?;
        let tlen = fields[8].parse::<i32>().map_err(|_| {
            Error::Parse(format!("Line {}: Invalid TLEN: {}", self.line_number, fields[8]))
        })?;

        // Sequence
        let seq = if fields[9] == "*" {
            Vec::new()
        } else {
            fields[9].as_bytes().to_vec()
        };

        // Quality
        let qual = if fields[10] == "*" {
            Vec::new()
        } else {
            // Convert ASCII to Phred+33 scores
            fields[10].bytes().map(|b| b.saturating_sub(33)).collect()
        };

        // Parse optional tags
        let mut tags = HashMap::new();
        for field in fields.iter().skip(11) {
            if let Some((tag, rest)) = field.split_once(':') {
                if let Some((type_str, value)) = rest.split_once(':') {
                    let tag_value = match type_str {
                        "A" => {
                            // Single character
                            TagValue::Char(value.bytes().next().unwrap_or(0))
                        }
                        "i" => {
                            // Integer
                            TagValue::Int(value.parse().unwrap_or(0))
                        }
                        "f" => {
                            // Float
                            TagValue::Float(value.parse().unwrap_or(0.0))
                        }
                        "Z" => {
                            // String
                            TagValue::String(value.to_string())
                        }
                        "H" => {
                            // Hex string
                            TagValue::String(value.to_string())
                        }
                        "B" => {
                            // Array (skip complex parsing for now)
                            TagValue::ByteArray(value.as_bytes().to_vec())
                        }
                        _ => continue,
                    };

                    tags.insert(tag.to_string(), tag_value);
                }
            }
        }

        Ok(SamRecord {
            qname,
            flag,
            rname,
            pos,
            mapq,
            cigar,
            rnext,
            pnext,
            tlen,
            seq,
            qual,
            tags,
        })
    }

    /// Compute statistics from this reader
    pub fn compute_stats(
        mut self,
        _config: Option<crate::parallel::ParallelConfig>,
    ) -> Result<super::stats::SamStats> {
        // Ensure header is read
        self.read_header()?;

        // For now, SAM uses sequential processing (parallel SAM is complex due to text format)
        super::stats::SamStats::compute(&mut self)
    }
}

impl GenomicRecordIterator for SamReader {
    type Record = SamRecord;

    fn next_raw(&mut self) -> Result<Option<Vec<u8>>> {
        self.line_buffer.clear();
        let bytes = self.reader.read_line(&mut self.line_buffer)?;

        if bytes == 0 {
            return Ok(None);
        }

        self.line_number += 1;

        let line = self.line_buffer.trim();

        // Skip empty lines
        if line.is_empty() {
            return self.next_raw();
        }

        // Skip header lines (shouldn't happen after read_header, but be safe)
        if line.starts_with('@') {
            return self.next_raw();
        }

        Ok(Some(line.as_bytes().to_vec()))
    }

    fn next_record(&mut self) -> Result<Option<Self::Record>> {
        self.line_buffer.clear();
        let bytes = self.reader.read_line(&mut self.line_buffer)?;

        if bytes == 0 {
            return Ok(None);
        }

        self.line_number += 1;

        let line = self.line_buffer.trim().to_string(); // Clone to avoid borrow issue

        // Skip empty lines
        if line.is_empty() {
            return self.next_record();
        }

        // Skip header lines (shouldn't happen after read_header, but be safe)
        if line.starts_with('@') {
            return self.next_record();
        }

        Ok(Some(self.parse_record(&line)?))
    }
}

impl GenomicReader for SamReader {
    type Metadata = SamHeader;

    fn metadata(&self) -> &Self::Metadata {
        self.header
            .as_ref()
            .expect("Header not read yet. Call read_header() first.")
    }
}

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

    #[test]
    fn test_cigar_parsing() {
        let cigar = SamReader::parse_cigar("10M2I5D").unwrap();
        assert_eq!(cigar.len(), 3);
        assert_eq!(cigar[0], CigarOp::Match(10));
        assert_eq!(cigar[1], CigarOp::Ins(2));
        assert_eq!(cigar[2], CigarOp::Del(5));
    }

    #[test]
    fn test_cigar_star() {
        let cigar = SamReader::parse_cigar("*").unwrap();
        assert!(cigar.is_empty());
    }

    #[test]
    fn test_sam_record_flags() {
        let record = SamRecord {
            qname: "read1".to_string(),
            flag: 0x1 | 0x2, // Paired and proper pair
            rname: "chr1".to_string(),
            pos: 100,
            mapq: 60,
            cigar: vec![CigarOp::Match(100)],
            rnext: "=".to_string(),
            pnext: 200,
            tlen: 150,
            seq: b"ACGT".to_vec(),
            qual: vec![30, 30, 30, 30],
            tags: HashMap::new(),
        };

        assert!(record.is_paired());
        assert!(record.is_proper_pair());
        assert!(!record.is_unmapped());
        assert!(!record.is_duplicate());
    }
}