fastx-io 0.1.0

Fast, streaming FASTA/FASTQ reader and writer for bioinformatics pipelines
Documentation

fastx

CI Crates.io Docs License

Fast, streaming FASTA/FASTQ I/O for bioinformatics pipelines — a Rust replacement for the everyday sequence work you would otherwise do in BioPython or BioPerl.

  • Streaming. One growable buffer and one record at a time: a 300 GB FASTQ costs the same memory as a 300 byte one.
  • Transparent gzip, and real BGZF. Detected from the extension and from the magic bytes, so a compressed file with a misleading name still works. Compressed output is BGZF, so it stays randomly accessible with a .gzi.
  • Honest parsing. Multi-line FASTA and FASTQ, CRLF, blank lines, missing final newline, @ inside quality strings, single-line chromosomes. Malformed input gives you an error with a line number, never a silently skipped record.
  • No allocation in the hot loop. read_into reuses one record.
  • Batteries included. Reverse complement, translation, k-mers, quality trimming, N50 and friends, samtools faidx-compatible random access, optional multicore processing, and a CLI.
  • No unsafe, two small dependencies (memchr, plus flate2 for gzip).
[dependencies]
fastx-io = "0.1"

The package is fastx-io because plain fastx was already taken on crates.io by an unrelated crate. The library itself is still called fastx, so every import reads use fastx::… and the command line tool is still fastx.

Quick start

// Reading: format and gzip are detected for you.
for record in fastx::open("reads.fq.gz")? {
    let record = record?;
    println!("{}\t{}\t{:?}", record.id, record.len(), record.gc_content());
}

In a pipeline stage, reuse the record so parsing does not allocate at all:

use fastx::{FastxReader, Sequence};

let mut reader = fastx::open("reads.fq.gz")?;
let mut record = Sequence::default();
let mut bases = 0u64;
while reader.read_into(&mut record)? {
    bases += record.len() as u64;
}

Writing, with the format and compression taken from the path:

use fastx::{Alphabet, Sequence, WriterBuilder};

let mut writer = WriterBuilder::new()
    .line_width(80)            // FASTA wrapping; 0 puts each record on one line
    .validate(Alphabet::Dna)   // refuse anything but ACGTN
    .create("contigs.fa.gz")?;

writer.write_record(&Sequence::fasta("contig1", b"ACGTACGT"))?;
writer.finish()?;              // always finish(): it surfaces gzip errors

What it does

Records

use fastx::Sequence;

let read = Sequence::fastq("read1", b"ACGTTT", b"IIIII!")?;

read.len();                       // 6
read.gc_content();                // Some(0.333…)
read.base_counts();               // BaseCounts { a: 1, c: 1, g: 1, t: 3, .. }
read.mean_quality();              // error-rate-weighted, not a naive average
read.expected_errors();           // what fastp/vsearch filter on
read.reverse_complement();        // reverses the quality string too
read.translate(0, true);          // standard genetic code, stop at the first stop
read.subseq(1..4)?;               // carries the matching quality slice
read.kmers(3).count();            // 4

Quality trimming

use fastx::qual::{self, PHRED33};

let quality = b"##IIIIIIIIII###";
qual::trim_mott(quality, PHRED33, 20);                  // Mott, as in seqtk trimfq
qual::trim_sliding_window(quality, PHRED33, 4, 20.0);   // as in trimmomatic
qual::trim_ends(quality, PHRED33, 20);                  // plain end trimming
qual::fraction_at_least(quality, PHRED33, 30);          // Q30 fraction

Statistics

let stats = fastx::stats_of("assembly.fa")?;
println!("{stats}");
// records      1_243
// total bases  4_812_991
// N50          182_004
// L50          9
// GC%          41.32

SeqStats is mergeable, so per-file or per-thread accumulators combine exactly.

Random access with a .fai index

use fastx::index::{FastaIndex, IndexedFasta};

FastaIndex::build_from_path("hg38.fa")?.write_to_path("hg38.fa")?; // hg38.fa.fai

let mut reference = IndexedFasta::open("hg38.fa")?;
let promoter = reference.fetch_region("chr7", 55_019_017, 55_019_365)?; // 0-based
let same = reference.fetch_locus("chr7:55019018-55019365")?;            // 1-based

The index format is byte-compatible with samtools faidx, in both directions.

Paired-end reads

R1 and R2 have to stay in lockstep. When they drift — one file filtered without the other, or the two sorted differently — the output is still perfectly well-formed FASTQ, just silently mispaired, and nothing downstream notices. PairedReader checks the mate names on every pair, so that becomes an error on the first record instead of a wrong answer at the end:

use fastx::paired::PairedReader;

let mut reader = PairedReader::open("reads_R1.fq.gz", "reads_R2.fq.gz")?;
reader.for_each_pair(|pair| {
    assert!(pair.names_match());
    Ok(())
})?;

Interleaved input works the same way (PairedReader::open_interleaved), and a mate file that ends early is an error rather than an early stop. On the command line:

fastx interleave R1.fq.gz R2.fq.gz -o both.fq.gz

BGZF: compressed and randomly accessible

Plain gzip is one deflate stream and can never be seeked. BGZF — what bgzip writes — is valid gzip split into independent 64 KiB members, so with a .gzi block index any offset is one seek away. Compressed output from this crate is BGZF by default for exactly that reason:

let mut reference = IndexedFasta::open("hg38.fa.gz")?;   // BGZF, .fai and .gzi
let exon = reference.fetch_locus("chr7:55019018-55019365)")?;

The same .fai describes the file whether or not it is compressed, because its offsets are uncompressed positions — so fastx faidx ref.fa.gz writes a .fai samtools can use, plus the .gzi beside it.

BgzfReader implements Read + Seek in uncompressed coordinates, so it drops into anything generic over those traits, and BgzfWriter hands you the index it built while writing:

use fastx::bgzf::BgzfWriter;

let mut writer = BgzfWriter::create("reads.fq.gz")?;
writer.write_all(b"@r\nACGT\n+\nIIII\n")?;
let (_file, index) = writer.finish_with_index()?;
index.write_to_path("reads.fq.gz")?;      // reads.fq.gz.gzi

Block boundaries cost a little compression — an 81 MiB FASTQ came out 4.9% larger than plain gzip at the same level — and buy random access, which plain gzip cannot offer at any price. Pass Compression::Gzip if you want a single deflate stream anyway.

Multicore (feature parallel)

fastx-io = { version = "0.1", features = ["parallel"] }

Two strategies, because they suit different problems:

use fastx::parallel;

// 1. One parser thread, parallel work per batch. Works with gzip and stdin.
parallel::par_for_each(&mut fastx::open("reads.fq.gz")?, 4096, |record| {
    expensive_analysis(record);
    Ok(())
})?;

// 2. Split an uncompressed file into byte ranges snapped to record boundaries,
//    so parsing itself scales. FASTQ boundaries are validated, not guessed.
let stats = parallel::par_stats_file("reads.fq", 8)?;

Command line tool

cargo install fastx-io
fastx stats reads_R1.fq.gz reads_R2.fq.gz
fastx stats --json reads.fq.gz | jq .stats.n50
fastx convert reads.fq.gz -o reads.fa.gz -w 80
fastx filter --min-len 200 --min-qual 20 reads.fq -o clean.fq
fastx head -n 1000 reads.fq.gz -o subsample.fq
fastx sample -n 10000 --seed 42 reads.fq.gz -o subset.fq.gz
fastx dedup --by-seq contigs.fa -o unique.fa
fastx rc contigs.fa | fastx translate --stop-at-stop -o proteins.faa
fastx faidx hg38.fa.gz                 # writes hg38.fa.gz.fai and .gzi
fastx faidx hg38.fa.gz chr1:1-60 chr2  # extracts regions
fastx interleave R1.fq.gz R2.fq.gz -o both.fq.gz
fastx deinterleave both.fq.gz --out1 R1.fq.gz --out2 R2.fq.gz

Every command reads plain, gzipped or BGZF FASTA/FASTQ, takes - (or nothing) for stdin, and writes stdout unless given -o. stats --json emits one JSON object per line, so a run over many files pipes straight into jq. sample is seeded, so the same seed always selects the same reads.

Performance

Measured with cargo bench on a laptop (Windows 11, x86-64, --release, single thread), median of 100 samples:

benchmark throughput
read FASTQ, 150 bp reads, read_into ~1.0 GiB/s
read FASTQ, 150 bp reads, iterator (allocates per record) ~620 MiB/s
read FASTA, 10 kb records wrapped at 60 ~1.8 GiB/s
write FASTQ ~1.4 GiB/s
write FASTA wrapped at 60 ~1.8 GiB/s
reverse complement ~1.3 GiB/s
translate ~330 MiB/s

Reproduce with cargo bench; absolute numbers on your hardware will differ, but the relative cost of the iterator versus read_into is the point worth remembering — the only difference between those two rows is one allocation per record.

Where the time actually goes

End-to-end on an 81 MiB FASTQ (250 000 × 150 bp) through the CLI, which is a more honest guide to tuning than any micro-benchmark:

stage time note
parse only 0.12–0.21 s 400–700 MiB/s — faster than most storage
parse + SeqStats 0.27–0.41 s composition and quality histogram
parse + write plain ~0.25 s
read gzip +0.6–0.8 s inflate, single-threaded
write gzip, level 6 10.7–18.6 s 38.9 MB out
write gzip, level 1 2.8 s 42.2 MB out

Run-to-run variance on this laptop is wide — up to 1.7× on the gzip rows — so treat the ranges as ranges. The level 1 versus level 6 pair was measured in a single run, which makes that 6.5× ratio the trustworthy number here.

Two conclusions worth acting on:

  1. The parser is not your bottleneck. At ~700 MiB/s it is ahead of most disks and nearly all networks. Optimising it further buys nothing real.
  2. gzip compression is the bottleneck, and the level is the lever. Level 1 is 6.5× faster than level 6 for 8% larger output — for anything another pipeline stage will read straight back, use -l 1 on the CLI or WriterBuilder::level(CompressionLevel::FAST).

If gzip is still the limit after that, swap flate2's backend for zlib-ng from your own manifest — Cargo's feature unification applies it to this crate too, with no code change here:

[dependencies]
fastx-io = "0.1"
flate2 = { version = "1", features = ["zlib-ng"] }  # needs cmake and a C compiler

(We have not benchmarked that combination ourselves; zlib-ng generally reports 2–3× on inflate/deflate over the pure-Rust backend.)

Cargo features

feature default effect
gzip yes transparent gzip/BGZF reading and writing via flate2
parallel no the parallel module, via rayon
full no both of the above

default-features = false gives you a pure-Rust, dependency-light parser with memchr as the only dependency.

Design decisions worth knowing

  • Auto-detection is explicit about failure. If neither the extension nor the first byte identifies the format, you get Error::UnknownFormat rather than a guess. You can always force it with ReaderBuilder::format.
  • FASTQ is parsed by length, not by counting lines. Quality is read until it matches the sequence length, which is what makes multi-line FASTQ and @ at the start of a quality line non-issues.
  • Nothing is invented. Writing a FASTQ record without quality scores is an error, not a string of fabricated Is. Writing a FASTQ record to a FASTA writer drops the quality, which is what fq2fa should do.
  • Phred offsets are always explicit in the qual module. QualityEncoding::detect returns None when a sample is compatible with both Phred+33 and Phred+64 instead of guessing.
  • The FASTA index refuses ragged records, exactly as samtools faidx does: only the last line of a record may be shorter, and none may be longer.

Comparison

fastx BioPython SeqIO seqkit
streaming yes yes yes
gzip in/out yes with a wrapper yes
allocation-free loop yes no n/a
.fai random access yes separate index class yes
parallel parsing yes no yes
embeddable as a library yes yes no (CLI)

Development

cargo test --all-features      # 164 tests: unit, integration, CLI, property, doc
cargo clippy --all-features --all-targets
cargo fmt --check
cargo bench
cargo +nightly fuzz run reader # fuzz targets live in fuzz/

CI also runs cargo deny (licences and RUSTSEC advisories), coverage, the MSRV check on three platforms, and a short fuzzing smoke run. See CONTRIBUTING.md for what belongs in which test layer, and SECURITY.md for the threat model — in short, input files are assumed hostile, and both panics and silently wrong data count as vulnerabilities.

The MSRV is 1.74. The optional parallel feature inherits rayon's MSRV of 1.80.

License

Licensed under either of

at your option. Unless you explicitly state otherwise, any contribution you intentionally submit for inclusion in this work shall be dual licensed as above, without any additional terms or conditions.