A fast program for classifying and separating paired-end (PE) and single-end (SE) FASTQ after removing adapters via adapter-index sequences.
- Anchored 5' adapter demultiplexing — every sequence within
--error-tolerancesubstitutions of a known adapter is precomputed into a lookup table, so classifying a read is one hash probe instead of a scan over every candidate adapter. cutadapt builds the same structure for anchored adapters and describes it as purely a speed optimisation; the same holds here. - 3' adapter trimming with an indel-tolerant fallback (banded edit distance).
- Index-sequence rescan inside the located adapter region, error tolerant.
- PE merging for pairs that read through into each other.
- Quality control — low-quality, N/non-ACGT and poly-A/C/G/T end trimming.
- Compressed input — plain, gzip, bzip2 and zstd, detected from magic bytes and read
as multi-member/multi-frame streams. BGZF (what
bgzipwrites) is recognised from its per-member size field and inflated block-parallel, which brings compressed input to within a hair of plain-file speed. - Multi-threaded, SIMD-accelerated (AVX2 / NEON), with a bounded memory budget.
# install Rust
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh
# or update; SegPE needs the toolchain named in Cargo.toml's rust-version
rustup update stable
cd SegPE
cargo build --releaseThere are no feature flags to pick per CPU — AVX2 and NEON kernels are chosen at runtime.
A C compiler is needed, because mimalloc compiles its own C sources: gcc/clang on
Linux, xcode-select --install on macOS, the MSVC build tools on Windows. Nothing else
does — gzip, bzip2 and zstd all use pure-Rust backends, and cmake is not involved.
On Apple Silicon, build natively for aarch64-apple-darwin (an x86-64 build under
Rosetta loses the vector path) and pass --max-mem explicitly, since the
available-memory query is implemented for Linux and Windows only. See INSTALL section
0x03.
Prebuilt binaries for Linux x86-64, Linux aarch64 and Windows x86-64 are in bin/;
build macOS from source.
cargo build --release --features gpuAdds a wgpu compute backend covering Vulkan (Linux, Windows), DX12 (Windows) and Metal (macOS). See Do I want the GPU? — for typical runs the answer is no.
# paired-end
RUST_LOG=INFO segpe --five-art-fa data/5_art.fa --three-art-fa data/3_art.fa \
--five-idx-fa data/idx.fa \
--pe1-fastq data/read1.fq.gz --pe2-fastq data/read2.fq.gz \
-n 16 -o ./output-pe
# single-end
RUST_LOG=INFO segpe --five-art-fa data/5_art.fa --three-art-fa data/3_art.fa \
--five-idx-fa data/idx.fa --pe1-fastq data/reads.fastq \
-n 16 -o ./output-se
# with quality control and PE merging
RUST_LOG=INFO segpe --five-art-fa data/5_art.fa --three-art-fa data/3_art.fa \
--five-idx-fa data/idx.fa \
--pe1-fastq data/read1.fq.gz --pe2-fastq data/read2.fq.gz \
--merge-pe --qual-trim 8 --n-trim --poly-trim 6 --length-offset 50 \
-n 16 -o ./output-peRun segpe --help for the full option list. Options added in v0.2:
| Option | Meaning |
|---|---|
--ambiguous <report|separate|unclassified> |
What to do when a read's 5' adapter is equally close to adapters carrying different index labels, so the assignment is a tie-break rather than evidence. report (default) keeps the existing behaviour and counts them; separate writes them to their own _ambiguous files; unclassified treats them as unmatched. All three log the totals. See Ambiguous demultiplexing. |
--max-mem <MiB> |
Ceiling on RAM held in flight. 0 (default) derives one from what the OS reports as available. The ceiling is honoured: queue depths and, if necessary, the block size shrink to fit. |
--device <auto|cpu|gpu> |
Compute device for anchored 5' matching. auto only reaches for the GPU when the adapter set is too large to index on the CPU. Requires --features gpu. |
2.5M read pairs (1.67 GB of FASTQ), 384 anchored 16 bp adapters, 196 index sequences. "baseline" is the v0.1.16 code on the same data and the same machine.
Speed, plain input
| Platform | baseline | v0.2 | speedup |
|---|---|---|---|
| Windows x86-64, 32 threads, NVMe SSD | 2.31 s | 0.68 s | 3.4× |
| Linux x86-64 (WSL2), 32 threads | 7.85 s | 0.41 s | 19.1× |
| Linux aarch64 (GB10, DGX Spark), 20 threads | 6.60 s | 0.34 s | 19.4× |
| macOS aarch64 (MacBook Air M5), 10 cores | — | 0.33 s | — |
The fanless 10-core M5 matches the 20-core GB10 and beats the 32-thread desktop: this workload is bound by memory bandwidth and allocator behaviour far more than by core count. Its output digest is identical to the other three platforms.
Speed, compressed input (Linux x86-64, 32 threads)
| Container | baseline | v0.2 |
|---|---|---|
| gzip (single member) | 8.51 s | 1.45 s |
BGZF (bgzip) |
8.5 s¹ | 0.48 s |
| — plain file, for reference | 7.85 s | 0.41 s |
¹ the baseline had no BGZF path; it read the file as ordinary gzip.
BGZF stores each member's compressed size in its header, so members can be located
without inflating and then inflated concurrently — which is why it lands within 17% of
an uncompressed file. Ordinary gzip carries no such index; locating its members requires
decompressing them, so that path stays single-threaded per file. If you control how
your FASTQ is compressed, use bgzip.
Windows, gzip input: 5.23 s → 1.68 s. The gap between Windows and Linux is the system allocator, not the code: glibc's malloc serialises this workload's allocations on its arena locks and stops scaling past about four threads, which the Windows heap does not. v0.2 ships mimalloc, so both converge.
Against cutadapt — same task (384 anchored adapters, one error, no indels), 2.5M reads, idle 32-thread x86-64 host:
| Tool | Throughput |
|---|---|
| cutadapt 5.2 | 1.95M reads/s |
| SegPE v0.2 | 7.43M reads/s (3.8×) |
Classifications were compared read by read, mapping cutadapt's adapter names through the same adapter → index rule SegPE uses:
| Reads | |
|---|---|
| Identical label | 245,586 (98.23%) |
| SegPE labelled, cutadapt did not | 4,414 |
| cutadapt labelled, SegPE did not | 0 |
| Both labelled, different label | 0 |
SegPE is a strict superset here: it never contradicts cutadapt, and the extra assignments come from routes cutadapt was not asked to use (the 3' index and the indel-tolerant fallback).
| Issue | Before | After |
|---|---|---|
Multi-member gzip (bgzip, pigz, concatenated .gz) |
Only the first member was read — a 250k-read two-member file yielded 125k reads, silently | All members read |
| 5' error-tolerant matching | Unreachable: the seed map was keyed by index ids but looked up in the adapter-id map, so the branch never ran. --error-tolerance and --seed-len had no effect on 5' classification |
Works; 7,198 more reads per 250k are classified |
--error-tolerance semantics |
A score cutoff of len - e, which with match +1 / mismatch −1 admits fewer than e/2 mismatches — at the default e = 1, only exact matches |
Edit distance ≤ e, matching what the flag name says and what cutadapt means by -e |
| Write failures | Logged at error! but the run still reported success and exited 0, so a full disk truncated output silently |
Propagated to the exit code |
Output vs --batch |
Varied with batch size | Identical for any --batch |
| Memory | Grew with --batch × --num-threads, unbounded |
Bounded by --max-mem |
| Ambiguous demultiplexing | Resolved silently by tie-break | Counted and warned about; --ambiguous can divert or reject them |
Output is byte-identical across Windows x86-64, Linux x86-64, Linux aarch64 and macOS aarch64.
Correcting e errors unambiguously requires the adapters to be at least 2e + 1 apart.
Barcode sets often are not, and when they are not, some read prefixes sit equally close to
adapters belonging to different libraries. Both SegPE and cutadapt then take the first
equally-good match — deterministic, but a coin toss rather than evidence.
SegPE now measures this instead of hiding it. On the 384-adapter set in test/data,
whose adapters are only 2 apart at the closest, at the default --error-tolerance 1:
WARN 1285 of 15094 indexed variants (8.51%) are equally close to adapters carrying
different index labels, so this adapter set cannot correct 1 error(s)
unambiguously; reads landing on them are resolved by tie-break.
WARN 838 reads (0.3352%) matched two or more adapters equally well but with different
index labels; kept (first equally-close adapter wins).
0.34% sounds small until you look at what it produces: of the 59 output libraries in
that run, 26 exist only because of tie-breaks — each holding 2 to 10 reads, every one
of them a guess. Those are barcodes that were most likely never on the machine, picking
up stray reads through single-substitution collisions with real ones. Run with
--ambiguous separate and the output drops to the 33 libraries that are actually
supported by evidence, with the guesses set aside for inspection.
The startup warning fires whenever the adapter set cannot correct the requested number of errors unambiguously, which is worth knowing before the run rather than after.
Probably not. The GPU backend is implemented, correct and verified against the CPU path on both an RTX 4090 and a GB10, but end to end it is slower:
| Platform | CPU | GPU |
|---|---|---|
| RTX 4090 (Vulkan) | 0.80 s | 1.19 s |
| GB10 (Vulkan) | 0.50 s | 3.60 s |
The kernel itself is far faster — brute-force nearest-adapter search over 1M reads runs 26–443× quicker on a 4090 than on one CPU core, results asserted identical. But that computation is not the bottleneck: this pipeline is dominated by parsing, buffer copies and IO, and when the CPU can build its variant index it never performs the search at all.
--device auto therefore stays on the CPU unless the adapter set is too large to index
(for 384 × 16 bp adapters that happens at --error-tolerance 4), which is the regime the
GPU exists for.
Files are named {prefix}_{index-name}.fastq, with _merged appended for merged pairs
and unclassified for reads that matched nothing. With --trim-name, trimming offsets
are appended to the read name.
--outdir is added to,
not replaced.
0.2.0, build-260728
Licensed under either of Apache License, Version 2.0 or MIT license at your option. Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in this crate by you, as defined in the Apache-2.0 license, shall be dual licensed as above, without any additional terms or conditions.