Conversation
A blank line after the header in an otherwise single-line file left single_line=true while seq_offset pointed at the newline, so the scan read misaligned bases and silently dropped SNPs. Detect any blank line in a record and route it through the newline-skipping scanner.
Under --include-gaps, a gap in the reference at a variable site wrote a literal '-' into the REF column, which is invalid VCF v4.2 (REF must be A/C/G/T/N) and inconsistent with the '-'→'*' mapping already used for ALT. Map REF the same way and add gap-in-REF / gap-in-ALT VCF tests.
An all-gap column ORs to BIT_GAP only, so ones==1 fell through the A/C/G/T tests without being tallied, breaking the reported variable+constant+ambiguous == length invariant. Extracted sites and fconst are unaffected; this corrects the stderr summary.
Path::parent() of a directory-less name returns Some(""), not None, so
the '.' fallback never fired and canonicalize failed with ENOENT — the
documented 'snpick -f in.fa -o snps.fasta' quickstart errored on every
run. Treat an empty parent as the current directory.
The genotype loop did one formatted write!() per field — num_var × num_samples formatted writes, the dominant cost of the VCF path on many-sample inputs. Build each row into a reused Vec<u8> with direct byte pushes (genotype index is always a single digit 0-4) and one write_all per row. Output is byte-identical; also raise the VCF writer buffer to IO_BUF.
Without it Rayon grabs every logical core, which oversubscribes shared HPC/SLURM nodes and makes wall-clock non-deterministic. -t/--threads caps the pool; rejects 0. The scan merges with a commutative OR, so thread count never affects the emitted FASTA/VCF, only speed.
The VCF contig was hardcoded to '1', forcing users to post-process every file to rename it to their reference (e.g. NC_000962.3). Make it a flag defaulting to '1' so existing output is unchanged, and reject values with whitespace that would break the tab-delimited columns.
Previously a --vcf run on an alignment with zero variable sites returned early and wrote no VCF at all, so a Snakemake/Nextflow rule declaring the .vcf as an output failed with a missing-file error. Write a valid header-only VCF (all samples listed) instead.
Route the [snpick] progress lines through a gate and move the pass-2 'Wrote N sequences' notice out of the extractor into the pipeline so it respects the flag too. Errors still print. All logs already go to stderr, so stdout stays clean for piping.
…smatch Extract scan_parallel and assert the chunk/merge path — the production path for whole-genome inputs, previously never exercised — matches a sequential scan byte-for-byte. Add tests for the ambiguous-reference REF fallback, lowercase normalization, the inconsistent-length error, and the CRLF single-line fast path.
The 'tag already exists' step wrote SKIP=true to $GITHUB_ENV, which is scoped to the release job and was never read — so merging any PR without bumping the version re-ran the whole matrix and re-published over the existing tag. Expose the check as a job output and gate the build and publish jobs on it.
The documented 'releases/latest/download/snpick' URL 404s — the release publishes per-platform assets (snpick-linux-x86_64, etc.). Point the README at the real asset names, list all four platforms and the checksum file, document the '*' gap encoding in the VCF, and correct the PR template title (was 'get_MNV').
The Project Description was copy-pasted from another project and described an MNV codon annotator. Describe what snpick actually does — extract variable sites from FASTA alignments — and its scope.
Full docs site (Home, Installation, Usage, Output formats, Benchmarks, Architecture, Contributing, Changelog) built with Material for MkDocs. Contributing and Changelog are pulled from the repo files via snippets. A Docs workflow deploys to GitHub Pages (gh-pages) on push to main.
Add navigation tabs/sections/instant + footer, light/dark/auto palette toggle with the brand colours, an announcement bar, content tabs on the install page, a Mermaid architecture diagram, code annotations, an abbreviations glossary with tooltips, tags, image lightbox, last-updated dates, HTML minification, and social cards (built in CI, where Cairo is installed). Update the Docs workflow with the native deps and full git history.
A hands-on notebook rendered via mkdocs-jupyter (execute: false, so the saved outputs ship as-is and CI needs no snpick install): build a toy alignment, extract variable sites, read the reduced FASTA/VCF, use the ASC fconst counts, plot the SNP matrix, and see gap handling.
mkdocs-jupyter does not rewrite .md links, so the Next-steps links (../usage.md etc.) resolved to /tutorials/usage.md — 404. Point them at the final page URLs (../../usage/ …).
Add a docs badge, a Documentation entry in the top nav, and a callout pointing to the GitHub Pages site (including the tutorial).
Fill in Cargo.toml (license, description, repository, keywords, MSRV 1.74). Add per-argument help text and a NOTES/EXAMPLES section to the CLI. Run the test suite on a Linux/macOS × x86_64/aarch64 matrix (the mmap offset arithmetic was previously never tested on aarch64) plus an MSRV build job.
--reference <ID> picks which record is the REF/polarity sequence (for allele polarization and the VCF ##reference line); default stays the first record. Empty sequence IDs are now rejected, and duplicate IDs error with a clear message (they yield spec-invalid VCFs and duplicate tree taxa) unless --allow-dup-ids is passed.
--stats-json <FILE> (or '-' for stdout) writes a flat JSON run summary — sequences, alignment length, variable/constant/ambiguous counts, and the fconst vector as a typed array — so pipelines stop scraping stderr for ASC parameters. --dry-run reports the statistics and exits without writing any FASTA/VCF (making -o optional).
Filter records before the scan (comma-separated IDs or @file), so site classification and fconst are computed for exactly the retained samples — a site variable only because of a dropped outlier correctly becomes constant and enters fconst, which a naive seqkit-grep + snp-sites gets wrong. Unknown IDs error; the two flags are mutually exclusive.
… --max-alleles) A sparse counting pass over the variable columns tallies per-site allele and missingness counts, then filters which variable sites reach the reduced FASTA/VCF. Filtered sites stay variable (never reclassified as constant), so fconst and ASC remain valid. --stats-json reports both the classified variable_sites and the written_sites. Counting runs only when a filter is requested.
--mask <BED> zeroes the given alignment columns before classification, so masked regions (PE/PPE, mobile elements, resistance genes) leave the output and the fconst counts entirely. --mask-ref reads the BED in reference coordinates, mapped to alignment columns via the reference's non-gap prefix positions.
--ref-coords sets VCF POS (and the contig length) to ungapped reference positions via a prefix sum over the reference's non-gap bases, instead of the raw alignment column. --sites-output writes a TSV mapping each variable site's alignment column to its reference position, REF and ALT.
--on-invalid {ignore,warn,error} guards against out-of-alphabet input (a
protein alignment or truncated download that still parses as equal-length
would otherwise silently yield zero variable sites). Default 'ignore'
skips the extra pass entirely. --check runs a read-only composition audit
(A/C/G/T, N, gap, IUPAC, invalid fractions) and exits.
-v prints thread count, layout and scan throughput; -vv adds an allele-class histogram. Level 0 stderr is unchanged. main() now exits 2 for bad input/data (InvalidData/InvalidInput) and 1 for I/O errors, so callers can distinguish user error from environment error.
pass 2 prints an in-place '[snpick] Extracting sequence N/M' line, gated on a TTY stderr, non-quiet, a real file output, and >= 500 samples, so piped/redirected runs and small inputs are unaffected.
Move the algorithm modules (fasta, scan, extract, vcf, filter, coords, audit, input, types) into a public library crate so they can be reused and documented on docs.rs; the binary now builds on top of it. Adds a crate-level doc example (run as a doctest). No behavior change.
proptest checks the classification invariant (variable+constant+ambiguous == length) and variable-position bounds over randomized alignments. A criterion micro-benchmark exercises the pass-1 scan (advisory, not a CI gate). A cargo-fuzz target fuzzes index_fasta so parser index arithmetic over untrusted bytes never panics. proptest/criterion are dev-only.
Multi-stage Dockerfile builds a static musl binary (snpick has no C dependencies) into a scratch image. Apptainer/Singularity definition bootstraps from the published image. A Container workflow builds and pushes to GHCR on version tags.
bindings/python is a maturin/pyo3 crate (its own workspace) exposing snpick.stats(path, include_gaps) — variable/constant/ambiguous counts and the fconst tuple — over the snpick library. Build with 'maturin develop'.
Gate rayon (parallel) and memmap2 (mmap) behind default-on features so the library also builds for wasm32 (no threads, no mmap) using the sequential scan. The binary declares them as required-features. bindings/wasm exposes variable_site_count() and fconst() via wasm-bindgen.
Rewrite the README usage reference and features, and update the docs site (usage, output, index, installation, architecture, CHANGELOG) to cover per-site filtering, BED masking, reference coordinates, sample selection, PHYLIP/NEXUS output, gzip/stdin-stdout I/O, --stats-json/--dry-run/--check, --on-invalid/--iupac-mode, verbosity/exit codes, and the library crate, Python/WASM bindings and containers.
These sidecar outputs were never path-checked, so pointing either at the input file truncated the in-place mmap mid-run — a SIGBUS crash for --sites-output and silent destruction of the alignment for --stats-json. They now get the same collision check as -o and the VCF, against the input and the other outputs.
The 2-byte magic peek used a single read(), which on a pipe can return fewer than 2 bytes, so gzip-over-stdin was mis-detected and parsed as FASTA. Read the magic in a fill loop. Also remove the spool temp file when decompression or mmap fails, instead of leaking one per failed run.
Render a gap reference base as 'N' (VCF v4.2 REF must be A/C/G/T/N; '*' is an ALT-only spanning-deletion allele that htslib/bcftools reject). Clamp the --sites-output ref_pos to 1 like the VCF POS, so the two agree at leading-gap reference columns. Cap MAX_SEQ_LENGTH below u32::MAX so POS / reference positions never truncate.
The sites TSV is a reporting sidecar like --stats-json, but it sat after the dry-run early return, so --dry-run --sites-output silently wrote nothing. Move it before the gate so both behave consistently.
Per-site filters recounted alleles with the strict lookup, so under --iupac-mode resolve an R/W/M was tallied as missing even though it had made the site variable — inverting filter decisions (e.g. --max-alleles 1 kept biallelic sites). Count over the same table pass 1 used, and credit every set base in a multi-flag mask (R = A|G -> both A and G).
IDs like 'iso1[London]' were written raw into the NEXUS MATRIX, so conformant readers dropped '[London]' as a comment (losing/colliding taxa) or failed to parse. Single-quote labels that contain anything other than [A-Za-z0-9_.], doubling embedded quotes.
- Reject --stats-json and --sites-output pointing at the same path (the TSV silently clobbered the JSON). - Validate --maf/--max-missing as finite fractions in [0,1] (NaN/Inf and out-of-range values were silently accepted, disabling or zeroing output). - parse_bed only skips header lines whose first token is exactly 'track'/'browser', not a chrom that starts with those letters. - --sites-output '-' now writes to stdout instead of a file named '-'.
write_vcf opened its destination with File::create unconditionally, so `--vcf-output -` created a file literally named "-" in the working directory instead of streaming the VCF to stdout — inconsistent with -o, --sites-output and --stats-json, which all treat "-" as stdout. Route the VCF through the shared open_sink helper so "-" means stdout there too. For the auto-derived case, `-o - --vcf` used to fabricate a file named "-.vcf" from the "-" stem; reject that combination with a clear message telling the user to pass --vcf-output explicitly.
…isions Two related output-routing fixes: - --check audits the alignment and exits without writing anything, yet it still validated the -o path, so `--check -o input.fa` failed with a spurious "paths resolve to same file" while `--dry-run -o input.fa` succeeded. Gate all output-path derivation and validation on a single writes_align flag so --check (like --dry-run) does not reject paths it never touches. - check_paths_differ exempts "-" because a file cannot collide with stdout, but that let several sinks target stdout at once — e.g. `--stats-json - --sites-output -` concatenated JSON and TSV onto one unparseable stream. Reject more than one output routed to "-".
parse_id_set skipped lines beginning with "#" in an @file list as comments. Sample IDs come from FASTA headers and may legitimately start with "#" (e.g. ">#iso1"), so such a line was silently dropped: the ID never entered the set, was never validated against the alignment, and the run proceeded over the wrong sample set with exit 0 — while the comma-separated form kept the same ID. Remove the comment handling so both forms behave identically; the existing "sample not found" check still catches typos.
To sniff the gzip magic, map_input read the first two bytes and then, on a hit, reopened the path to decode from the start. Reopening rewinds a regular file, but a non-seekable path — process substitution `<(gzip -c ...)` or a named FIFO — yields another handle to the same pipe whose magic bytes are already consumed, so the decoder started mid-stream and aborted with "invalid gzip header" (the identical stream over stdin decoded fine). Chain the peeked bytes back onto the same handle instead, mirroring the stdin path; this works for both seekable and non-seekable inputs.
build_mask cast BED interval bounds to u32 in the --mask-ref branch. With overflow checks off in release, an out-of-range coordinate (>= 2^32) or a degenerate empty interval at usize::MAX wrapped/truncated into the valid range and masked the wrong columns — silently dropping SNPs and corrupting the ASC fconst counts — while the alignment-coordinate branch clamps and is unaffected. Compare in usize with a saturating +1, matching that branch, so a pathological interval is a no-op instead of masking real sites.
Under --ref-coords, POS is floored to 1, but the declared ##contig length came straight from the last ungapped reference position. An all-gap reference has ungapped length 0, so the header declared length=0 while the data rows sat at POS=1 — a self-contradictory VCF (POS > contig length) that htslib rejects, written silently with exit 0. Clamp the contig length to >= 1 to match the POS clamp.
- POS is the 1-based alignment column by default (reference coordinates are opt-in via --ref-coords); the feature list wrongly implied POS was reference-anchored by default. - The VCF example header is ##reference=<ref ID>, not the literal "first_sequence"; note that ##reference carries the reference record ID. - A gap in REF renders as N (VCF v4.2 forbids -/* in REF); only ALT gaps render as *. The docs and CHANGELOG claimed a gap REF becomes *.
`-` writes to stdout for --output, --vcf-output, --sites-output and --stats-json (and reads stdin for --fasta), with at most one stdout sink per run. The usage/README notes and the --vcf-output/--sites-output table rows only listed some of these, implying the VCF and sites TSV could not stream to stdout.
The published table (250 seqs 0.9 s / 1000 seqs ~3 s, ~140 MB; snp-sites 520 MB / 3+ GB) did not match the committed results.tsv (250 seqs 1.72 s / 105 MB, snp-sites 9.38 s / 213 MB; 1000 seqs 10.27 s / 217 MB, snp-sites killed). Use the recorded numbers, note their source and hardware/thread dependence, and replace the "O(L) memory regardless of sequence count" claim (contradicted by the measured 39->217 MB curve) with accurate sub-linear wording. Fixed in benchmarks.md, index.md and README.
"No copies of the sequence data are ever made" overstated it — compressed/ piped input is spooled to a temp file and the reference sequence is copied once (both O(L)). The accurate claim is that the full N×L matrix is never materialized. Likewise, "O(L) memory" holds only for the reduced-alignment path; --vcf allocates a variable_sites × samples genotype matrix (up to O(L×N), bounded by the 4 GB guard).
- From-source install now names the minimum toolchain (Rust 1.74, matching Cargo.toml rust-version), not just "edition 2021". - README quick example printed "Mapped 63 bytes" for an input that is 90 bytes (the binary prints 90); every other line already matched. - Add release dates to the 1.0.1 (2026-03-31) and 1.0.0 (2024-11-16) changelog headers for Keep a Changelog consistency with 1.0.2.
The "Include gaps" cell said a gap is "rendered as * in the VCF" unconditionally; that holds only for ALT gaps. A gap in the reference is written N (* is not a legal REF allele).
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Bumps
mainfrom 1.0.1 → 1.0.2. This release turns snpick from a focused variable-siteextractor into a pipeline-ready tool: per-site/sample filtering, reference-anchored coordinates
and BED masking, PHYLIP/NEXUS/VCF output, gzip/stdin/stdout streaming, a machine-readable JSON
summary, a full documentation site, Python/WASM bindings and container recipes — plus an
extensive round of correctness hardening. Every addition is opt-in, so existing
snpick -f in.fa -o out.fa [--vcf]runs are byte-for-byte unchanged.62 commits · 56 files · +4,760/−202 · 58 tests green ·
clippy --all-targets --all-featuresclean.New features
Filtering & selection
--max-missing,--mac,--maf,--min-samples,--max-alleles. Filteredsites never re-enter
fconst, so ASC correction stays valid.--keep-samples/--exclude-samples(comma list or@file), applied before the scan sofconstis recomputed for the retained panel.Coordinates & references
--reference <ID>chooses the REF/polarity sequence and the VCF##reference.--ref-coordswrites VCFPOSas ungapped reference positions;--mask <BED>(+--mask-ref)excludes regions from both the output and
fconst;--sites-output <TSV>maps each site.I/O & formats
-f -) and stdout (-o -,--vcf-output -,--sites-output -,--stats-json -, at most one at a time).--format {fasta,phylip,nexus};--stats-jsontyped run summary with thefconstarray.QC & robustness
--dry-run,--check,--on-invalid {ignore,warn,error},--iupac-mode {missing,resolve},--allow-dup-ids,-v/-vv, distinct exit codes (2data,1I/O).Threading & VCF
-t/--threads,--chrom,-q/--quiet, and a valid header-only VCF when--vcfis requestedbut there are no variable sites.
Performance
table lookup, asserted by a test); VCF rows assembled in a reused buffer.
Correctness fixes
20 real bugs were found and fixed across three adversarial testing rounds — each candidate was
reproduced by running the release binary, not just by review. Highlights:
REFfor a gap reference (nowN; ALT gaps stay*).--iupac-mode resolve.fconst.-),--check, and@filesample-list edge cases.The itemized list is in
CHANGELOG.md.Documentation
Full MkDocs Material site under
docs/(home, installation, usage, output formats,architecture, benchmarks, Jupyter tutorial, changelog, tags), deployed by
.github/workflows/docs.yml. Every documented claim was re-verified against this branch's binary.Packaging & ecosystem
Split into a reusable Rust library crate + thin binary; Python (pyo3/maturin) and
WebAssembly bindings; Docker (static musl →
scratch, published to GHCR by.github/workflows/docker.yml) and Apptainer (snpick.def) recipes.Tests & CI
58 tests (unit +
tests/cli.rsend-to-end + proptest + doctest), a criterion benchmark and acargo-fuzz scaffold. CI matrix (Linux/macOS, x86_64/aarch64) with clippy
-D warnings,cargo fmt --checkand an MSRV (1.74) job.Compatibility
Backward compatible — all additions are opt-in flags; default output is unchanged. MSRV 1.74.
After merge
release.ymlbuilds the cross-platform binaries and creates tag1.0.2(novprefix).docker.ymlpublishesghcr.io/pathogenomics-lab/snpick:1.0.2and:latest.docs.ymlruns onmainand creates thegh-pagesbranch → then enableSettings → Pages → Deploy from a branch →
gh-pages/root.latestrelease match the docs.