fastx-io 0.2.0

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

[![CI](https://github.com/ScioFuturum/fastx/actions/workflows/ci.yml/badge.svg)](https://github.com/ScioFuturum/fastx/actions/workflows/ci.yml)
[![Crates.io](https://img.shields.io/crates/v/fastx-io.svg)](https://crates.io/crates/fastx-io)
[![Docs](https://docs.rs/fastx-io/badge.svg)](https://docs.rs/fastx-io)
[![License](https://img.shields.io/badge/license-MIT%20OR%20Apache--2.0-blue.svg)](#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).

```toml
[dependencies]
fastx-io = "0.2"
```

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

```rust
// 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:

```rust
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:

```rust
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

```rust
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

```rust
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

```rust
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

```rust
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:

```rust
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:

```bash
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:

```rust
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:

```rust
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`)

```toml
fastx-io = { version = "0.2", features = ["parallel"] }
```

Two strategies, because they suit different problems:

```rust
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:

```toml
[dependencies]
fastx-io = "0.2"
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` |
| `zstd` | no | Zstandard in and out; pulls in the reference C library, so it needs a C compiler |
| `full` | no | all of the above |

Zstandard is worth it for intermediate files — smaller and faster than gzip at
comparable settings — but the surrounding bioinformatics toolchain does not read
it, and a zstd frame cannot be randomly accessed. A reference genome you intend
to index still wants BGZF.

`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 `I`s. 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

```bash
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](CONTRIBUTING.md) for what belongs in which test layer, and
[SECURITY.md](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 for the library and the CLI. Optional features inherit the MSRV
of what they pull in, and CI checks each at its own version rather than claiming
one number covers everything: `parallel` needs 1.80 (rayon), and `zstd` needs
1.85 (its C build chain, via `cc` → `jobserver`).

## License

Licensed under either of

- Apache License, Version 2.0 ([LICENSE-APACHE]LICENSE-APACHE)
- MIT license ([LICENSE-MIT]LICENSE-MIT)

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.