# rsomics-seq-grep
Filter FASTA/FASTQ records by ID, full header, or sequence content — a Rust
port of `seqkit grep`.
## Install
```
cargo install rsomics-seq-grep
```
## Usage
```
# Keep records by exact ID (default target is the ID — the header up to the
# first space/tab, not the whole header line)
rsomics-seq-grep -p seq1,seq2 in.fa
# Match against the full header line instead of just the ID
rsomics-seq-grep -n -p "seq1 description text" in.fa
# 10k-ID allow-list from a file (a pattern file always overrides -p entirely)
rsomics-seq-grep -f ids.txt in.fq.gz
# Regex against ID/name (unanchored — partial match, unlike the default
# exact-match mode)
rsomics-seq-grep -r -p 'sample_[0-9]+$' in.fa
# Search the sequence itself; both strands are searched by default
rsomics-seq-grep -s -p GAATTC in.fa
# IUPAC degenerate motif (implies -s)
rsomics-seq-grep -s -d -p RYSW in.fa
# Restrict sequence search to a 1-based region (implies -s)
rsomics-seq-grep -s -R 1:12 -p ATG in.fa
# Invert: drop matching records, keep the rest
rsomics-seq-grep -f ids.txt --invert-match in.fa
# Count matches without writing records
rsomics-seq-grep -C -f ids.txt in.fa
# stdin
## Semantics (verified against `seqkit` v2.13.0, both by reading
`shenwei356/seqkit`/`shenwei356/bio` source and by black-box comparison)
- **Default target is the ID**, defined the same way as `seqkit`'s
`DefaultIDRegexp` (`^(\S+)\s?`): the header up to the first space or tab,
whichever comes first. `-n/--by-name` switches the target to the entire
header line (ID + description) instead.
- **Default matching is whole-string equality**, not substring — `-p seq1`
does **not** match a record named `seq1_dup`. This is `seqkit`'s own
documented deviation from POSIX/GNU grep. `-r/--use-regexp` switches to an
**unanchored** (partial-match) regex — the regex source is used exactly as
given, not escaped.
- `-s/--by-seq` switches the target to the sequence itself and switches
matching to **substring** search (not equality). Both the plus and minus
(reverse-complement) strand are searched by default —
`-P/--only-positive-strand` restricts to the plus strand only. Strand
search short-circuits: the first strand that hits wins, matching
`seqkit`'s `for strand := range strands { if hit { break } }`.
- `-d/--degenerate` expands an IUPAC nucleotide motif (`ACGTURYMKSWHBVDN`,
case-preserved) into the exact regex character classes from `seqkit`'s
`bio/seq/seq.go` `DegenerateBaseMapNucl` table, then matches as an
unanchored regex. Implies `-s`. Protein degenerate codes
(`DegenerateBaseMapProt`) are not implemented — `-d` here is
nucleotide-only, the overwhelmingly common use (restriction sites, primer
motifs, degenerate probes).
- `-R/--region <start:end>` restricts `-s` search to a 1-based window,
applied independently per strand (so `-R 1:12` on the minus strand means
the first 12 bases of the *reverse complement*). An out-of-range region
skips that strand rather than erroring. Negative indices count from the
end (`-12:-1` = last 12 bases); mixed-sign spans (`start<0, end>0`) are
rejected as ambiguous, matching `seqkit`. Implies `-s`.
- An **empty pattern** (`-p ""`) is legal and specific: in ID/name mode it
matches only records with an empty ID/name; in `-s` mode
(`bytes.Contains` semantics upstream) it matches only records with an
**empty sequence** rather than "matches everywhere" — `seqkit` special-cases
this rather than relying on the substring-search library's usual
empty-needle-always-matches behaviour, and so do we. In regex/degenerate
mode an empty pattern compiles to `^$` (exact-empty match) rather than the
empty regex (which would match every position).
- `-i/--ignore-case` lowercases both sides for equality/substring search and
prefixes `(?i)` for regex/degenerate search — same as upstream.
- **A pattern file (`-f`) overrides `-p` entirely** when both are given —
this is `seqkit`'s actual behaviour (verified empirically, not just from
reading the help text): the two sources are never merged. Pattern-file
lines are read verbatim except for a trailing `\r`/`\n` strip; blank lines
are skipped.
- **FASTA output line width** defaults to 60 (`-w/--line-width`, 0 = no
wrap), matching `seqkit`'s default. **FASTQ output always ignores
`-w`** — `seqkit grep` forces line width to 0 for FASTQ regardless of the
flag (`fastx.ForcelyOutputFastq` / the per-file `LineWidth = 0` override in
`grep.go`), and so do we.
- Header lines are **reproduced verbatim** — grep never reconstructs or
reflows the header, only the sequence/quality body is wrapped.
- Input: plain or gzip FASTA/FASTQ, content-sniffed (not by file extension);
`-` reads stdin. Output order always matches input order.
## Flags implemented vs. deferred
Implemented to real feature-completeness: `-p/--pattern`, `-f/--pattern-file`,
`-n/--by-name`, `-r/--use-regexp`, `-i/--ignore-case`, `-s/--by-seq`,
`-d/--degenerate`, `--invert-match`, `-R/--region`,
`-P/--only-positive-strand`, `-w/--line-width`, `-C/--count`, `-o/--output`.
Deferred — **not defined on the CLI at all**, so clap itself rejects them
with "unexpected argument" (exit 2) rather than silently ignoring them:
| `-m/--max-mismatch` | Upstream's mismatch search builds an FM-index (BWT) per record via a separate goroutine-parallel code path — a materially different algorithm, high effort for a rarely-used flag (upstream's own docs suggest alignment tools instead for large genomes). |
| `-c/--circular` | Niche (plasmid/circular-genome wraparound search); not exercised by the common FASTA/FASTQ filtering use case this crate targets. |
| `-D/--allow-duplicated-patterns` | Edge-case interaction with duplicate pattern counting; no behavioural difference for the normal (non-duplicated) pattern case this crate is built for. |
| `--delete-matched` | Upstream performance-only optimisation for very large pattern sets under `-r`; does not change output for the common case. |
| `-I/--immediate-output` | Buffering behaviour only, no output difference. |
| `-t/--seq-type`, `--alphabet-guess-seq-length`, `--id-ncbi`, `--id-regexp`, `--compress-level`, `-X/--infile-list`, `--skip-file-check` | Alphabet-guessing / ID-parsing / IO knobs outside this crate's scope; the default ID/alphabet behaviour (documented above) already matches upstream's default (`auto`/`DefaultIDRegexp`) path. |
One further documented divergence: upstream auto-detects the input
alphabet and forces `-P` (positive-strand-only) automatically for
protein/unclassified sequences, since reverse-complementing a protein
sequence is meaningless. This crate always honours `-s`'s "both strands by
default" rule regardless of alphabet — it only matters if `-s` is used
without `-P` on protein input, an unusual combination for a nucleotide-first
tool.
`-v` is not used for invert-match here (unlike upstream) because `-v` is
already `CommonFlags`' `--verbose` short flag across every `rsomics-*` tool
in this campaign; `--invert-match` is long-only by the established
conflict-resolution convention.
## Origin
Independent Rust reimplementation of `seqkit grep`, informed by:
- The seqkit paper: Shen, W. et al. *SeqKit2: A Swiss Army Knife for
Sequence and Alignment Processing.* iMeta (2024)
[doi:10.1002/imt2.191].
- Reading the upstream source (seqkit is MIT-licensed, so reading and citing
is allowed and standard practice in this project): `seqkit/cmd/grep.go`
(flag wiring, matching precedence, strand loop, region handling),
`shenwei356/bio` `seqio/fastx/reader.go` (`parseHeadIDAndDesc`,
`DefaultIDRegexp`), `seqio/fastx/records.go` (`Format`/`FormatToWriter` —
the exact FASTA wrap / FASTQ-always-unwrapped output rule),
`seq/seq.go` (`DegenerateBaseMapNucl`, `RevCom`/`Complement`), and
`seq/alphabet.go` (the `DNAredundant` complement pairing table).
- Black-box behaviour comparison against `seqkit` v2.13.0 across ID/name
exact match, regex, ignore-case, invert-match, by-seq substring, both
strands, only-positive-strand, IUPAC degenerate motifs, region
restriction, line-width variants, gzip input, stdin, pattern files, and
count-only mode (`tests/compat.rs`, frozen goldens in `tests/golden/`).
License: MIT OR Apache-2.0.
Upstream credit: [seqkit](https://github.com/shenwei356/seqkit) (MIT).
### External-dep quadrant classification
- `needletail` — Quadrant ① (pure Rust FASTA/FASTQ parser; content-sniffed
gzip via its `flate2`/`zlib-rs` backend).
- `regex` — Quadrant ① (pure Rust, SIMD-accelerated literal prefilter
internally).
- `memchr` — Quadrant ① (pure Rust, SIMD substring search for `-s` literal
matching).
- `rayon` — Quadrant ① (pure Rust parallelism; drives the by-seq path's
per-record matching across threads, honouring `-t/--threads` via
`rsomics-common`'s default-on rayon-pool wiring).
- `rsomics-common`, `rsomics-help`, `clap`, `serde` — Quadrant ④.
## JSON output schema (`--json`)
```jsonc
{
"schema_version": "1.0",
"tool": "rsomics-seq-grep",
"tool_version": "0.1.0",
"status": "ok",
"result": {
"input": "in.fa",
"output": "-",
"pattern_count": 3,
"total_records": 200000,
"matched": 51994
}
}
```