rsomics-fasta-utils 0.1.1

FASTA utility toolkit — count, chroms, len, revcomp, rename, tab, wrap, unique, convert, and more
Documentation
#![allow(clippy::cast_precision_loss)]

use std::io::{BufWriter, Write};
use std::path::Path;

use needletail::parse_fastx_file;
use rsomics_common::{Result, RsomicsError};

pub fn window_stats(
    input: &Path,
    output: &mut dyn Write,
    window_size: usize,
    step: usize,
) -> Result<u64> {
    if std::fs::metadata(input).is_ok_and(|m| m.len() == 0) {
        return Err(RsomicsError::InvalidInput("empty file".into()));
    }

    let mut reader = parse_fastx_file(input)
        .map_err(|e| RsomicsError::InvalidInput(format!("{}: {e}", input.display())))?;

    let mut out = BufWriter::with_capacity(256 * 1024, output);
    writeln!(out, "chrom\tstart\tend\tgc_pct\tlen").map_err(RsomicsError::Io)?;

    let mut count: u64 = 0;

    while let Some(record) = reader.next() {
        let record = record.map_err(|e| RsomicsError::InvalidInput(format!("reading: {e}")))?;
        let id = std::str::from_utf8(record.id()).unwrap_or("unknown");
        let seq = record.seq();
        let seq_len = seq.len();

        let mut pos = 0;
        while pos + window_size <= seq_len {
            let window = &seq[pos..pos + window_size];
            let gc = window
                .iter()
                .filter(|&&b| b == b'G' || b == b'g' || b == b'C' || b == b'c')
                .count();
            let gc_pct = gc as f64 / window_size as f64 * 100.0;

            writeln!(
                out,
                "{id}\t{pos}\t{}\t{gc_pct:.2}\t{window_size}",
                pos + window_size
            )
            .map_err(RsomicsError::Io)?;
            count += 1;
            pos += step;
        }
    }

    out.flush().map_err(RsomicsError::Io)?;
    Ok(count)
}