# fastx
[](https://github.com/ScioFuturum/fastx/actions/workflows/ci.yml)
[](https://crates.io/crates/fastx-io)
[](https://docs.rs/fastx-io)
[](#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.3"
```
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")?;
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.3", 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.
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 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:
| read FASTQ, 150 bp reads, `read_into` | ~1.6 GiB/s |
| read FASTQ, 150 bp reads, iterator (allocates per record) | ~590 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:
Measured floors on the same file, so the gap to each is known rather than
guessed — all in one process, so conditions are identical:
| `read()` into a 128 KiB buffer, bytes untouched | ~17 ms | ~5000 MB/s |
| + `memchr` over every byte | ~19 ms | ~4400 MB/s |
| + copying each line out — the floor for owned records | ~30 ms | ~2840 MB/s |
| `read_into` — owned records | ~45 ms | ~1900 MB/s |
| **`read_ref` — borrowed records** | **~30 ms** | **~2860 MB/s** |
| `memchr` over a buffer already in RAM (no I/O) | ~6 ms | ~13000 MB/s |
A full pass with `read_ref` now costs about 30 ms, of which roughly 17 ms is the
`read()` itself. The parser adds ~13 ms on top of that, spread thin — bisecting it
found no hotspot left worth naming: header splitting, for instance, measures as
noise. This is the end of the road without either `unsafe` or giving up on
streaming.
`read_ref` is within 1.65× of "read the bytes and find the newlines", so most of
what is left is the `read()` itself — and that cost is the copy out of the page
cache, not the syscalls: buffer sizes from 256 KiB to 512 KiB win about 1.5% over
the 128 KiB default, and 16 MiB is measurably worse.
If you want the in-memory number, read the file yourself and parse the slice —
`FastxReader::new(&bytes[..])` measured 23.9 ms. That is the whole trick, and it
needs no `unsafe`; it just costs you the memory.
### Two optimisations that failed
Recorded because a plausible optimisation that does not work is more useful to
the next person than one that does.
**A newline index instead of per-line `memchr`: 4% slower.** The
micro-benchmark that predicted +18% only *counted* newlines. The real parser has
to materialise a million `u32` offsets into a `Vec` and read them back, which
costs more than the `memchr` calls it saves.
**Memory mapping: 1.8–2.5× slower.** The received wisdom is that `mmap` beats
`read()`, and the ceiling looked like 1.68×. Measured on Windows, one pass per
process, interleaved: ~35 ms for `read()` against 62–89 ms for a mapping. The
cause is soft page faults — 85 MB is some 21 000 pages, and each fault costs
microseconds here. Linux fields them far more cheaply and has `MAP_POPULATE`, so
the result may well invert there; it has not been measured, so no claim is made.
Either way the crate keeps `#![forbid(unsafe_code)]`, since mapping would have
traded that for a regression on the one platform actually measured.
The rows below are end-to-end through the CLI, so each includes about 90 ms of
process startup on this machine — worth subtracting before comparing them with
the in-process figures above.
| parse only | ~0.15 s | — (already ahead of storage) |
| parse + `SeqStats` | ~0.26 s | — |
| **read BGZF** | ~0.95 s | **~0.56 s** (1.7×) |
| **write BGZF, level 6** | ~9.4 s | **~2.7 s** (3.4×) |
| **write BGZF, level 6 + `libdeflate`** | — | **~1.5 s** (6.3×, and 6% smaller) |
| **write BGZF, level 1** | ~1.6 s | **~0.72 s** |
A note on method, because it changes what these numbers mean: the *first*
measurement in a session on this laptop is consistently the fastest, by up to
35%, so an earlier version of this table understated compression by exactly that
much. Every pair above was measured by interleaving the two variants in one run,
which is the only way to get a ratio worth quoting. Absolute values will differ
on your hardware; the ratios should not.
Compression scales to about 3.4× and then flattens — the increments (1→2 gives
1.65×, 4→8 only 1.15×) are the signature of four physical cores behind eight
logical ones, since deflate gains little from hyperthreading. A machine with real
cores will go further.
Three conclusions worth acting on:
1. **The parser is not your bottleneck.** At ~550 MiB/s it is ahead of most disks
and nearly all networks. Optimising it further buys nothing real.
2. **Compression is the bottleneck, and the level is the first lever.** Level 1
is several times faster than level 6 for 8% larger output — for anything
another pipeline stage will read straight back, use `-l 1` or
`WriterBuilder::level(CompressionLevel::FAST)`.
3. **Cores are the second lever, and they are free.** BGZF blocks are
independent, so batching them across cores is on by default with the
`parallel` feature and produces byte-identical output. `-@ N` on the CLI or
`WriterBuilder::blocks_per_batch` tunes it; 1 forces single-threaded.
4. **The deflate implementation is the third, and the largest single one.**
`--features libdeflate` nearly halves compression time again and gives up
nothing: 1.95× faster with 6% smaller output, still plain gzip.
Together, levels 2 and 3 take level-6 compression of this file from ~9.4 s to
~1.5 s. What is left is dominated by deflate itself, which is why pipelining the
parse stage behind it is not worth building: parsing and I/O are ~440 ms of a
~2.7 s run, so hiding them entirely would buy under 17%, against 2× for changing
the compressor.
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.3"
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
| `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 |
| `libdeflate` | no | swaps BGZF's deflate for libdeflate — measured 1.95× faster *and* 6% smaller. Same C compiler requirement |
| `full` | no | all of the above |
`libdeflate` is the rare optimisation that costs nothing in output quality: it
compresses faster and tighter, and the result is still ordinary gzip that any
tool reads. Files written by either backend are byte-different but mutually
readable, so switching is safe mid-pipeline.
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
| 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) |
### Measured against BioPython
Same 81 MiB FASTQ (250 000 × 150 bp), same machine, same tasks, best of three,
process startup excluded. All four implementations returned identical answers —
250 000 records, 37 500 000 bases, GC 49.99% — so this is like for like.
BioPython 1.88 on CPython 3.14.
| count records + bases | 2678 ms | 893 ms | 555 ms | **88 ms** |
| GC content | 8838 ms | 1022 ms | — | **136 ms** |
| reverse complement | 3038 ms | 1134 ms | — | **112 ms** |
| count from gzip | — | 1648 ms | — | **218 ms** |
So 6–10× faster than the best Python approach and 20–65× faster than the
idiomatic one. Two things in that table deserve saying out loud rather than being
quietly enjoyed:
- **Most of BioPython's cost is its object model, not its parsing.** Its own
`FastqGeneralIterator` — which yields plain strings — is 3× faster than
`SeqIO.parse`, and `str.count("G")` is 8.6× faster than `gc_fraction`. If you
are staying in Python, that alone is most of the available win.
- **A hand-written `readline` loop beats BioPython at counting.** 555 ms against
893 ms. For genuinely trivial passes, the library is overhead.
`fastx`'s own allocating iterator, the convenient one, takes 169 ms — still 3×
faster than the plain loop, so you do not have to reach for `read_into` to come
out ahead.
Startup matters in the other direction, for scripts over small files: importing
BioPython costs about 650 ms before any work happens, so a 3-record file takes
~770 ms end to end against ~90 ms for the `fastx` CLI. If you run a tool once per
sample across a few thousand samples, that difference is the whole runtime.
## 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.