gwseq_io/

lib.rs

//! # gwseq-io
//!
//! Reading and writing bigWig, bigBed, BAM, CRAM and HiC files.
//!
//! See `ARCHITECTURE.md` at the workspace root for the module map and the
//! reasoning behind the layering.
//!
//! ## Shape
//!
//! Four readers and one writer, each opened on a path or a URL:
//!
//! ```no_run
//! # fn main() -> Result<(), Box<dyn std::error::Error>> {
//! use gwseq_io::bbi::{ValuesRequest, Zoom};
//! use gwseq_io::genomic::Locs;
//! use gwseq_io::{open, Reader};
//!
//! let chr_ids = vec!["chr1".to_string(), "chr2".to_string()];
//! let starts = vec![1_000_000, 2_000_000];
//! let ends = vec![1_010_000, 2_010_000];
//!
//! let Reader::Bbi(bw) = open("track.bigwig", Default::default())? else {
//!     panic!("not a bigwig")
//! };
//! let values = bw.read_values(
//!     &ValuesRequest::new(Locs::spans(&chr_ids, &starts, &ends)?)
//!         .bin_size(10.0)
//!         .zoom(Zoom::Auto),
//! )?; // Array2<f32>, (loci, bins)
//! # Ok(())
//! # }
//! ```
//!
//! Requests are builders because the API they mirror has up to fourteen
//! defaulted keyword arguments; a builder is where each default is written down
//! once, for the Python layer, the CLI and Rust callers alike.
//!
//! ## Threads and handles
//!
//! A reader owns `parallel` worker threads and one file handle for its
//! lifetime. `close()` gives both back and is idempotent; a closed reader
//! fails every read with [`Error::Closed`], while the headers it read at open
//! stay available. A reader that is never closed gives everything back when it
//! is dropped.

#![forbid(unsafe_code)]

// The crate README is the crates.io landing page, so its examples are compiled
// and run as doctests rather than left to rot. Not rendered into these docs —
// it says the same things this module does, for a reader who has not arrived
// here yet.
#[cfg(doctest)]
#[doc = include_str!("../README.md")]
struct ReadmeExamples;

pub mod align;
pub mod arrays;
pub mod bam;
pub mod bbi;
pub mod bytes;
pub mod cram;
pub mod error;
#[cfg(test)]
mod fuzz;
pub mod genomes;
pub mod genomic;
pub mod hic;
#[cfg(feature = "npz")]
pub mod npz;
pub mod parallel;
pub mod progress;
pub mod source;

pub use align::{AlignmentWalk, Alignments};
pub use error::{Error, Result};
pub use genomes::get_chr_sizes;

/// What [`open`] returns. The format is sniffed from the file's magic number,
/// not from its extension.
///
/// The variants differ in size — a `HiCReader` carries three memo tables a
/// `BbiReader` does not — and boxing them is not worth the ergonomics: one of
/// these is constructed per file opened and then held for the reader's life,
/// so the padding is paid once against a file handle and a thread pool.
#[allow(clippy::large_enum_variant)]
pub enum Reader {
    Bbi(bbi::BbiReader),
    Bam(bam::BamReader),
    Cram(cram::CramReader),
    HiC(hic::HiCReader),
}

/// Everything [`open`] takes beyond the path.
///
/// `None` on a buffer field means "recommended", which is what `-1` spells in
/// the Python API — and the recommendation differs by source: 32 KiB blocks for
/// a local file, 1 MiB for a URL, 128 blocks either way.
#[derive(Debug, Clone)]
pub struct OpenOptions {
    /// Worker threads. Zero or less means one per core, capped at 12.
    pub parallel: i64,
    /// bigWig only: scaling factor for automatic zoom selection.
    pub zoom_correction: f64,
    pub block_size: Option<u64>,
    pub max_blocks: Option<usize>,
    /// BAM and CRAM. Defaults to the file's path with `.bai` or `.crai`
    /// appended.
    pub index_path: Option<String>,
    /// CRAM only: a FASTA holding the reference the file was compressed
    /// against. `None` looks at the `UR` field of the header's `@SQ` lines and
    /// then at `REF_CACHE`/`REF_PATH`; a CRAM that finds none still opens and
    /// reads everything but `sequence`.
    pub reference: Option<String>,
}

impl Default for OpenOptions {
    fn default() -> Self {
        Self {
            parallel: -1,
            zoom_correction: 1.0 / 3.0,
            block_size: None,
            max_blocks: None,
            index_path: None,
            reference: None,
        }
    }
}

/// What a file's first bytes say it is.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum FileKind {
    BigWig,
    BigBed,
    Bam,
    Cram,
    HiC,
}

/// Sniff the format of an already-open source.
///
/// A BAM carries its magic *behind* the gzip one, being a BGZF file, so that
/// read is only reached once the gzip magic matches. A CRAM says `CRAM` in
/// clear at offset zero, so it is settled before any of that. A byte-swapped bbi magic
/// is recognised and refused — these readers do not swap rather than swapping
/// silently, and saying which it is beats "unrecognised file".
pub fn sniff_source(source: &dyn source::ByteSource) -> Result<FileKind> {
    let head = source.read_at(0, 4)?;
    if head.len() < 4 {
        return Err(Error::format(
            source.path(),
            "file is too short to carry a magic number",
        ));
    }
    if &head[..3] == b"HIC" {
        return Ok(FileKind::HiC);
    }
    if head[..4] == cram::container::MAGIC {
        return Ok(FileKind::Cram);
    }
    let magic = u32::from_le_bytes([head[0], head[1], head[2], head[3]]);
    match magic {
        bbi::BIGWIG_MAGIC => return Ok(FileKind::BigWig),
        bbi::BIGBED_MAGIC => return Ok(FileKind::BigBed),
        bbi::header::BIGWIG_MAGIC_SWAPPED | bbi::header::BIGBED_MAGIC_SWAPPED => {
            return Err(Error::format(source.path(), "incompatible endianness"))
        }
        _ => {}
    }
    if head[0] == 0x1F && head[1] == 0x8B {
        return Ok(FileKind::Bam);
    }
    Err(Error::format(
        source.path(),
        "not a bigwig, bigbed, bam, cram or hic file",
    ))
}

/// Sniff the format of a file from its magic number.
pub fn sniff(path: &str) -> Result<FileKind> {
    // A tiny cache: this reads four bytes and is thrown away.
    let source = source::open(path, Some(4096), Some(1))?;
    sniff_source(source.as_ref())
}

/// Open a file for reading, dispatching on its magic number.
///
/// The source is opened once and handed to whichever reader the magic names, so
/// sniffing does not cost a second open — which over HTTP would be a second
/// round trip.
pub fn open(path: &str, options: OpenOptions) -> Result<Reader> {
    let source = source::open(path, options.block_size, options.max_blocks)?;
    match sniff_source(source.as_ref())? {
        FileKind::BigWig | FileKind::BigBed => Ok(Reader::Bbi(bbi::BbiReader::from_source(
            source,
            path,
            options.parallel,
            options.zoom_correction,
        )?)),
        FileKind::Bam => Ok(Reader::Bam(bam::BamReader::from_source(
            source,
            path,
            options.index_path.as_deref(),
            options.parallel,
            options.block_size,
            options.max_blocks,
        )?)),
        FileKind::Cram => Ok(Reader::Cram(cram::CramReader::from_source(
            source,
            path,
            options.index_path.as_deref(),
            options.reference.as_deref(),
            options.parallel,
            options.block_size,
            options.max_blocks,
        )?)),
        FileKind::HiC => Ok(Reader::HiC(hic::HiCReader::from_source(
            source,
            path,
            options.parallel,
        )?)),
    }
}