diff --git a/.dockerignore b/.dockerignore new file mode 100644 index 0000000..d4652ef --- /dev/null +++ b/.dockerignore @@ -0,0 +1,10 @@ +target/ +fuzz/target/ +site/ +docs/ +benchmarks/ +logo/ +.git/ +.github/ +.cache/ +bindings/*/target/ diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md index 8b2ad00..198e3c0 100644 --- a/.github/CONTRIBUTING.md +++ b/.github/CONTRIBUTING.md @@ -15,10 +15,12 @@ Thank you for your interest in contributing to **snpick**! We welcome contributi - [Code of Conduct](#code-of-conduct) ## Project Description -`snpick` is a tool designed to identify **Multi-Nucleotide Variants (MNVs)** within the same codon in genomic sequences. MNVs occur when multiple Single Nucleotide Variants (SNVs) are present within the same codon, leading to the translation of a different amino acid. This tool addresses limitations in current annotation programs like **ANNOVAR** or **SnpEff**, which are primarily designed to work with individual SNVs and might overlook the actual amino acid changes resulting from MNVs. +`snpick` is a fast, memory-efficient tool for extracting variable (SNP) sites from whole-genome **FASTA alignments**. It produces reduced alignments ready for phylogenetic inference with ascertainment-bias correction (ASC) in **IQ-TREE** and **RAxML**, and can optionally emit a VCF file. Its zero-copy, memory-mapped architecture scales to thousands of genomes in seconds with minimal RAM, where matrix-in-memory tools such as **snp-sites** struggle. -### Current Limitations -**IMPORTANT**: This script currently works only with **SNVs** against a reference genome. Insertions and deletions that modify the reading frame are not supported yet. +### Scope +- Input is a **FASTA alignment**: every sequence must have the same length. The first sequence is used as the reference for `REF`/`ALT` polarization. +- Standard bases `A/C/G/T` (case-insensitive) define variability; IUPAC ambiguous bases (N, R, Y, …) are treated as missing data rather than alleles. +- Gaps (`-`) are ignored by default and can be included as a 5th character with `-g`. ## How to Contribute We appreciate all contributions, whether it’s fixing bugs, proposing new features, improving the documentation, or suggesting a new direction for the tool. diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index 6baf038..e1230f6 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -1,4 +1,4 @@ -# Pull Request Template for get_MNV +# Pull Request Template for SNPick ## Description Please include a summary of the changes made in this pull request. Mention any related issues or features that are being addressed. diff --git a/.github/workflows/docker.yml b/.github/workflows/docker.yml new file mode 100644 index 0000000..89881e5 --- /dev/null +++ b/.github/workflows/docker.yml @@ -0,0 +1,40 @@ +name: Container + +on: + push: + tags: [ "[0-9]+.[0-9]+.[0-9]+" ] + workflow_dispatch: + +permissions: + contents: read + packages: write + +jobs: + docker: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Log in to GHCR + uses: docker/login-action@v3 + with: + registry: ghcr.io + username: ${{ github.actor }} + password: ${{ secrets.GITHUB_TOKEN }} + + - name: Docker metadata + id: meta + uses: docker/metadata-action@v5 + with: + images: ghcr.io/${{ github.repository_owner }}/snpick + tags: | + type=ref,event=tag + type=raw,value=latest + + - name: Build and push + uses: docker/build-push-action@v6 + with: + context: . + push: true + tags: ${{ steps.meta.outputs.tags }} + labels: ${{ steps.meta.outputs.labels }} diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 0000000..9479fcf --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,43 @@ +name: Docs + +on: + push: + branches: [ main ] + paths: + - 'docs/**' + - 'overrides/**' + - 'mkdocs.yml' + - 'CHANGELOG.md' + - '.github/CONTRIBUTING.md' + - '.github/workflows/docs.yml' + workflow_dispatch: + +permissions: + contents: write + +jobs: + deploy: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + with: + # Full history so git-revision-date-localized can read commit dates. + fetch-depth: 0 + + - uses: actions/setup-python@v5 + with: + python-version: '3.x' + cache: pip + + - name: Install Cairo/Pango (social cards) + run: | + sudo apt-get update + sudo apt-get install -y libcairo2-dev libfreetype6-dev libffi-dev \ + libjpeg-dev libpng-dev libz-dev pngquant + + - name: Install MkDocs Material and plugins + run: pip install -r docs/requirements.txt + + - name: Build and deploy to gh-pages + # GitHub Actions sets CI=true, which enables the social-cards plugin. + run: mkdocs gh-deploy --force diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 1e96c29..d4d7046 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -20,6 +20,7 @@ jobs: outputs: version: ${{ steps.version.outputs.version }} tag: ${{ steps.version.outputs.tag }} + skip: ${{ steps.tagcheck.outputs.skip }} steps: - uses: actions/checkout@v4 @@ -32,16 +33,17 @@ jobs: echo "πŸ“¦ Version: $VERSION" - name: Check if tag already exists + id: tagcheck run: | if git ls-remote --tags origin | grep -q "refs/tags/${{ steps.version.outputs.tag }}$"; then echo "⚠️ Tag ${{ steps.version.outputs.tag }} already exists β€” skipping release." - echo "SKIP=true" >> "$GITHUB_ENV" + echo "skip=true" >> "$GITHUB_OUTPUT" fi # Build matrix for cross-platform binaries build: needs: release - if: needs.release.outputs.version != '' + if: needs.release.outputs.version != '' && needs.release.outputs.skip != 'true' strategy: fail-fast: false matrix: @@ -94,6 +96,7 @@ jobs: # Create GitHub release with all binaries publish: needs: [release, build] + if: needs.release.outputs.skip != 'true' runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index 9fc0c18..0361062 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -2,7 +2,7 @@ name: Rust on: push: - branches: [ "main", "1.0.1" ] + branches: [ "main", "1.0.2" ] pull_request: branches: [ "main" ] @@ -10,19 +10,34 @@ env: CARGO_TERM_COLOR: always jobs: - build: - runs-on: ubuntu-latest + test: + name: Test (${{ matrix.os }}) + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ ubuntu-latest, ubuntu-24.04-arm, macos-14, macos-15-intel ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v4 + + - uses: dtolnay/rust-toolchain@stable + with: + components: clippy - - name: Clippy - run: cargo clippy -- -D warnings + - name: Clippy + run: cargo clippy -- -D warnings - - name: Build - run: cargo build --verbose + - name: Test + run: cargo test --verbose - - name: Test - run: cargo test --verbose + - name: Build (release) + run: cargo build --release - - name: Build (release) - run: cargo build --release + msrv: + name: MSRV 1.74 + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: dtolnay/rust-toolchain@1.74.0 + - name: Build on the minimum supported Rust version + run: cargo build --verbose diff --git a/.gitignore b/.gitignore index e7035df..18f414b 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,20 @@ test_*.fasta test_*.vcf *.vcf !src/ + +# MkDocs build output +site/ + +# MkDocs Material social-cards / plugin cache +.cache/ + +# cargo-fuzz +fuzz/target/ +fuzz/corpus/ +fuzz/artifacts/ +fuzz/Cargo.lock + +# binding crates +bindings/*/target/ +bindings/*/Cargo.lock +bindings/wasm/pkg/ diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..7c78d3f --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,82 @@ +# Changelog + +All notable changes to SNPick are documented here. The format is based on +[Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project +adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + +## [1.0.2] - 2026-07-18 + +### Added + +**Filtering & selection** +- Per-site filtering: `--max-missing`, `--mac`, `--maf`, `--min-samples`, + `--max-alleles`. Filtered sites stay variable (never re-enter `fconst`), so ASC + correction remains valid. +- `--keep-samples` / `--exclude-samples` (comma-separated IDs or `@file`) subset + the panel before the scan, so `fconst` is recomputed for the retained samples. +- `--mask ` (and `--mask-ref`) excludes regions from both the output and `fconst`. + +**Coordinates & references** +- `--reference ` chooses the REF/polarity sequence and the VCF `##reference`. +- `--ref-coords` writes VCF `POS` as ungapped reference positions. +- `--sites-output ` maps each site's alignment column to its reference position. + +**I/O & formats** +- Transparent gzip/bgzip input, plus stdin (`-f -`) and stdout (`-o -`) streaming. +- `--format {fasta,phylip,nexus}` for the reduced alignment. +- `--stats-json ` (or `-`) writes a typed JSON run summary with the `fconst` array. + +**Robustness & QC** +- `--dry-run` (stats only) and `--check` (composition audit) exit without writing. +- `--on-invalid {ignore,warn,error}` guards against non-nucleotide input. +- `--iupac-mode resolve` optionally resolves IUPAC codes to their bases. +- Empty/duplicate sequence IDs are rejected (`--allow-dup-ids` to permit). +- `-v`/`-vv` diagnostics; distinct exit codes (`2` bad input, `1` I/O). + +**Threading & VCF** +- `-t, --threads ` pins the Rayon thread pool (deterministic wall-clock; output + is unaffected by thread count). +- `--chrom ` sets the VCF `CHROM` / `##contig` name (default `1`). +- `-q, --quiet` silences the `[snpick]` progress logs (errors still print). +- A valid header-only VCF is written when `--vcf` is requested but there are no + variable sites, so pipelines declaring the `.vcf` output no longer break. + +**Packaging** +- Published as a reusable Rust library crate, with Python (pyo3/maturin) and + WebAssembly bindings, and Docker/Apptainer container recipes. + +### Fixed +- Silent SNP loss when a single-line FASTA record contained a blank line: such + records were read with misaligned byte offsets. They now use the + newline-skipping scanner. +- Gap reference bases produced an invalid VCF `REF` of `-` under `--include-gaps`; + the `REF` is now rendered as `N` (a valid VCF base). ALT gaps stay `*` (the snp-sites + convention), so `REF` and `ALT` gap encodings differ by design. +- All-gap columns were tallied in no category under `--include-gaps`, breaking the + reported `variable + constant + ambiguous == length` invariant. +- Bare output filenames (`-o snps.fasta`, no directory) were rejected during path + resolution, breaking the documented quickstart. +- The release workflow's "tag already exists" guard was never read, so merges + without a version bump re-published over existing tags. +- The README's pre-built binary download command pointed at an asset name the + release never publishes (404). + +### Performance +- VCF data rows are assembled in a reused byte buffer instead of one formatted + write per genotype field β€” a large speedup on many-sample inputs. Output is + byte-identical. + +## [1.0.1] - 2026-03-31 + +- First Bioconda release. +- Cross-platform Build & Release CI workflow (Linux/macOS, x86_64/aarch64). +- README overhaul and benchmarks. + +## [1.0.0] - 2024-11-16 + +- Initial release: zero-copy memory-mapped extraction of variable sites from + FASTA alignments, optional VCF v4.2 output, and ASC `fconst` reporting. + +[1.0.2]: https://github.com/PathoGenOmics-Lab/snpick/compare/1.0.1...1.0.2 +[1.0.1]: https://github.com/PathoGenOmics-Lab/snpick/compare/1.0.0...1.0.1 +[1.0.0]: https://github.com/PathoGenOmics-Lab/snpick/releases/tag/1.0.0 diff --git a/Cargo.lock b/Cargo.lock index b413335..1242ac7 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2,6 +2,27 @@ # It is not intended for manual editing. version = 4 +[[package]] +name = "adler2" +version = "2.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "320119579fcad9c21884f5c4861d16174d0e06250625266f50fe6898340abefa" + +[[package]] +name = "aho-corasick" +version = "1.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ddd31a130427c27518df266943a5308ed92d4b226cc639f5a8f1002816174301" +dependencies = [ + "memchr", +] + +[[package]] +name = "anes" +version = "0.1.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4b46cbb362ab8752921c97e041f5e366ee6297bd428a31275b9fcf1e380f7299" + [[package]] name = "anstream" version = "0.6.15" @@ -51,6 +72,78 @@ dependencies = [ "windows-sys", ] +[[package]] +name = "autocfg" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" + +[[package]] +name = "bit-set" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08807e080ed7f9d5433fa9b275196cfc35414f66a0c79d864dc51a0d825231a3" +dependencies = [ + "bit-vec", +] + +[[package]] +name = "bit-vec" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5e764a1d40d510daf35e07be9eb06e75770908c27d411ee6c92109c9840eaaf7" + +[[package]] +name = "bitflags" +version = "2.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b588b76d00fde79687d7646a9b5bdf3cc0f655e0bbd080335a95d7e96f3587da" + +[[package]] +name = "bumpalo" +version = "3.20.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72f5acc6cb2ba439de613abc23857ec3d78374d8ed5ac84e9d11336e87da8649" + +[[package]] +name = "cast" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5" + +[[package]] +name = "cfg-if" +version = "1.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" + +[[package]] +name = "ciborium" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "42e69ffd6f0917f5c029256a24d0161db17cea3997d185db0d35926308770f0e" +dependencies = [ + "ciborium-io", + "ciborium-ll", + "serde", +] + +[[package]] +name = "ciborium-io" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "05afea1e0a06c9be33d539b876f1ce3692f4afea2cb41f740e7743225ed1c757" + +[[package]] +name = "ciborium-ll" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "57663b653d948a338bfb3eeba9bb2fd5fcfaecb9e199e87e1eda4d9e8b240fd9" +dependencies = [ + "ciborium-io", + "half", +] + [[package]] name = "clap" version = "4.5.60" @@ -97,6 +190,51 @@ version = "1.0.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d3fd119d74b830634cea2a0f58bbd0d54540518a14397557951e79340abc28c0" +[[package]] +name = "crc32fast" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9481c1c90cbf2ac953f07c8d4a58aa3945c425b7185c9154d67a65e4230da511" +dependencies = [ + "cfg-if", +] + +[[package]] +name = "criterion" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2b12d017a929603d80db1831cd3a24082f8137ce19c69e6447f54f5fc8d692f" +dependencies = [ + "anes", + "cast", + "ciborium", + "clap", + "criterion-plot", + "is-terminal", + "itertools", + "num-traits", + "once_cell", + "oorandom", + "plotters", + "rayon", + "regex", + "serde", + "serde_derive", + "serde_json", + "tinytemplate", + "walkdir", +] + +[[package]] +name = "criterion-plot" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6b50826342786a51a89e2da3a28f1c32b06e387201bc2d19791f622c673706b1" +dependencies = [ + "cast", + "itertools", +] + [[package]] name = "crossbeam-deque" version = "0.8.6" @@ -122,30 +260,181 @@ version = "0.8.21" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d0a5c400df2834b80a4c3327b3aad3a4c4cd4de0629063962b03235697506a28" +[[package]] +name = "crunchy" +version = "0.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "460fbee9c2c2f33933d720630a6a0bac33ba7053db5344fac858d4b8952d77d5" + [[package]] name = "either" version = "1.15.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719" +[[package]] +name = "errno" +version = "0.3.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "39cab71617ae0d63f51a36d69f866391735b51691dbda63cf6f96d042b63efeb" +dependencies = [ + "libc", + "windows-sys", +] + +[[package]] +name = "fastrand" +version = "2.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f1f227452a390804cdb637b74a86990f2a7d7ba4b7d5693aac9b4dd6defd8d6" + +[[package]] +name = "flate2" +version = "1.1.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "843fba2746e448b37e26a819579957415c8cef339bf08564fe8b7ddbd959573c" +dependencies = [ + "crc32fast", + "miniz_oxide", +] + +[[package]] +name = "fnv" +version = "1.0.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f9eec918d3f24069decb9af1554cad7c880e2da24a9afd88aca000531ab82c1" + +[[package]] +name = "futures-core" +version = "0.3.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2cd50c473c80f6d7c3670a752354b8e569b1a7cbfdc0419ec88e5edad85e0dc7" + +[[package]] +name = "futures-task" +version = "0.3.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b231ed28831efb4a61a08580c4bc233ec56bc009f4cd8f52da2c3cb97df0c109" + +[[package]] +name = "futures-util" +version = "0.3.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a77a90a256fce34da66415271e30f94ee91c57b04b8a2c042d9cf3220179deaa" +dependencies = [ + "futures-core", + "futures-task", + "pin-project-lite", + "slab", +] + +[[package]] +name = "getrandom" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "899def5c37c4fd7b2664648c28120ecec138e4d395b459e5ca34f9cce2dd77fd" +dependencies = [ + "cfg-if", + "libc", + "r-efi 5.3.0", + "wasip2", +] + +[[package]] +name = "getrandom" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "300e883d756b2e4ec94e02791f39b04b522276138852cfc41d9fb7e904106099" +dependencies = [ + "cfg-if", + "libc", + "r-efi 6.0.0", +] + +[[package]] +name = "half" +version = "2.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6ea2d84b969582b4b1864a92dc5d27cd2b77b622a8d79306834f1be5ba20d84b" +dependencies = [ + "cfg-if", + "crunchy", + "zerocopy", +] + [[package]] name = "heck" version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" +[[package]] +name = "hermit-abi" +version = "0.5.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc0fef456e4baa96da950455cd02c081ca953b141298e41db3fc7e36b1da849c" + +[[package]] +name = "is-terminal" +version = "0.4.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3640c1c38b8e4e43584d8df18be5fc6b0aa314ce6ebf51b53313d4306cca8e46" +dependencies = [ + "hermit-abi", + "libc", + "windows-sys", +] + [[package]] name = "is_terminal_polyfill" version = "1.70.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7943c866cc5cd64cbc25b2e01621d07fa8eb2a1a23160ee81ce38704e97b8ecf" +[[package]] +name = "itertools" +version = "0.10.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b0fd2260e829bddf4cb6ea802289de2f86d6a7a690192fbe91b3f46e0f2c8473" +dependencies = [ + "either", +] + +[[package]] +name = "itoa" +version = "1.0.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f42a60cbdf9a97f5d2305f08a87dc4e09308d1276d28c869c684d7777685682" + +[[package]] +name = "js-sys" +version = "0.3.103" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "53b44bfcdb3f8d5837a46dae1ca9660a837176eee74a28b229bc626816589102" +dependencies = [ + "cfg-if", + "futures-util", + "wasm-bindgen", +] + [[package]] name = "libc" version = "0.2.183" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b5b646652bf6661599e1da8901b3b9522896f01e736bad5f723fe7a3a27f899d" +[[package]] +name = "linux-raw-sys" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32a66949e030da00e8c7d4434b251670a91556f4144941d37452769c25d58a53" + +[[package]] +name = "memchr" +version = "2.8.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cf8baf1c55e62ffcace7a9f06f4bd9cd3f0c4beb022d3b367256b91b87513d98" + [[package]] name = "memmap2" version = "0.9.10" @@ -155,6 +444,80 @@ dependencies = [ "libc", ] +[[package]] +name = "miniz_oxide" +version = "0.8.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fa76a2c86f704bdb222d66965fb3d63269ce38518b83cb0575fca855ebb6316" +dependencies = [ + "adler2", + "simd-adler32", +] + +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", +] + +[[package]] +name = "once_cell" +version = "1.21.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50" + +[[package]] +name = "oorandom" +version = "11.1.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d6790f58c7ff633d8771f42965289203411a5e5c68388703c06e14f24770b41e" + +[[package]] +name = "pin-project-lite" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a89322df9ebe1c1578d689c92318e070967d1042b512afbe49518723f4e6d5cd" + +[[package]] +name = "plotters" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5aeb6f403d7a4911efb1e33402027fc44f29b5bf6def3effcc22d7bb75f2b747" +dependencies = [ + "num-traits", + "plotters-backend", + "plotters-svg", + "wasm-bindgen", + "web-sys", +] + +[[package]] +name = "plotters-backend" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "df42e13c12958a16b3f7f4386b9ab1f3e7933914ecea48da7139435263a4172a" + +[[package]] +name = "plotters-svg" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "51bae2ac328883f7acdfea3d66a7c35751187f870bc81f94563733a154d7a670" +dependencies = [ + "plotters-backend", +] + +[[package]] +name = "ppv-lite86" +version = "0.2.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85eae3c4ed2f50dcfe72643da4befc30deadb458a9b590d720cde2f2b1e97da9" +dependencies = [ + "zerocopy", +] + [[package]] name = "proc-macro2" version = "1.0.87" @@ -164,15 +527,90 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "proptest" +version = "1.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4b45fcc2344c680f5025fe57779faef368840d0bd1f42f216291f0dc4ace4744" +dependencies = [ + "bit-set", + "bit-vec", + "bitflags", + "num-traits", + "rand", + "rand_chacha", + "rand_xorshift", + "regex-syntax", + "rusty-fork", + "tempfile", + "unarray", +] + +[[package]] +name = "quick-error" +version = "1.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a1d01941d82fa2ab50be1e79e6714289dd7cde78eba4c074bc5a4374f650dfe0" + [[package]] name = "quote" -version = "1.0.37" +version = "1.0.46" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b5b9d34b8991d19d98081b46eacdd8eb58c6f2b201139f7c5f643cc155a633af" +checksum = "dfbc457d0c7a0759a614551b11a6409e5951f6c7537be1f1b7682b9ae9230368" dependencies = [ "proc-macro2", ] +[[package]] +name = "r-efi" +version = "5.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "69cdb34c158ceb288df11e18b4bd39de994f6657d83847bdffdbd7f346754b0f" + +[[package]] +name = "r-efi" +version = "6.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f8dcc9c7d52a811697d2151c701e0d08956f92b0e24136cf4cf27b57a6a0d9bf" + +[[package]] +name = "rand" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9ef1d0d795eb7d84685bca4f72f3649f064e6641543d3a8c415898726a57b41" +dependencies = [ + "rand_chacha", + "rand_core", +] + +[[package]] +name = "rand_chacha" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d3022b5f1df60f26e1ffddd6c66e8aa15de382ae63b3a0c1bfc0e4d3e3f325cb" +dependencies = [ + "ppv-lite86", + "rand_core", +] + +[[package]] +name = "rand_core" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "76afc826de14238e6e8c374ddcc1fa19e374fd8dd986b0d2af0d02377261d83c" +dependencies = [ + "getrandom 0.3.4", +] + +[[package]] +name = "rand_xorshift" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "513962919efc330f829edb2535844d1b912b0fbe2ca165d613e4e8788bb05a5a" +dependencies = [ + "rand_core", +] + [[package]] name = "rayon" version = "1.11.0" @@ -193,12 +631,139 @@ dependencies = [ "crossbeam-utils", ] +[[package]] +name = "regex" +version = "1.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f020237b6c8eed93db2e2cb53c00c60a8e1bc73da7d073199a1180401450218d" +dependencies = [ + "aho-corasick", + "memchr", + "regex-automata", + "regex-syntax", +] + +[[package]] +name = "regex-automata" +version = "0.4.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8fcfdb36bda0c880c5931cdc7a2bcdc8ba4556847b9d912bca70bc94708711ad" +dependencies = [ + "aho-corasick", + "memchr", + "regex-syntax", +] + +[[package]] +name = "regex-syntax" +version = "0.8.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d6f6ff9a378485b298a5286656da665ba74413d36db0979633275d2e708145d4" + +[[package]] +name = "rustix" +version = "1.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6fe4565b9518b83ef4f91bb47ce29620ca828bd32cb7e408f0062e9930ba190" +dependencies = [ + "bitflags", + "errno", + "libc", + "linux-raw-sys", + "windows-sys", +] + +[[package]] +name = "rustversion" +version = "1.0.23" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cf54715a573b99ac80df0bc206da022bcd442c974952c7b9720069370852e21f" + +[[package]] +name = "rusty-fork" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cc6bf79ff24e648f6da1f8d1f011e9cac26491b619e6b9280f2b47f1774e6ee2" +dependencies = [ + "fnv", + "quick-error", + "tempfile", + "wait-timeout", +] + +[[package]] +name = "same-file" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93fc1dc3aaa9bfed95e02e6eadabb4baf7e3078b0bd1b4d7b6b0b68378900502" +dependencies = [ + "winapi-util", +] + +[[package]] +name = "serde" +version = "1.0.228" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9a8e94ea7f378bd32cbbd37198a4a91436180c5bb472411e48b5ec2e2124ae9e" +dependencies = [ + "serde_core", + "serde_derive", +] + +[[package]] +name = "serde_core" +version = "1.0.228" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "41d385c7d4ca58e59fc732af25c3983b67ac852c1a25000afe1175de458b67ad" +dependencies = [ + "serde_derive", +] + +[[package]] +name = "serde_derive" +version = "1.0.228" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d540f220d3187173da220f885ab66608367b6574e925011a9353e4badda91d79" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "serde_json" +version = "1.0.150" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e8014e44b4736ed0538adeecded0fce2a272f22dc9578a7eb6b2d9993c74cfb9" +dependencies = [ + "itoa", + "memchr", + "serde", + "serde_core", + "zmij", +] + +[[package]] +name = "simd-adler32" +version = "0.3.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3a219298ac11a56ea9a6d2120044824d6f01aeb034955e7af7bc16858527deea" + +[[package]] +name = "slab" +version = "0.4.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0c790de23124f9ab44544d7ac05d60440adc586479ce501c1d6d7da3cd8c9cf5" + [[package]] name = "snpick" -version = "1.0.1" +version = "1.0.2" dependencies = [ "clap", + "criterion", + "flate2", "memmap2", + "proptest", "rayon", ] @@ -219,6 +784,35 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "tempfile" +version = "3.27.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32497e9a4c7b38532efcdebeef879707aa9f794296a4f0244f6f69e9bc8574bd" +dependencies = [ + "fastrand", + "getrandom 0.4.3", + "once_cell", + "rustix", + "windows-sys", +] + +[[package]] +name = "tinytemplate" +version = "1.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "be4d6b5f19ff7664e8c98d03e2139cb510db9b0a60b55f8e8709b689d939b6bc" +dependencies = [ + "serde", + "serde_json", +] + +[[package]] +name = "unarray" +version = "0.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eaea85b334db583fe3274d12b4cd1880032beab409c0d774be044d4480ab9a94" + [[package]] name = "unicode-ident" version = "1.0.13" @@ -231,6 +825,98 @@ version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "06abde3611657adf66d383f00b093d7faecc7fa57071cce2578660c9f1010821" +[[package]] +name = "wait-timeout" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "09ac3b126d3914f9849036f826e054cbabdc8519970b8998ddaf3b5bd3c65f11" +dependencies = [ + "libc", +] + +[[package]] +name = "walkdir" +version = "2.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "29790946404f91d9c5d06f9874efddea1dc06c5efe94541a7d6863108e3a5e4b" +dependencies = [ + "same-file", + "winapi-util", +] + +[[package]] +name = "wasip2" +version = "1.0.4+wasi-0.2.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b67efb37e106e55ce722a510d6b5f9c17f083e5fc79afc2badeb12cc313d9487" +dependencies = [ + "wit-bindgen", +] + +[[package]] +name = "wasm-bindgen" +version = "0.2.126" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4b067c0c11094aef6b7a801c1e34a26affafdf3d051dba08456b868789aaf9a4" +dependencies = [ + "cfg-if", + "once_cell", + "rustversion", + "wasm-bindgen-macro", + "wasm-bindgen-shared", +] + +[[package]] +name = "wasm-bindgen-macro" +version = "0.2.126" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "167ce5e579f6bcf889c4f7175a8a5a585de84e8ff93976ce393efa5f2837aab1" +dependencies = [ + "quote", + "wasm-bindgen-macro-support", +] + +[[package]] +name = "wasm-bindgen-macro-support" +version = "0.2.126" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f3997c7839262f4ef12cf90b818d6340c18e80f263f1a94bf157d0ec4420380e" +dependencies = [ + "bumpalo", + "proc-macro2", + "quote", + "syn", + "wasm-bindgen-shared", +] + +[[package]] +name = "wasm-bindgen-shared" +version = "0.2.126" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dc1b4cb0cc549fcf58d7dfc081778139b3d283a081644e833e84682ad71cea24" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "web-sys" +version = "0.3.103" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8622dcb61c0bcc9fffa6938bed81210af2da9a7e4a1a834b2e37a59b6dfb6141" +dependencies = [ + "js-sys", + "wasm-bindgen", +] + +[[package]] +name = "winapi-util" +version = "0.1.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c2a7b1c03c876122aa43f3020e6c3c3ee5c05081c9a00739faf7503aeba10d22" +dependencies = [ + "windows-sys", +] + [[package]] name = "windows-sys" version = "0.52.0" @@ -303,3 +989,35 @@ name = "windows_x86_64_msvc" version = "0.52.6" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "589f6da84c646204747d1270a2a5661ea66ed1cced2631d546fdfb155959f9ec" + +[[package]] +name = "wit-bindgen" +version = "0.57.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ebf944e87a7c253233ad6766e082e3cd714b5d03812acc24c318f549614536e" + +[[package]] +name = "zerocopy" +version = "0.8.54" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b7cbbc0a705a0fd05cc3676525980d2bf5a9bc4adac6d6475209a7887cf59d19" +dependencies = [ + "zerocopy-derive", +] + +[[package]] +name = "zerocopy-derive" +version = "0.8.54" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e2e817b7b52d0c7358d3246da9d69935ebb18116b2b102b4230dac079b4862f5" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "zmij" +version = "1.0.23" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "29666d0abbfad1e3dc4dcf6144730dd3a3ab225bbbdac83319345b1b44ccfc1b" diff --git a/Cargo.toml b/Cargo.toml index 257d500..3395b03 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,12 +1,48 @@ [package] name = "snpick" -version = "1.0.1" +version = "1.0.2" edition = "2021" +rust-version = "1.74" +authors = [ + "Paula Ruiz-Rodriguez ", + "Mireia Coscolla", +] +description = "Fast, memory-efficient extraction of variable sites from FASTA alignments." +license = "GPL-3.0-or-later" +repository = "https://github.com/PathoGenOmics-Lab/snpick" +homepage = "https://pathogenomics-lab.github.io/snpick/" +documentation = "https://pathogenomics-lab.github.io/snpick/" +readme = "README.md" +keywords = ["bioinformatics", "phylogenetics", "snp", "fasta", "vcf"] +categories = ["science", "command-line-utilities"] +exclude = ["benchmarks/", "docs/", "logo/", "site/", ".github/"] + +[features] +default = ["parallel", "mmap"] +# Parallel scan via Rayon (disable to build for targets without threads, e.g. wasm). +parallel = ["dep:rayon"] +# Memory-mapped / gzip / stdin input (disable for wasm, which has no mmap). +mmap = ["dep:memmap2"] [dependencies] clap = { version = "4.5.21", features = ["derive"] } -memmap2 = "0.9.10" -rayon = "1.11.0" +memmap2 = { version = "0.9.10", optional = true } +rayon = { version = "1.11.0", optional = true } +# Pure-Rust (miniz_oxide) backend keeps the build C-toolchain-free for Bioconda. +flate2 = "1.0" + +[[bin]] +name = "snpick" +path = "src/main.rs" +required-features = ["parallel", "mmap"] + +[dev-dependencies] +proptest = "1" +criterion = "0.5" + +[[bench]] +name = "scan" +harness = false [profile.release] opt-level = 3 diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..a029dd1 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,21 @@ +# syntax=docker/dockerfile:1 + +# ---- Build stage: static musl binary (snpick has no C dependencies) ---- +FROM rust:1-slim AS build +RUN apt-get update \ + && apt-get install -y --no-install-recommends musl-tools \ + && rm -rf /var/lib/apt/lists/* +RUN rustup target add x86_64-unknown-linux-musl +WORKDIR /src +COPY . . +RUN cargo build --release --target x86_64-unknown-linux-musl --bin snpick \ + && strip target/x86_64-unknown-linux-musl/release/snpick + +# ---- Runtime stage: minimal scratch image ---- +FROM scratch +LABEL org.opencontainers.image.title="snpick" \ + org.opencontainers.image.description="Fast extraction of variable sites from FASTA alignments" \ + org.opencontainers.image.source="https://github.com/PathoGenOmics-Lab/snpick" \ + org.opencontainers.image.licenses="GPL-3.0-or-later" +COPY --from=build /src/target/x86_64-unknown-linux-musl/release/snpick /usr/local/bin/snpick +ENTRYPOINT ["/usr/local/bin/snpick"] diff --git a/README.md b/README.md index 905337a..a19186e 100644 --- a/README.md +++ b/README.md @@ -10,10 +10,11 @@ [![Anaconda-Version Badge](https://anaconda.org/bioconda/snpick/badges/version.svg)](https://anaconda.org/bioconda/snpick) [![Anaconda-Downloads](https://img.shields.io/conda/dn/bioconda/snpick.svg?style=flat-square&label=downloads)](https://anaconda.org/bioconda/snpick) [![install with bioconda](https://img.shields.io/badge/install%20with-bioconda-brightgreen.svg?style=flat-square)](http://bioconda.github.io/recipes/snpick/README.html) +[![Documentation](https://img.shields.io/badge/docs-online-%23af64d1?style=flat-square)](https://pathogenomics-lab.github.io/snpick/) **Fast, memory-efficient extraction of variable sites from FASTA alignments.** -[Quick Start](#-quick-start) Β· [Features](#-features) Β· [Usage](#-usage) Β· [Benchmarks](#-benchmarks) Β· [Citation](#-citation) +**[πŸ“– Documentation](https://pathogenomics-lab.github.io/snpick/)** Β· [Quick Start](#-quick-start) Β· [Features](#-features) Β· [Usage](#-usage) Β· [Benchmarks](#-benchmarks) Β· [Citation](#-citation) Β· [Changelog](CHANGELOG.md) @@ -30,13 +31,15 @@ SNPick extracts variable (SNP) sites from whole-genome FASTA alignments. It prod **Why not snp-sites?** snp-sites works well for small datasets but struggles with large alignments β€” it loads everything into memory and scales poorly. SNPick uses a zero-copy memory-mapped architecture that handles thousands of genomes in seconds with minimal RAM. +> πŸ“– **Full documentation β€” [pathogenomics-lab.github.io/snpick](https://pathogenomics-lab.github.io/snpick/)** β€” installation, the complete usage reference, output formats, benchmarks, architecture, and a hands-on [tutorial](https://pathogenomics-lab.github.io/snpick/tutorials/tutorial/). + ### SNPick vs snp-sites | | **SNPick** | **snp-sites** | |---|---|---| | Architecture | Zero-copy mmap, parallel scan | Full matrix in memory | -| 250 seqs Γ— 4.4 Mbp | **0.9 s**, 105 MB | 9.5 s, 520 MB | -| 1000 seqs Γ— 4.4 Mbp | **~3 s**, ~140 MB | >26 min (killed), 3+ GB | +| 250 seqs Γ— 4.4 Mbp | **1.72 s**, 105 MB | 9.38 s, 213 MB | +| 1000 seqs Γ— 4.4 Mbp | **10.27 s**, 217 MB | killed (OOM) | | ASC fconst output | βœ… Built-in | ❌ Not supported | | VCF output | βœ… Optional | βœ… Default | | Gap handling | βœ… Optional (`-g`) | βœ… Default | @@ -53,11 +56,18 @@ conda install -c bioconda snpick # Extract variable sites snpick -f alignment.fasta -o snps.fasta -# With VCF output -snpick -f alignment.fasta -o snps.fasta --vcf +# With a VCF (reference-anchored POS) and a machine-readable summary +snpick -f alignment.fasta.gz -o snps.fasta --vcf --chrom NC_000962.3 \ + --ref-coords --stats-json stats.json + +# Drop singletons, keep biallelic sites, mask repetitive regions +snpick -f alignment.fasta -o snps.fasta --mac 2 --max-alleles 2 --mask exclude.bed -# Include gaps as informative -snpick -f alignment.fasta -o snps.fasta -g +# NEXUS output on 8 threads, quietly +snpick -f alignment.fasta -o snps.nex --format nexus -t 8 -q + +# Just the statistics, no files written +snpick -f alignment.fasta --dry-run --stats-json - ``` --- @@ -83,16 +93,54 @@ iqtree2 -s snps.fasta -m GTR+ASC -fconst 744123,1382922,1382180,743556 ### VCF generation -Optional VCF v4.2 output with per-sample genotypes. Reference allele taken from the first sequence. Ambiguous bases reported as missing (`.`). +Optional VCF v4.2 output with per-sample genotypes. Reference allele taken from the first sequence. Ambiguous bases reported as missing (`.`). `POS` is the 1-based alignment column (not an ungapped reference coordinate) and `CHROM` defaults to `1` β€” set it with `--chrom` (e.g. `NC_000962.3`) to match your reference. ### IUPAC and gap handling - **Ambiguous bases** (N, R, Y, etc.): not counted as alleles β€” positions are only variable if they have β‰₯2 standard bases (A, C, G, T) -- **Gaps** (`-`): ignored by default, included as a 5th character with `-g` +- **Gaps** (`-`): ignored by default, included as a 5th character with `-g`. In the VCF, gap alleles are written as `*` (an alignment-gap convention shared with snp-sites; note some downstream tools read `*` as a spanning deletion) ### Parallel processing -Automatic multi-threaded scanning via Rayon when the dataset is large enough. Falls back to single-threaded for small inputs to avoid overhead. +Automatic multi-threaded scanning via Rayon when the dataset is large enough. Falls back to single-threaded for small inputs to avoid overhead. Cap the thread count with `-t/--threads` (e.g. to match a SLURM allocation); the thread count never changes the output, only the wall-clock time. + +### Per-site filtering + +Drop low-quality or uninformative sites in the same pass β€” no round-trip through `vcftools`: + +- `--max-missing ` β€” maximum fraction of missing genotypes per site +- `--mac ` / `--maf ` β€” minimum minor-allele count / frequency (drop singletons, sequencing-error alleles) +- `--min-samples ` β€” minimum samples with data +- `--max-alleles ` β€” e.g. `2` keeps only biallelic sites + +Filtered sites are **variable** sites removed from the output only β€” they never re-enter `fconst`, so ASC stays valid. + +### Sample selection + +`--keep-samples` / `--exclude-samples` (comma-separated IDs or `@file`) subset the panel **before** the scan, so `fconst` and site classification are recomputed for exactly the retained samples β€” a site variable only because of a dropped outlier correctly becomes constant. + +### Region masking & reference coordinates + +- `--mask ` excludes regions (PE/PPE, mobile elements, resistance genes) from both the output **and** `fconst`. `--mask-ref` reads the BED in reference coordinates. +- `--reference ` chooses the REF/polarity sequence (and the VCF `##reference`). +- `--ref-coords` writes VCF `POS` as ungapped reference positions; `--sites-output ` maps each site's alignment column β†’ reference position β†’ REF/ALT. + +### Output formats + +`--format {fasta,phylip,nexus}` β€” FASTA (default), relaxed PHYLIP (IQ-TREE / RAxML) or a NEXUS DATA block (MrBayes / PAUP* / SplitsTree). + +### Compressed & streaming I/O + +Reads plain or **gzip/bgzip** FASTA transparently, from a file or **stdin** (`-f -`); writes the reduced FASTA to a file or **stdout** (`-o -`) for piping. + +### Machine-readable output & QC + +- `--stats-json ` (or `-` for stdout) β€” a typed JSON summary (counts + the `fconst` array), so pipelines never scrape stderr. +- `--dry-run` β€” report statistics without writing any output. +- `--check` β€” audit the alignment composition (A/C/G/T, N, gap, IUPAC, invalid fractions) and exit. +- `--on-invalid {ignore,warn,error}` β€” guard against non-nucleotide input (e.g. a protein alignment). +- `--iupac-mode resolve` β€” optionally resolve IUPAC codes (R = A|G, …) to their bases when classifying. +- `-v`/`-vv` β€” extra diagnostics and an allele-class histogram. --- @@ -115,28 +163,84 @@ cargo build --release # Binary at target/release/snpick ``` -### Pre-built binary (Linux) +### Pre-built binary + +Grab the binary for your platform from the [latest release](https://github.com/PathoGenOmics-Lab/snpick/releases/latest) β€” Linux (`x86_64`, `aarch64`) and macOS (`x86_64`, `aarch64`), with `SHA256SUMS.txt` published for verification: ```bash -wget https://github.com/PathoGenOmics-Lab/snpick/releases/latest/download/snpick -chmod +x snpick +# choose: snpick-linux-x86_64 | snpick-linux-aarch64 | snpick-macos-x86_64 | snpick-macos-aarch64 +curl -LO https://github.com/PathoGenOmics-Lab/snpick/releases/latest/download/snpick-linux-x86_64 +chmod +x snpick-linux-x86_64 +./snpick-linux-x86_64 --help ``` +### Container + +```bash +docker run --rm -v "$PWD:/data" ghcr.io/pathogenomics-lab/snpick \ + -f /data/alignment.fasta -o /data/snps.fasta +``` + +Also available as an [Apptainer/Singularity](snpick.def) image. + +### Library & bindings + +snpick is also a Rust **library crate** ([docs](https://pathogenomics-lab.github.io/snpick/)), with **Python** (`bindings/python`, via [maturin](https://www.maturin.rs)) and **WebAssembly** (`bindings/wasm`) bindings for use in notebooks, pipelines and the browser. + --- ## πŸ—ƒοΈ Usage ``` -snpick [OPTIONS] --fasta --output +snpick [OPTIONS] --fasta ``` -| Argument | Required | Description | -|---|---|---| -| `-f, --fasta ` | βœ… | Input FASTA alignment | -| `-o, --output ` | βœ… | Output FASTA (variable sites only) | -| `-g, --include-gaps` | | Treat gaps (`-`) as a 5th character | -| `--vcf` | | Generate VCF file (derived from output name) | -| `--vcf-output ` | | Custom VCF output path | +`-f/--fasta` is always required; `-o/--output` is required unless `--dry-run` or `--check`. +Use `-` for `--fasta` (stdin) or for any one of `--output`, `--vcf-output`, `--sites-output` or `--stats-json` (stdout; at most one at a time). See `snpick --help` for the authoritative list. + +**Core** + +| Argument | Description | +|---|---| +| `-f, --fasta ` | Input FASTA alignment (`-` = stdin; gzip/bgzip auto-detected) | +| `-o, --output ` | Output FASTA of variable sites (`-` = stdout) | +| `-g, --include-gaps` | Treat gaps (`-`) as a 5th character | +| `--format ` | `fasta` (default), `phylip`, or `nexus` | +| `-t, --threads ` | Threads for the parallel scan (default: all cores) | +| `-q, --quiet` Β· `-v`/`-vv` | Silence logs Β· increase detail | + +**VCF & coordinates** + +| Argument | Description | +|---|---| +| `--vcf` Β· `--vcf-output ` | Write a VCF (derived name, or a custom path) | +| `--chrom ` | CHROM / contig name (default `1`) | +| `--reference ` | REF/polarity sequence (default: first) | +| `--ref-coords` | VCF `POS` as ungapped reference positions | +| `--sites-output ` | Per-site alignmentβ†’reference coordinate map | + +**Filtering & masking** + +| Argument | Description | +|---|---| +| `--max-missing ` | Max fraction of missing genotypes per site | +| `--mac ` Β· `--maf ` | Min minor-allele count / frequency | +| `--min-samples ` | Min samples with data | +| `--max-alleles ` | Max distinct alleles (`2` = biallelic) | +| `--keep-samples` / `--exclude-samples ` | Subset samples (list or `@file`) | +| `--mask ` Β· `--mask-ref` | Mask regions (alignment or reference coords) | + +**QC, reporting & robustness** + +| Argument | Description | +|---|---| +| `--stats-json ` | Machine-readable JSON summary (`-` = stdout) | +| `--dry-run` Β· `--check` | Stats only Β· composition audit, then exit | +| `--on-invalid ` | `ignore` (default) Β· `warn` Β· `error` on non-nucleotides | +| `--iupac-mode ` | `missing` (default) Β· `resolve` ambiguity codes | +| `--allow-dup-ids` | Permit duplicate sequence IDs | + +**Exit codes:** `0` success Β· `1` I/O error Β· `2` bad input/data. ### Example @@ -167,7 +271,7 @@ A **stderr:** ``` -[snpick] Mapped 63 bytes. 3 sequences Γ— 18 positions. +[snpick] Mapped 90 bytes. 3 sequences Γ— 18 positions. [snpick] 1 variable, 17 constant (A:4 C:4 G:4 T:5), 0 ambiguous-only, 18 total. [snpick] ASC fconst: 4,4,4,5 [snpick] Done in 0.00s. 1 vars from 3 seqs Γ— 18 pos. @@ -191,7 +295,7 @@ Simulated *M. tuberculosis*-like genomes (4.4 Mbp, ~65% GC, 3.6% variable sites) Benchmark: length scaling

-SNPick maintains **O(L)** memory regardless of sequence count, while snp-sites requires **O(NΓ—L)**. +SNPick's memory grows far more slowly than snp-sites' **O(NΓ—L)** β€” β‰ˆ39 MB at 10 sequences up to β‰ˆ217 MB at 1000, versus snp-sites holding the full matrix in memory until it is killed on large inputs. --- diff --git a/benches/scan.rs b/benches/scan.rs new file mode 100644 index 0000000..beabc84 --- /dev/null +++ b/benches/scan.rs @@ -0,0 +1,32 @@ +//! Criterion micro-benchmark for the pass-1 bitmask scan. Advisory only β€” run locally with +//! `cargo bench`; not a CI gate (shared runners produce noisy thresholds). + +use criterion::{black_box, criterion_group, criterion_main, Criterion}; +use snpick::fasta::index_fasta; +use snpick::scan::pass1_scan; +use snpick::types::build_lookup; + +fn make_alignment(nseq: usize, len: usize) -> Vec { + let mut fa = Vec::new(); + for i in 0..nseq { + fa.extend_from_slice(format!(">s{}\n", i).as_bytes()); + let mut seq = vec![b'A'; len]; + seq[i % len] = b'C'; + seq[(i * 7) % len] = b'G'; + fa.extend_from_slice(&seq); + fa.push(b'\n'); + } + fa +} + +fn bench_scan(c: &mut Criterion) { + let fa = make_alignment(500, 20_000); + let lookup = build_lookup(false); + let (recs, sl, layout) = index_fasta(&fa).unwrap(); + c.bench_function("pass1_scan_500x20000", |b| { + b.iter(|| pass1_scan(black_box(&fa), &recs, sl, layout, &lookup)) + }); +} + +criterion_group!(benches, bench_scan); +criterion_main!(benches); diff --git a/bindings/python/Cargo.toml b/bindings/python/Cargo.toml new file mode 100644 index 0000000..2aad707 --- /dev/null +++ b/bindings/python/Cargo.toml @@ -0,0 +1,18 @@ +[package] +name = "snpick-python" +version = "1.0.2" +edition = "2021" +publish = false +description = "Python bindings for snpick" +license = "GPL-3.0-or-later" + +# Own workspace so it does not affect the parent crate's build. +[workspace] + +[lib] +name = "snpick" +crate-type = ["cdylib"] + +[dependencies] +pyo3 = { version = "0.22", features = ["extension-module", "abi3-py38"] } +snpick_core = { path = "../..", package = "snpick" } diff --git a/bindings/python/README.md b/bindings/python/README.md new file mode 100644 index 0000000..b059915 --- /dev/null +++ b/bindings/python/README.md @@ -0,0 +1,24 @@ +# snpick (Python bindings) + +Python bindings for [snpick](https://github.com/PathoGenOmics-Lab/snpick), built on the Rust +library crate with [pyo3](https://pyo3.rs) and [maturin](https://www.maturin.rs). + +## Build + +```bash +cd bindings/python +pip install maturin +maturin develop --release # or: maturin build --release +``` + +## Usage + +```python +import snpick + +s = snpick.stats("alignment.fasta") # gzip and "-" (stdin) accepted +print(s["variable_sites"], s["fconst"]) +``` + +`stats(path, include_gaps=False)` returns a dict with `sequences`, `alignment_length`, +`variable_sites`, `constant_sites`, `ambiguous_sites`, and the `fconst` tuple `(A, C, G, T)`. diff --git a/bindings/python/pyproject.toml b/bindings/python/pyproject.toml new file mode 100644 index 0000000..177d038 --- /dev/null +++ b/bindings/python/pyproject.toml @@ -0,0 +1,24 @@ +[build-system] +requires = ["maturin>=1.5,<2"] +build-backend = "maturin" + +[project] +name = "snpick" +version = "1.0.2" +description = "Fast extraction of variable sites from FASTA alignments" +readme = "README.md" +requires-python = ">=3.8" +license = { text = "GPL-3.0-or-later" } +authors = [{ name = "Paula Ruiz-Rodriguez" }] +classifiers = [ + "Programming Language :: Rust", + "Topic :: Scientific/Engineering :: Bio-Informatics", +] + +[project.urls] +Homepage = "https://github.com/PathoGenOmics-Lab/snpick" +Documentation = "https://pathogenomics-lab.github.io/snpick/" + +[tool.maturin] +module-name = "snpick" +features = ["pyo3/extension-module"] diff --git a/bindings/python/src/lib.rs b/bindings/python/src/lib.rs new file mode 100644 index 0000000..8f98af7 --- /dev/null +++ b/bindings/python/src/lib.rs @@ -0,0 +1,51 @@ +//! Python bindings for snpick, built on the `snpick` library crate. +//! +//! Build a wheel with maturin: +//! cd bindings/python && maturin develop --release +//! +//! Then in Python: +//! >>> import snpick +//! >>> snpick.stats("alignment.fasta") +//! {'sequences': 250, 'alignment_length': 4411532, 'variable_sites': 158934, ...} + +use pyo3::exceptions::PyIOError; +use pyo3::prelude::*; +use pyo3::types::PyDict; + +use snpick_core::fasta::{get_ref_seq, index_fasta}; +use snpick_core::input::map_input; +use snpick_core::scan::{analyze, pass1_scan}; +use snpick_core::types::build_lookup; + +/// Compute variable-site statistics for a FASTA alignment (gzip and '-' for stdin are accepted). +/// +/// Returns a dict with `sequences`, `alignment_length`, `variable_sites`, `constant_sites`, +/// `ambiguous_sites`, and the `fconst` tuple (A, C, G, T). +#[pyfunction] +#[pyo3(signature = (path, include_gaps = false))] +fn stats(py: Python<'_>, path: &str, include_gaps: bool) -> PyResult> { + let input = map_input(path).map_err(|e| PyIOError::new_err(e.to_string()))?; + let data = &input.mmap[..]; + let lookup = build_lookup(include_gaps); + let (records, seq_len, layout) = + index_fasta(data).map_err(|e| PyIOError::new_err(e.to_string()))?; + let bm = pass1_scan(data, &records, seq_len, layout, &lookup); + let rs = get_ref_seq(data, &records[0], seq_len, layout); + let (v, sc) = analyze(&bm, &rs, &lookup, include_gaps); + + let d = PyDict::new_bound(py); + d.set_item("sequences", records.len())?; + d.set_item("alignment_length", seq_len)?; + d.set_item("variable_sites", v.len())?; + d.set_item("constant_sites", sc.constant.total())?; + d.set_item("ambiguous_sites", sc.ambiguous)?; + d.set_item("fconst", (sc.constant.a, sc.constant.c, sc.constant.g, sc.constant.t))?; + Ok(d.unbind()) +} + +#[pymodule] +fn snpick(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add("__version__", env!("CARGO_PKG_VERSION"))?; + m.add_function(wrap_pyfunction!(stats, m)?)?; + Ok(()) +} diff --git a/bindings/wasm/Cargo.toml b/bindings/wasm/Cargo.toml new file mode 100644 index 0000000..99d4323 --- /dev/null +++ b/bindings/wasm/Cargo.toml @@ -0,0 +1,22 @@ +[package] +name = "snpick-wasm" +version = "1.0.2" +edition = "2021" +publish = false +description = "WebAssembly bindings for snpick" +license = "GPL-3.0-or-later" + +# Own workspace so it does not affect the parent crate's build. +[workspace] + +[lib] +crate-type = ["cdylib"] + +[dependencies] +wasm-bindgen = "0.2" +# The library builds without rayon/memmap for wasm (no threads, no mmap). +snpick_core = { path = "../..", package = "snpick", default-features = false } + +[profile.release] +opt-level = "s" +lto = true diff --git a/bindings/wasm/README.md b/bindings/wasm/README.md new file mode 100644 index 0000000..a44c659 --- /dev/null +++ b/bindings/wasm/README.md @@ -0,0 +1,25 @@ +# snpick (WebAssembly bindings) + +Run snpick's variable-site classification in the browser or Node via WebAssembly. + +The `snpick` library is compiled with `default-features = false`, which drops rayon (WASM has no +threads) and memmap (no mmap), using the sequential in-memory scan instead. + +## Build + +```bash +cargo install wasm-pack +cd bindings/wasm +wasm-pack build --target web +``` + +## Usage + +```js +import init, { variable_site_count, fconst } from "./pkg/snpick_wasm.js"; + +await init(); +const fasta = new TextEncoder().encode(">a\nACGT\n>b\nACGA\n"); +console.log(variable_site_count(fasta, false)); // 1 +console.log(fconst(fasta, false)); // Uint32Array [1, 1, 1, 0] +``` diff --git a/bindings/wasm/src/lib.rs b/bindings/wasm/src/lib.rs new file mode 100644 index 0000000..99b68b1 --- /dev/null +++ b/bindings/wasm/src/lib.rs @@ -0,0 +1,41 @@ +//! WebAssembly bindings for snpick. The `snpick` library is compiled with +//! `default-features = false`, so it drops rayon (no threads) and memmap (no mmap) and runs the +//! sequential in-memory scan β€” suitable for the browser / Node. +//! +//! Build with wasm-pack: +//! cd bindings/wasm && wasm-pack build --target web + +use wasm_bindgen::prelude::*; + +use snpick_core::fasta::{get_ref_seq, index_fasta}; +use snpick_core::scan::{analyze, pass1_scan}; +use snpick_core::types::build_lookup; + +/// Count the variable (SNP) sites in an in-memory FASTA alignment. Returns 0 on parse error. +#[wasm_bindgen] +pub fn variable_site_count(fasta: &[u8], include_gaps: bool) -> usize { + let lookup = build_lookup(include_gaps); + let (records, seq_len, layout) = match index_fasta(fasta) { + Ok(x) => x, + Err(_) => return 0, + }; + let bitmask = pass1_scan(fasta, &records, seq_len, layout, &lookup); + let ref_seq = get_ref_seq(fasta, &records[0], seq_len, layout); + let (variable, _) = analyze(&bitmask, &ref_seq, &lookup, include_gaps); + variable.len() +} + +/// The four ASC `fconst` constant-site counts (A, C, G, T) for an in-memory FASTA alignment. +#[wasm_bindgen] +pub fn fconst(fasta: &[u8], include_gaps: bool) -> Vec { + let lookup = build_lookup(include_gaps); + let (records, seq_len, layout) = match index_fasta(fasta) { + Ok(x) => x, + Err(_) => return vec![0, 0, 0, 0], + }; + let bitmask = pass1_scan(fasta, &records, seq_len, layout, &lookup); + let ref_seq = get_ref_seq(fasta, &records[0], seq_len, layout); + let (_, counts) = analyze(&bitmask, &ref_seq, &lookup, include_gaps); + let c = &counts.constant; + vec![c.a as u32, c.c as u32, c.g as u32, c.t as u32] +} diff --git a/docs/architecture.md b/docs/architecture.md new file mode 100644 index 0000000..5732a3a --- /dev/null +++ b/docs/architecture.md @@ -0,0 +1,72 @@ +--- +tags: + - internals +--- + +# Architecture + +```mermaid +flowchart LR + A([Input FASTA]) -->|mmap once| B[Index records] + B --> C[Pass 1: bitmask scan] + C -->|parallel Β· OR-merge| D[Analyze / classify] + D --> E{Variable columns} + B --> F[Pass 2: sparse extract] + E --> F + F --> G([Reduced FASTA]) + F --> H([VCF v4.2]) +``` + +SNPick memory-maps the input once and shares it, read-only, across two passes β€” the **full +alignment matrix is never materialized in memory**; sequence bytes are read straight from the +shared mmap. (Compressed or piped input is first spooled to a temp file so it can be mapped, and +the single reference sequence is copied out once for polarity β€” both O(L), not O(NΓ—L).) + +## Indexing + +The FASTA is scanned once to record, for each sequence, the byte offset of its data and its +length as zero-copy slices into the mmap. All sequences must have the same length (an alignment +requirement); a mismatch is a hard error. The indexer detects whether the file is single-line or +wrapped (multi-line) so later passes can pick the fastest access strategy. + +## Pass 1 β€” bitmask scan + +Each position accumulates a 5-bit mask (`A`, `C`, `G`, `T`, and gap) by OR-ing every sequence's +base at that column. Once built, each column is classified: + +- **variable** β€” more than one allele present, +- **constant** β€” exactly one standard base, +- **ambiguous-only** β€” no standard base (all N/gap without `-g`). + +For large inputs the scan runs in parallel: records are split into disjoint chunks, each thread +fills its own bitmask, and the partials are merged with OR. Because OR is commutative and +associative, the parallel result is **byte-for-byte identical** to a sequential scan β€” thread +count never affects output. Small inputs fall back to a single-threaded scan to avoid overhead. + +Before the hot loop, SNPick prefaults the mapped pages and hints the OS that access is sequential, +eliminating soft page faults on multi-gigabyte files. + +## Pass 2 β€” sparse extraction + +Only the variable columns are read back, via sparse random access into the mmap, and written to +the reduced FASTA. When a VCF is requested the same pass fills a compact genotype matrix, which +is then emitted as VCF rows assembled in a reused byte buffer (avoiding a formatted write per +genotype). + +## Design notes + +- **Lookup tables** β€” 256-byte arrays give O(1) nucleotide classification and case conversion. +- **Auto-vectorized scan** β€” the single-line hot loop uses a branchless A/C/G/T(+gap) kernel that + LLVM vectorizes (SSE2/AVX2/NEON); a test asserts it is byte-identical to the table lookup. +- **Zero-copy** β€” IDs, descriptions, and sequence bytes are slices into the single mmap. +- **O(L) memory (reduced-alignment path)** β€” for FASTA/PHYLIP/NEXUS output the working set scales + with alignment length, not the number of sequences, which is what lets SNPick handle thousands + of genomes with a small footprint. Requesting a VCF (`--vcf`) additionally holds a genotype + matrix of `variable_sites Γ— samples` bytes (so it grows with sequence count, up to O(LΓ—N)), + bounded by a 4 GB guard above which SNPick errors instead of allocating. + +The core is a **library crate** (`snpick`) of focused modules β€” `fasta` (indexing), `scan` +(pass 1 + classification), `extract` (pass 2), `vcf`, `filter` (per-site counting), `coords` +(reference coordinates & BED masking), `audit` (composition), `input` (gzip/stdin) and `types` β€” +with the CLI as a thin binary on top. `rayon` and `memmap2` sit behind default-on `parallel` and +`mmap` features, so the library also builds for `wasm32`. diff --git a/docs/assets/benchmark.png b/docs/assets/benchmark.png new file mode 100644 index 0000000..0834b60 Binary files /dev/null and b/docs/assets/benchmark.png differ diff --git a/docs/assets/benchmark_length.png b/docs/assets/benchmark_length.png new file mode 100644 index 0000000..a48e7ce Binary files /dev/null and b/docs/assets/benchmark_length.png differ diff --git a/docs/assets/logo.png b/docs/assets/logo.png new file mode 100644 index 0000000..20af2c8 Binary files /dev/null and b/docs/assets/logo.png differ diff --git a/docs/assets/logo.svg b/docs/assets/logo.svg new file mode 100644 index 0000000..c1ac613 --- /dev/null +++ b/docs/assets/logo.svg @@ -0,0 +1,75 @@ + + + +SNPcki diff --git a/docs/benchmarks.md b/docs/benchmarks.md new file mode 100644 index 0000000..f58183e --- /dev/null +++ b/docs/benchmarks.md @@ -0,0 +1,36 @@ +--- +tags: + - performance +--- + +# Benchmarks + +Benchmarks use simulated *M. tuberculosis*-like genomes (4.4 Mbp, ~65% GC, 3.6% variable sites). +The figures below are the committed run in +[`benchmarks/results.tsv`](https://github.com/PathoGenOmics-Lab/snpick/blob/main/benchmarks/results.tsv); +absolute wall-clock and peak memory depend on the host CPU and thread count. + +## Scaling by number of sequences + +

+ Benchmark: sequence scaling +

+ +## Scaling by sequence length + +

+ Benchmark: length scaling +

+ +SNPick's memory grows far more slowly than snp-sites' **O(N Γ— L)**: snp-sites holds the full +matrix in memory and is eventually killed on large inputs, while SNPick's footprint rises only +gently with sequence count (β‰ˆ39 MB at 10 sequences up to β‰ˆ217 MB at 1000). + +| Dataset | SNPick | snp-sites | +|---|---|---| +| 250 seqs Γ— 4.4 Mbp | **1.72 s**, 105 MB | 9.38 s, 213 MB | +| 1000 seqs Γ— 4.4 Mbp | **10.27 s**, 217 MB | killed (OOM) | + +!!! tip "Reproducing" + Wall-clock depends on core count; pin it with [`--threads`](usage.md#core) + for comparable runs. The extracted sites and VCF are identical regardless of thread count. diff --git a/docs/changelog.md b/docs/changelog.md new file mode 100644 index 0000000..786b75d --- /dev/null +++ b/docs/changelog.md @@ -0,0 +1 @@ +--8<-- "CHANGELOG.md" diff --git a/docs/contributing.md b/docs/contributing.md new file mode 100644 index 0000000..0c7bc9c --- /dev/null +++ b/docs/contributing.md @@ -0,0 +1 @@ +--8<-- ".github/CONTRIBUTING.md" diff --git a/docs/includes/abbreviations.md b/docs/includes/abbreviations.md new file mode 100644 index 0000000..9469645 --- /dev/null +++ b/docs/includes/abbreviations.md @@ -0,0 +1,16 @@ +*[SNP]: Single-Nucleotide Polymorphism +*[SNPs]: Single-Nucleotide Polymorphisms +*[SNV]: Single-Nucleotide Variant +*[ASC]: Ascertainment-bias Correction +*[VCF]: Variant Call Format +*[WGA]: Whole-Genome Alignment +*[IUPAC]: International Union of Pure and Applied Chemistry β€” nucleotide ambiguity codes +*[mmap]: memory-mapped file +*[HPC]: High-Performance Computing +*[SLURM]: Simple Linux Utility for Resource Management (job scheduler) +*[LTO]: Link-Time Optimization +*[CLI]: Command-Line Interface +*[REF]: Reference allele (VCF column) +*[ALT]: Alternate allele(s) (VCF column) +*[GT]: Genotype (VCF FORMAT field) +*[NS]: Number of Samples with data (VCF INFO field) diff --git a/docs/index.md b/docs/index.md new file mode 100644 index 0000000..4848eca --- /dev/null +++ b/docs/index.md @@ -0,0 +1,115 @@ +--- +tags: + - overview +--- + +# SNPick + +

+ ![SNPick logo](assets/logo.png){ width="600" } +

+ +**Fast, memory-efficient extraction of variable sites from FASTA alignments.** + +[Get started :material-rocket-launch:](installation.md){ .md-button .md-button--primary } +[Usage reference :material-console:](usage.md){ .md-button } +[View on GitHub :fontawesome-brands-github:](https://github.com/PathoGenOmics-Lab/snpick){ .md-button } + +SNPick extracts variable (SNP) sites from whole-genome FASTA alignments. It produces reduced +alignments ready for phylogenetic inference with ascertainment-bias correction (ASC) in +**IQ-TREE** and **RAxML**, and optionally generates VCF files. + +!!! question "Why not snp-sites?" + snp-sites works well for small datasets but struggles with large alignments β€” it loads the + whole matrix into memory and scales poorly. SNPick uses a zero-copy, memory-mapped + architecture that handles thousands of genomes in seconds with minimal RAM. + +## What you get + +
+ +- :material-dna:{ .lg .middle } __Variable sites only__ + + --- + + A reduced FASTA with just the informative columns β€” a drop-in, much smaller input for + phylogenetics. + +- :material-tree:{ .lg .middle } __ASC-ready__ + + --- + + Constant-site counts (`fconst`) printed for IQ-TREE's `+ASC` models, so branch lengths stay + unbiased. + +- :material-file-table:{ .lg .middle } __Optional VCF__ + + --- + + VCF v4.2 with per-sample genotypes, a configurable contig name, alignment-column `POS` + (or ungapped reference coordinates with `--ref-coords`), and gap/ambiguity handling. + +- :material-filter:{ .lg .middle } __Filter & mask__ + + --- + + Per-site missingness / MAC / MAF / allele filters, BED region masking and sample + selection β€” without breaking `fconst`. + +- :material-pipe:{ .lg .middle } __Pipeline-native__ + + --- + + gzip & stdin/stdout streaming, PHYLIP/NEXUS output, and a `--stats-json` sidecar so nothing + scrapes stderr. + +- :material-lightning-bolt:{ .lg .middle } __Built for scale__ + + --- + + Zero-copy mmap + parallel, auto-vectorized scan: **O(L)** memory, thousands of genomes in + seconds. + +
+ +## SNPick vs snp-sites + +| | **SNPick** | **snp-sites** | +|---|---|---| +| Architecture | Zero-copy mmap, parallel scan | Full matrix in memory | +| 250 seqs Γ— 4.4 Mbp | **1.72 s**, 105 MB | 9.38 s, 213 MB | +| 1000 seqs Γ— 4.4 Mbp | **10.27 s**, 217 MB | killed (OOM) | +| ASC `fconst` output | :material-check:{ .snpick-yes } Built-in | :material-close: Not supported | +| VCF output | :material-check:{ .snpick-yes } Optional | :material-check: Default | +| Gap handling | :material-check:{ .snpick-yes } Optional (`-g`) | :material-check: Default | +| IUPAC ambiguous | :material-check:{ .snpick-yes } Tracked as ambiguous | :material-alert: Treated as variant | + +## Quick start + +```bash +# Install from Bioconda +conda install -c bioconda snpick + +# Extract variable sites +snpick -f alignment.fasta -o snps.fasta + +# With a VCF, on 8 threads, quietly +snpick -f alignment.fasta -o snps.fasta --vcf -t 8 -q +``` + +## Citation + +If you use SNPick in your research, please cite: + +```bibtex +@software{snpick, + author = {Ruiz-Rodriguez, Paula and Coscolla, Mireia}, + title = {SNPick: Fast extraction of variable sites from FASTA alignments}, + url = {https://github.com/PathoGenOmics-Lab/snpick}, + doi = {10.5281/zenodo.14191809}, + license = {GPL-3.0} +} +``` + +Paula Ruiz-Rodriguez and Mireia Coscolla β€” Institute for Integrative Systems Biology, +IΒ²SysBio, University of Valencia-CSIC, Valencia, Spain. diff --git a/docs/installation.md b/docs/installation.md new file mode 100644 index 0000000..2ef9fa0 --- /dev/null +++ b/docs/installation.md @@ -0,0 +1,75 @@ +--- +tags: + - install +--- + +# Installation + +=== ":simple-anaconda: Bioconda" + + [![Bioconda version](https://anaconda.org/bioconda/snpick/badges/version.svg)](https://anaconda.org/bioconda/snpick) + + The recommended way to install SNPick: + + ```bash + conda install -c bioconda snpick + # or + mamba install -c bioconda snpick + ``` + +=== ":material-download: Pre-built binary" + + Grab the binary for your platform from the + [latest release](https://github.com/PathoGenOmics-Lab/snpick/releases/latest) β€” Linux + (`x86_64`, `aarch64`) and macOS (`x86_64`, `aarch64`), with `SHA256SUMS.txt` published for + verification: + + ```bash + # choose one: + # snpick-linux-x86_64 | snpick-linux-aarch64 | snpick-macos-x86_64 | snpick-macos-aarch64 + curl -LO https://github.com/PathoGenOmics-Lab/snpick/releases/latest/download/snpick-linux-x86_64 + chmod +x snpick-linux-x86_64 + ./snpick-linux-x86_64 --help + ``` + + !!! tip "Verify the download" + Fetch `SHA256SUMS.txt` from the same release and compare the checksum of your binary + against the matching line. + +=== ":material-language-rust: From source" + + Requires a [Rust toolchain](https://rustup.rs/) (edition 2021, Rust 1.74 or newer). + + ```bash + git clone https://github.com/PathoGenOmics-Lab/snpick.git + cd snpick + cargo build --release + # Binary at target/release/snpick + ``` + + !!! note "Release builds are slow on purpose" + The release profile enables fat LTO and a single codegen unit for maximum throughput, so + a release build takes noticeably longer than a debug build. + +=== ":material-docker: Container" + + ```bash + docker run --rm -v "$PWD:/data" ghcr.io/pathogenomics-lab/snpick \ + -f /data/alignment.fasta -o /data/snps.fasta + ``` + + An Apptainer/Singularity image is available too β€” see the + [`snpick.def`](https://github.com/PathoGenOmics-Lab/snpick/blob/main/snpick.def) recipe. + +## Library & bindings + +snpick is also a Rust **library crate** (these docs' API surface), with **Python** +([maturin](https://www.maturin.rs)/[pyo3](https://pyo3.rs), under `bindings/python`) and +**WebAssembly** (under `bindings/wasm`) bindings for notebooks, pipelines and the browser. + +## Check the install + +```bash +snpick --version +snpick --help +``` diff --git a/docs/output.md b/docs/output.md new file mode 100644 index 0000000..daee3e5 --- /dev/null +++ b/docs/output.md @@ -0,0 +1,173 @@ +--- +tags: + - vcf + - output + - phylogenetics +--- + +# Output formats + +## Reduced FASTA + +The primary output (`-o`) is a FASTA where every sequence keeps only the **variable columns**, +in alignment order. Headers (ID and description) are preserved verbatim; bases are upper-cased. + +Constant and ambiguous-only columns are dropped, so the reduced alignment is a drop-in input for +phylogenetic tools while staying orders of magnitude smaller. + +!!! info "No variable sites" + If the alignment has no variable columns, SNPick writes a valid FASTA with each record's + header and an empty sequence line, and exits `0`. + +### Alternative formats + +`--format` writes the reduced alignment as something other than FASTA: + +=== "PHYLIP" + + ```bash + snpick -f alignment.fasta -o snps.phy --format phylip + ``` + + Relaxed PHYLIP (`ntax nchar` header, then `name sequence`), read by IQ-TREE and RAxML. + +=== "NEXUS" + + ```bash + snpick -f alignment.fasta -o snps.nex --format nexus + ``` + + A NEXUS `DATA` block (`DIMENSIONS` / `FORMAT DATATYPE=DNA` / `MATRIX`) for MrBayes, PAUP* + and SplitsTree. + +## ASC `fconst` for ascertainment-bias correction + +Removing invariant sites biases branch lengths unless the model is told how many constant sites +of each base were dropped. SNPick reports these counts on stderr, ready for IQ-TREE's `+ASC` +models: + +```text +[snpick] ASC fconst: 744123,1382922,1382180,743556 +``` + +The four numbers are the constant-site counts for **A, C, G, T**. Feed them to IQ-TREE: + +```bash +iqtree2 -s snps.fasta -m GTR+ASC -fconst 744123,1382922,1382180,743556 +``` + +## VCF v4.2 + +With `--vcf` (or `--vcf-output`), SNPick also writes a VCF v4.2 file with one row per variable +site and per-sample genotypes. + +```text +##fileformat=VCFv4.2 +##source=snpick v1.0.2 +##reference=ref +##contig= +##INFO= +##FORMAT= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ref s1 s2 +1 2 . T C . PASS NS=3 GT 0 0 1 +1 4 . C T . PASS NS=3 GT 0 1 0 +``` + +| Field | Meaning | +|---|---| +| `CHROM` | Contig name β€” `1` by default, override with [`--chrom`](usage.md#vcf-coordinates) | +| `POS` | **1-based alignment column** by default, or the ungapped reference position with `--ref-coords` | +| `REF` | Base of the reference sequence (the first, or the one set by `--reference`); if that base is ambiguous, the first observed base in A, C, G, T order; if it is a gap (under `-g`), `N` | +| `ALT` | The other observed alleles, comma-separated | +| `INFO=NS` | Number of samples with data (a called base; gaps count only under `-g`) | +| `FORMAT=GT` | Per-sample allele index: `0` = REF, `1..` = the *n*-th ALT, `.` = missing/ambiguous | + +The `##reference` header carries the **ID** of the reference sequence (the first record, or the +one chosen with `--reference `), and `##contig` its length. + +!!! warning "POS is an alignment coordinate by default" + Without `--ref-coords`, `POS` is the column index in the alignment, so when the reference + contains gaps it diverges from the true genomic position and `##contig` length is the + alignment length. `--ref-coords` maps `POS` (and the contig length) onto ungapped reference + positions; `--reference ` picks which sequence is that reference. + + Under `--ref-coords`, columns where the reference has a **gap** (insertions relative to the + reference) have no reference coordinate of their own, so they take the position of the + preceding reference base. Several such variable columns therefore share one `POS` (and any + leading-gap column clamps to `POS 1`), producing duplicate-`POS` records β€” expected for + insertions, but some downstream tools may need `bcftools norm` or a sort. + +### Genotype matrix guard + +The genotype matrix is `variants Γ— samples` bytes. To avoid accidental multi-gigabyte VCFs, +SNPick refuses to emit a VCF whose matrix would exceed **4 GB** and tells you to drop `--vcf` or +reduce the input. The reduced FASTA is unaffected by this guard. + +### Header-only VCF + +If `--vcf` is requested but there are no variable sites, SNPick still writes a valid VCF +containing just the header and sample columns (no data rows), so a Snakemake/Nextflow rule that +declares the `.vcf` as an output does not break. + +## Gaps and ambiguous bases + +- **Ambiguous bases** (N, R, Y, …) are treated as **missing data**, never as alleles. A column + is variable only if it has β‰₯2 standard bases (A, C, G, T). Ambiguous genotypes are written as + `.` in the VCF. +- **Gaps** (`-`) are **ignored by default**. With `-g` a gap becomes a 5th allele. In the VCF a + gap in **ALT** is written as `*`, but a gap in **REF** (the reference sequence has a gap at a + variable column) is written as **`N`**, since VCF v4.2 forbids `-`/`*` in the REF field. The + sample whose base is that gap still gets genotype `0`, so `REF=N` there marks a gap, not an N. + +!!! note "The `*` gap encoding" + Writing ALT gaps as `*` is an alignment convention shared with snp-sites. Note that in strict + VCF v4.2, `*` denotes a spanning deletion, so some downstream tools may interpret gap sites + accordingly. A gap in REF is written as `N` instead (`*` is not a legal REF allele). + +With `--iupac-mode resolve`, IUPAC codes (R = A|G, …) are resolved to their bases when +classifying, so a column that is all-A plus one `R` becomes variable. Resolution applies only to +the presence scan; an `R` is still excluded from `NS` and genotyped as missing. + +## Machine-readable summary (`--stats-json`) + +`--stats-json ` (or `-` for stdout) writes a flat JSON run summary β€” so pipelines get the +`fconst` array as a typed field instead of scraping the `[snpick] ASC fconst:` stderr line: + +```json +{ + "snpick_version": "1.0.2", + "input": "alignment.fasta", + "reference": "H37Rv", + "sequences": 250, + "alignment_length": 4411532, + "variable_sites": 158934, + "written_sites": 141002, + "constant_sites": 4252598, + "constant_by_base": { "A": 744123, "C": 1382922, "G": 1382180, "T": 743556 }, + "ambiguous_sites": 0, + "fconst": [744123, 1382922, 1382180, 743556], + "include_gaps": false, + "threads": 8 +} +``` + +`variable_sites` is the classified count; `written_sites` is what remains after per-site +filtering. Pair with `--dry-run` to compute the summary without writing any alignment. + +## Variable-site coordinate map (`--sites-output`) + +`--sites-output ` writes one row per variable site mapping its alignment column to the +ungapped reference position and alleles β€” useful for cross-referencing gene coordinates or a +masking BED: + +```text +alignment_pos ref_pos ref alt +2 1 A T +8 7 C A +``` + +## Composition audit (`--check`) + +`--check` prints an A/C/G/T, N, gap, IUPAC and invalid-byte breakdown, then exits β€” a quick +pre-flight to catch a mis-supplied protein alignment or truncated download (see also +`--on-invalid`). diff --git a/docs/requirements.txt b/docs/requirements.txt new file mode 100644 index 0000000..853e02c --- /dev/null +++ b/docs/requirements.txt @@ -0,0 +1,5 @@ +mkdocs-material[imaging]>=9.5 +mkdocs-git-revision-date-localized-plugin>=1.2 +mkdocs-glightbox>=0.4 +mkdocs-minify-plugin>=0.8 +mkdocs-jupyter>=0.24 diff --git a/docs/stylesheets/extra.css b/docs/stylesheets/extra.css new file mode 100644 index 0000000..d14db35 --- /dev/null +++ b/docs/stylesheets/extra.css @@ -0,0 +1,73 @@ +/* ------------------------------------------------------------------ * + * SNPick β€” brand palette & polish + * ------------------------------------------------------------------ */ + +:root { + --md-primary-fg-color: #5e35b1; + --md-primary-fg-color--light: #7e57c2; + --md-primary-fg-color--dark: #4527a0; + --md-accent-fg-color: #ff0077; + --snpick-alt: #6ee36a; +} + +[data-md-color-scheme="slate"] { + --md-accent-fg-color: #ff4d97; +} + +/* Brand gradient on the primary heading of the landing page */ +.md-typeset h1#snpick { + background: linear-gradient(90deg, var(--md-primary-fg-color--light), var(--snpick-alt)); + -webkit-background-clip: text; + background-clip: text; + color: transparent; + font-weight: 800; + letter-spacing: -0.02em; +} + +/* Cards: subtle lift on hover */ +.md-typeset .grid.cards > ul > li, +.md-typeset .grid > .card { + transition: border-color 120ms ease, box-shadow 120ms ease, transform 120ms ease; +} +.md-typeset .grid.cards > ul > li:hover { + border-color: var(--md-accent-fg-color); + box-shadow: 0 4px 18px rgba(0, 0, 0, 0.12); + transform: translateY(-2px); +} + +/* Centre the hero logo and give it breathing room */ +.md-typeset h1#snpick + p img, +.md-typeset p > img[alt="SNPick logo"] { + display: block; + margin: 1.2rem auto 0.6rem; +} + +/* Tighten the comparison / option tables */ +.md-typeset table:not([class]) th { + background: var(--md-primary-fg-color); + color: #fff; +} + +/* "Last updated" footer meta a touch quieter */ +.md-source-file { + opacity: 0.85; +} + +/* Green "supported" checks in comparison tables */ +.md-typeset .snpick-yes { + color: #2e9e2e; +} +[data-md-color-scheme="slate"] .md-typeset .snpick-yes { + color: var(--snpick-alt); +} + +/* Larger brand logo in the header */ +.md-header__button.md-logo { + padding: 0.1rem 0.2rem; +} +.md-header__button.md-logo img, +.md-header__button.md-logo svg { + height: 2.1rem; + width: auto; +} + diff --git a/docs/tags.md b/docs/tags.md new file mode 100644 index 0000000..c177011 --- /dev/null +++ b/docs/tags.md @@ -0,0 +1,10 @@ +--- +hide: + - toc +--- + +# Tags + +Browse the documentation by topic. + + diff --git a/docs/tutorials/tutorial.ipynb b/docs/tutorials/tutorial.ipynb new file mode 100644 index 0000000..514c130 --- /dev/null +++ b/docs/tutorials/tutorial.ipynb @@ -0,0 +1,424 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "988158c2", + "metadata": {}, + "source": [ + "# Tutorial: from an alignment to a SNP matrix\n", + "\n", + "This hands-on tutorial walks through the full SNPick workflow on a tiny toy alignment: extracting variable sites, reading the reduced FASTA and VCF, using the ASC `fconst` counts, visualising the SNP matrix, and handling gaps.\n", + "\n", + "It assumes `snpick` is on your `PATH` (e.g. `conda install -c bioconda snpick`). Every cell below is real, executed output." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "01e857c7", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-18T16:42:27.061441Z", + "iopub.status.busy": "2026-07-18T16:42:27.061315Z", + "iopub.status.idle": "2026-07-18T16:42:27.601688Z", + "shell.execute_reply": "2026-07-18T16:42:27.601051Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "snpick 1.0.2\r\n" + ] + } + ], + "source": [ + "!snpick --version" + ] + }, + { + "cell_type": "markdown", + "id": "4cecf9ff", + "metadata": {}, + "source": [ + "## 1. Build a small example alignment\n", + "\n", + "We write six short sequences of equal length (an alignment is required). Most columns are constant; a handful carry SNPs, plus one ambiguous base (`N`)." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7e55b187", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-18T16:42:27.603110Z", + "iopub.status.busy": "2026-07-18T16:42:27.603021Z", + "iopub.status.idle": "2026-07-18T16:42:27.605921Z", + "shell.execute_reply": "2026-07-18T16:42:27.605522Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + ">ref\n", + "ACGTACGTACGTACGTACGTACGTACGTAC\n", + ">sample1\n", + "ACGTTCGTACGCACGTACGTACGAACGTAC\n", + ">sample2\n", + "ACGAACGTACGCACGTAGGTACGTACGTAC\n", + ">sample3\n", + "ACGTTCGTACGTACGTACGTACGAACGTAC\n", + ">sample4\n", + "ACGTACGTNCGCACGTAGGTACGTACGTAC\n", + ">sample5\n", + "ACGAACGTACGCACGTACGTACGTACGTAC\n", + "\n" + ] + } + ], + "source": [ + "seqs = {\n", + " \"ref\": \"ACGTACGTACGTACGTACGTACGTACGTAC\",\n", + " \"sample1\": \"ACGTTCGTACGCACGTACGTACGAACGTAC\",\n", + " \"sample2\": \"ACGAACGTACGCACGTAGGTACGTACGTAC\",\n", + " \"sample3\": \"ACGTTCGTACGTACGTACGTACGAACGTAC\",\n", + " \"sample4\": \"ACGTACGTNCGCACGTAGGTACGTACGTAC\",\n", + " \"sample5\": \"ACGAACGTACGCACGTACGTACGTACGTAC\",\n", + "}\n", + "with open(\"example.fasta\", \"w\") as fh:\n", + " for name, seq in seqs.items():\n", + " fh.write(f\">{name}\\n{seq}\\n\")\n", + "\n", + "print(open(\"example.fasta\").read())" + ] + }, + { + "cell_type": "markdown", + "id": "01b73957", + "metadata": {}, + "source": [ + "## 2. Extract the variable sites\n", + "\n", + "The reduced FASTA keeps only the columns that vary across samples β€” a much smaller, phylogenetics-ready alignment. Progress (including the ASC `fconst` counts) is printed to **stderr**." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "6ebde791", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-18T16:42:27.607322Z", + "iopub.status.busy": "2026-07-18T16:42:27.607227Z", + "iopub.status.idle": "2026-07-18T16:42:27.817427Z", + "shell.execute_reply": "2026-07-18T16:42:27.816853Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[snpick] Mapped 236 bytes. 6 sequences Γ— 30 positions.\r\n", + "[snpick] 5 variable, 25 constant (A:7 C:7 G:7 T:4), 0 ambiguous-only, 30 total.\r\n", + "[snpick] ASC fconst: 7,7,7,4\r\n", + "[snpick] Pass 1 took 0.00s.\r\n", + "[snpick] Pass 2: Wrote 6 sequences to snps.fasta.\r\n", + "[snpick] Done in 0.00s. 5 vars from 6 seqs Γ— 30 pos.\r\n" + ] + } + ], + "source": [ + "!snpick -f example.fasta -o snps.fasta 2>&1" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "2c5608da", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-18T16:42:27.819040Z", + "iopub.status.busy": "2026-07-18T16:42:27.818913Z", + "iopub.status.idle": "2026-07-18T16:42:27.821360Z", + "shell.execute_reply": "2026-07-18T16:42:27.820987Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + ">ref\n", + "TATCT\n", + ">sample1\n", + "TTCCA\n", + ">sample2\n", + "AACGT\n", + ">sample3\n", + "TTTCA\n", + ">sample4\n", + "TACGT\n", + ">sample5\n", + "AACCT\n", + "\n" + ] + } + ], + "source": [ + "print(open(\"snps.fasta\").read())" + ] + }, + { + "cell_type": "markdown", + "id": "e66b596e", + "metadata": {}, + "source": [ + "The line `ASC fconst: A,C,G,T` reports how many constant sites of each base were dropped. Feed it straight to IQ-TREE so branch lengths stay unbiased:\n", + "\n", + "```bash\n", + "iqtree2 -s snps.fasta -m GTR+ASC -fconst \n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "78a668f1", + "metadata": {}, + "source": [ + "## 3. Generate a VCF\n", + "\n", + "Add `--vcf` for a VCF v4.2 with per-sample genotypes. `--chrom` sets the contig name so it lines up with your reference." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "52b82c84", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-18T16:42:27.822650Z", + "iopub.status.busy": "2026-07-18T16:42:27.822550Z", + "iopub.status.idle": "2026-07-18T16:42:28.026198Z", + "shell.execute_reply": "2026-07-18T16:42:28.025443Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[snpick] Mapped 236 bytes. 6 sequences Γ— 30 positions.\r\n", + "[snpick] 5 variable, 25 constant (A:7 C:7 G:7 T:4), 0 ambiguous-only, 30 total.\r\n", + "[snpick] ASC fconst: 7,7,7,4\r\n", + "[snpick] Pass 1 took 0.00s.\r\n", + "[snpick] Pass 2: Wrote 6 sequences to snps.fasta.\r\n", + "[snpick] VCF written to snps.vcf.\r\n", + "[snpick] Done in 0.00s. 5 vars from 6 seqs Γ— 30 pos.\r\n" + ] + } + ], + "source": [ + "!snpick -f example.fasta -o snps.fasta --vcf --chrom NC_000962.3 2>&1" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "1d893854", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-18T16:42:28.028015Z", + "iopub.status.busy": "2026-07-18T16:42:28.027866Z", + "iopub.status.idle": "2026-07-18T16:42:28.030400Z", + "shell.execute_reply": "2026-07-18T16:42:28.030036Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tref\tsample1\tsample2\tsample3\tsample4\tsample5\n", + "NC_000962.3\t4\t.\tT\tA\t.\tPASS\tNS=6\tGT\t0\t0\t1\t0\t0\t1\n", + "NC_000962.3\t5\t.\tA\tT\t.\tPASS\tNS=6\tGT\t0\t1\t0\t1\t0\t0\n", + "NC_000962.3\t12\t.\tT\tC\t.\tPASS\tNS=6\tGT\t0\t1\t1\t0\t1\t1\n", + "NC_000962.3\t18\t.\tC\tG\t.\tPASS\tNS=6\tGT\t0\t0\t1\t0\t1\t0\n", + "NC_000962.3\t24\t.\tT\tA\t.\tPASS\tNS=6\tGT\t0\t1\t0\t1\t0\t0\n" + ] + } + ], + "source": [ + "vcf = open(\"snps.vcf\").read()\n", + "# Show the column header and the data rows (skip the ## metadata lines)\n", + "for line in vcf.splitlines():\n", + " if not line.startswith(\"##\"):\n", + " print(line)" + ] + }, + { + "cell_type": "markdown", + "id": "54439878", + "metadata": {}, + "source": [ + "Note how `sample4`'s `N` at a variable site becomes a missing genotype (`.`) and is excluded from `NS` (number of samples with data)." + ] + }, + { + "cell_type": "markdown", + "id": "378aaa91", + "metadata": {}, + "source": [ + "## 4. Visualise the SNP matrix\n", + "\n", + "The reduced FASTA is small enough to plot directly: samples on the y-axis, variable positions on the x-axis, coloured by allele." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "294d8c96", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-18T16:42:28.032260Z", + "iopub.status.busy": "2026-07-18T16:42:28.032114Z", + "iopub.status.idle": "2026-07-18T16:42:28.486735Z", + "shell.execute_reply": "2026-07-18T16:42:28.486257Z" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArEAAAE1CAYAAADwE7h+AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOoxJREFUeJzt3QmcjXX///HPGIxl7Pu+3PY9S2UpbvuapSiSNZQ1bt1u+1ZUKFKEW0T1lxSSjK3sSyYUkkKF7BSGjIbzf3w+v985v3NmYUZjzrnG6/l4nNucM9d1ne85h+73+Vyf7/cKcrlcLgEAAAAcJIW/BwAAAAAkFCEWAAAAjkOIBQAAgOMQYgEAAOA4hFgAAAA4DiEWAAAAjkOIBQAAgOMQYgEAAOA4hFgAAAA4DiEWgOnSpYsULlw4wfsFBQVJ3759JZBs2LDBxqV/Bip9r/U9v1fHc8J7AAB/ByEWCBD79u2TJ554QgoVKiRp0qSRfPnySYMGDWT69OkxwoqGk379+sU4hju4LFmyxPPY/Pnz7TH3TY9dokQJC55nzpxJktcGJKUPP/xQpk6d6u9hALjHUt7rJwBwZ9u2bZN//vOfUrBgQenRo4fkzp1bjh8/Ljt27JBp06bFGljnzJkjQ4cOlbx588brOcaNGydFihSR69evy5YtW2TmzJnyxRdfyP79+yVdunR2vFu3bt2DVwd/ePTRR+XPP/+U1KlTy/0YYvXv9QsvvODvoQC4hwixQAB4+eWXJVOmTLJr1y7JnDmzz+/Onj0bY/uyZcvKoUOH5JVXXpE333wzXs/RpEkTqVq1qv387LPPSrZs2eT111+X5cuXS/v27SVVqlTibxqib9y4YdVi/D0pUqTgfQSQrNFOAASAI0eOWDCNHmBVzpw5YzymLQWdOnWy6unJkyfv6jnr1q1rf/78889x9sRqqNRKcPny5S0Q5ciRQxo3bizh4eG3PfZLL71kISp6K0Rc/bQffPCBvf6QkBAJCwuz3/3222/SrVs3yZUrlz2uv3/33XdjHOPEiRPSqlUrSZ8+vb1XAwcOlMjIyHj3oNapU8du3rRaPWbMGGu70NedJ08eadOmjX1O3u+NnrLWcek2Os5evXrJ77//7nMsl8tl70f+/Pmt4q0V9wMHDkh8TZ48WWrUqGFfOtKmTStVqlTxaReJS1w9sW+//bYULVrUjvXggw/K5s2bY7wH7n0XL15sX7B07Poa69WrJ4cPH47x/pUrV06+++47qV27tr3GYsWKeca4ceNGeeihh+z5SpYsKevWrYsx1vh81vEdk45n5cqV8uuvv3paaO6m1xtA4KMSCwQA7YPdvn27nQLVQBAfw4cPlwULFiSoGuvNHcg0HMWle/fu1lOrVVyt3kZFRVno0TYHd1U3uhEjRsiECRNk1qxZ1hpxJ19++aUFEw2z2bNnt8ChvboPP/ywJ+RqeF61apWN5/Lly57TxHq6XEPMsWPHpH///tZasXDhQjvm3bp586Y0b95c1q9fL0899ZQMGDBArly5ImvXrrXP5x//+Idtp4FV35uuXbvac+uXgbfeekv27NkjW7du9VS2R40aZSG2adOmdtu9e7c0bNjQKs7xoV8iHnvsMXn66adtn0WLFknbtm3l888/l2bNmiXotWkLib6fjzzyiIX9X375xb4AZMmSxUJhdPp3S7+MDB48WC5duiSvvfaajWPnzp0+22lw1/dM3y8dmz6P/qxfTvSzeu6556RDhw4yadIk6/vWVpkMGTLYvvH9rOM7Jv13oY/rl5s33njDHgsNDU3Q+wTAIVwA/G7NmjWu4OBgu1WvXt3173//27V69WrXjRs3YmxbqFAhV7Nmzeznrl27utKkSeM6efKk3f/qq69c+s/6448/9mw/b948e2zdunWuc+fOuY4fP+5atGiRK1u2bK60adO6Tpw4Ydt17tzZju325Zdf2n79+/ePMYZbt255ftZt+vTpYz//61//cqVIkcI1f/78eL1u3Ve3P3DggM/j3bt3d+XJk8d1/vx5n8efeuopV6ZMmVzXrl2z+1OnTrVjLF682LPN1atXXcWKFbPH9f3wft/0NUZXu3Ztu7m9++67tu/rr78e5+vevHmzbfPBBx/4/D4sLMzn8bNnz7pSp05tn5f3ezZs2DDbLrbxROd+rW76d6JcuXKuunXr+jwe/fW5/y6434PIyEj7zKtVq+b666+/PNvpZ6Xbeb8H7n1Lly5t+7lNmzbNHt+3b5/P+6ePffjhh57HfvjhB89nu2PHDs/j+ndaH9e/kwn9rBMyJn2/vf8uA0ieaCcAAoCuQqCVWK24ffvtt1ZdatSoka1Q8Nlnn8W5n1Y9tTqq1ak7qV+/vlW5ChQoYFUyrU4tXbrUniM2n3zyiVXHRo8eHeN3+rg3zaNaRdOq4fvvvy+dO3eW+NJT0GXKlPE5lj53ixYt7Ofz5897bvqeaJVNq5lKJ6bpqX6t7rnp6eyePXvK3dLn1opwbJPp3K/7448/th5m/dy8x6en+vV9/eqrr2w7PXWu1VM9lvd7lpAJR3oa3rviqa9fK6nu9yC+tAXkwoULVh1PmfL/TsJpFVMrsbHRKrP3xDB9XnX06FGf7fQ1698pN20b0NaY0qVLWyuBm/tn9/4J+awTOiYAyR/tBECAqFatmnz66acWejTIasDU06Ea0Pbu3esT9Ny0t/GZZ56R2bNny3/+85/bHl97IbXHUwOM9h5q0NDTsrdrN9DT81mzZr3j2LWtISIiwk4j6ySxhNAVE7ydO3dO/vjjD3tNeouNe7Kb9j1q/2X0UK2v7W7p69b9vYNedD/99JMFrNj6laOPTxUvXtzn9/plIq7gGJ22DWg7gv4d8O71jf6a78Q9Fn2/vOnrjKtnVFfL8OYec/S+X21FiD4eDfn6hSn6Y977J+SzTuiYACR/hFggwGiVSQOt3jR0auVJK3+xVUTdPYDaB/rqq69af2NcdBJPXH2sf1fNmjUtZGlPaLt27eIVfGOrNCr3Ml8dO3aMs6JboUKFBI8xrtCnPbDBwcEJOpaOUQOs9nzGRkNqYtD+Y63O63JZM2bMsKqz9trOmzfPlpG61+J6X/6nE+TO291p/7v5rOM7JgDJHyEWCGDu0Hnq1Kk4t9GJRhoCdCKV96nbv0uPu3r1arl48eIdQ6lW97QFQmeG6+oFOinKPXEnoTQA6r4aLrUF4k4T4nSylQYY75Cqy49FpxU7rfrFVqHUirb369ZJQn/99Vecy47pNtoqoOE9egiPPj535db7ObQCGZ/KoZ5q1xn4+jnorH03DbEJ5R6LzuTXFRLctB1FJ3jdzReDvyshn3VCJLRKDcCZ6IkFAoD2UMZWSdKez/icHtfeWA1dGiQTy+OPP25jGjt2bIzfxTZWDUE63oMHD1qPo64ccDe00qbPrQFOA2p0GgDddLa/LjHmveTUtWvXYj01rcFTV1XwXhVAT9XrTHlv+tzak6lV5bhet1abNXiNHz8+xjYaCt1hWYOZBmFdasz7PYvv1aT0vdBAps/lpoFz2bJlcjdfiHQlCl2WTcfoptVkf52KT8hnnRC63Jq2ewBI3qjEAgFAJ/5o+GrdurWUKlXKgpZexeujjz6yfkVtKbgddzX2vffeS7QxabVO+211+S6tJGqFVU//6ilu/Z1O5IpOl0rSiydouNReXg1bd3MRBZ2opsFeK8s6EUn7gbUirJN8tAKqPyv9nYZNXTP3m2++sdPt2lqhk7ui0yXCNOzq69AQqr2vOgnNvWSWmx5Le3wHDRokX3/9tU0cunr1qj1v7969pWXLljYZTZfYmjhxorVR6JJZ+jr1fdLWD53gpq9fK426FJRup0tQ6fuiS3DpElI6eexOdAktvSCFjlmXqNL+UO1t1sq3rsua0DYVXftW/67pGsH6Hmgg1mXC9D3wV/Uyvp91QugEO/23o5+htuXoxDP9YgUgmfH38ggAXK5Vq1a5unXr5ipVqpQrNDTUlmXSZaL69evnOnPmTJxLbHn76aefbImuuJbY2rVr123HEH2JLRUVFeWaNGmSjUvHlCNHDleTJk1c33zzTaxLbLktX77clTJlSteTTz7punnzZpzPGdu+bvq69XcFChRwpUqVypU7d25XvXr1XLNnz/bZ7tdff3U99thjrnTp0rmyZ8/uGjBggGepK+8lttSUKVNc+fLlc4WEhLhq1qzpCg8Pj7HEltJlnYYPH+4qUqSI57mfeOIJ15EjR3y207FUqVLFlirLkCGDq3z58rY8mnvJM6Wvf+zYsbaMlG5Xp04d1/79++Nc8iu6uXPnuooXL25j1s9BP8/Ro0fb60vIEltub775pm2rx3vwwQddW7dutdfQuHHjGPt6/z1SP//8c4wlsvS9K1u2bIxxx/X3NLbPPD6fdULGFBER4erQoYMrc+bM9juW2wKSpyD9H38HaQCAf2h1XSvGekUybTUAAKegJxYA7hN6Od3odQttndBT9tEvvQsAgY5KLADcJzZs2GCXm9VLw+okL+07nTt3rl2UQHuKvS8iAACBjoldAHCf0EmCegECnaznXjpNJ7Lp5CoCLACnoRILAAAAx6EnFgAAAI5DO8HfmNGri6zr1Wa4OgwAAMmbnri+cuWK5M2bV1KkoAYYCAixd0kDrPaWAQCA+4de5S9//vz+HgYIsXfPfV343fWqSmhK3sZA0npoZX8PAbFYOnG3v4eAWDzZeIa/h4BYfBTW299DQDQRUVFSeX245///4X+kr7vkbiHQAJshFW9jIAkODfH3EBAL/p0EpuC0/B9yIOLfS+CihTBw0NQBAAAAxyHEAgAAwHEIsQAAAHAcQiwAAAAchxALAAAAxyHEAgAAwHEIsQAAAHAcQiwAAAAchxALAAAAxyHEAgAAwHEIsQAAAHAcLs4MAACQhCp+MyjJnuvbKq9LckUlVkRcLpf07NlTsmbNKkFBQbJ3715/DwkAAMCvtm/fLsHBwdKsWTMJRIRYEQkLC5P58+fL559/LqdOnZJy5cr5e0gAAAB+NXfuXOnXr59s2rRJTp48KYEm2bcT3LhxQ1KnTn3bbY4cOSJ58uSRGjVqJNm4AAAAAlVERIR89NFHEh4eLqdPn7Zi37BhwySQJLtKbJ06daRv377ywgsvSPbs2aVRo0ayf/9+adKkiYSGhkquXLnkmWeekfPnz9v2Xbp0sW8Zx44ds1aCwoULx3rcyMhIuXz5ss8NAAAgOVq8eLGUKlVKSpYsKR07dpR3333X2i8DSbILseq9996z6uvWrVvllVdekbp168oDDzxg3ya0deDMmTPSrl0723batGkybtw4yZ8/v7US7Nq1K9ZjTpw4UTJlyuS5FShQIIlfFQAAQNK1EnTs2NF+bty4sVy6dEk2btwogSRZhtjixYvLa6+9Zt8e1q5dawF2woQJ9o1Cf9ZvE1999ZX8+OOPFkgzZMhgjcu5c+eWHDlyxHrMoUOH2gfovh0/fjzJXxcAAMC9dujQIfn666+lffv2dj9lypTy5JNPWrANJMmyJ7ZKlSqen7/99lsLrNpKEFsvbIkSJeJ1zJCQELsBAAAkZ3PnzpWoqCjJmzev5zFtJdAc9NZbb1kBMBAkyxCbPn16n8bkFi1ayKuvvhpjO53MBQAAgP+h4XXBggUyZcoUadiwoXhr1aqV/L//9//kueeek0CQLEOst8qVK8snn3xiE7a0HA4AAIDY6XKjv//+u3Tv3j1GxfXxxx+3Ki0hNon06dNH5syZY30d//73v+2CBocPH5ZFixbJf//7X+uFBQAASCqBfBWtuXPnSv369WNtGdAQq3OOvvvuO6lQoYL4W7IPsdrPoasUDBkyxMriulRWoUKFbKZdihTJcl4bAADAXVmxYkWcv3vwwQcDapmtZBdiN2zYEOtqBZ9++mmc++iasnoDAACAM1CKBAAAgOMQYgEAAOA4hFgAAAA4DiEWAAAAjkOIBQAAgOMQYgEAAOA4hFgAAAA4DiEWAAAAjpPsLnYAAAAQyKrNOpxkz7WrVzFJrqjEAgAAwMfp06elX79+UrRoUQkJCZECBQpIixYtZP369RIoqMQi2Vk9Zpe/h4BYtGg+z99DAByj0Zhq/h4CorkZESmyeofcD3755RepWbOmZM6cWSZNmiTly5eXv/76S1avXi19+vSRH374QQIBIRYAAAAevXv3lqCgIPn6668lffr0nsfLli0r3bp1k0BBOwEAAADMxYsXJSwszCqu3gHWTauzgYIQCwAAAHP48GFxuVxSqlQpCXSEWAAAABgNsE5BiAUAAIApXry49cMGyuSt2yHEAgAAwGTNmlUaNWokb7/9tly9elWi++OPPyRQEGIBAADgoQH25s2b8uCDD8onn3wiP/30kxw8eFDefPNNqV69ugQKltgCAABIQoF+Fa2iRYvK7t275eWXX5Z//etfcurUKcmRI4dUqVJFZs6cKYGCEAsAAAAfefLkkbfeestugYp2AgAAADgOIRYAAACOQ4gFAACA4xBiAQAA4DiEWAAAADhOsgixXbp0kVatWvl7GAAAAEgiySLE3o3+/fvbemchISFSqVIlfw8HAAAACXDfhljVrVs3efLJJ/09DAAAANzrELtkyRIpX768pE2bVrJlyyb169e3a+vu2rVLGjRoINmzZ5dMmTJJ7dq17WoP3oKCgmTWrFnSvHlzSZcunZQuXVq2b98uhw8fljp16kj69OmlRo0acuTIEc8+Y8aMsUqp7legQAHbr127dnLp0qU4x3jr1i2ZOHGiFClSxMZZsWJFG7c3vXRanz597KoU8REZGSmXL1/2uQEAAMABV+zSy461b99eXnvtNWndurVcuXJFNm/eLC6Xy37u3LmzTJ8+3e5PmTJFmjZtatfbzZAhg+cY48ePl9dff91uQ4YMkQ4dOliQHDp0qBQsWNCqo3379pVVq1Z59tGQu3jxYlmxYoWFx+7du0vv3r3lgw8+iHWcGmDff/99eeedd6R48eKyadMm6dixo10yTcP13dBjjh079q72BQAAcLu8NXOSPVfGmn9IcpXgEBsVFSVt2rSRQoUK2WNalVV169b12Xb27NmSOXNm2bhxo1Ve3bp27WqVVKUhtnr16jJy5Ehp1KiRPTZgwADbxtv169dlwYIFki9fPruvQblZs2YWlHPnzh2jYjphwgRZt26dHVtpSN6yZYtVc+82xGrIHjRokOe+hmmtDAMAACQ3p0+ftgLeypUr5cSJE3aWvVixYlYU1KKlnhl3VIjV0/L16tWz4Kqhs2HDhvLEE09IlixZ5MyZMzJixAjZsGGDnD17Vm7evCnXrl2TY8eO+RyjQoUKnp9z5crlE4Tdj2lo1ZCYMWNGe0wrtO4AqzScasvAoUOHYoRYrdrq82prg7cbN27IAw88IHdLJ4DpDQAAIDk7evSo1KxZ04qRWhjUnKYZaN++fVak1Ez22GOPOSvEBgcHy9q1a2Xbtm2yZs0aq4gOHz5cdu7cKc8//7xcuHBBpk2bZlVafbEaNjU8ekuVKpVPj2xcj2lIvRsRERH2p35z8A6+ihAKAABwe9qymTJlSgkPD7f5Sm56Zrtly5bWNhoIEhRi3SFT07neRo0aZYF16dKlsnXrVpkxY4b1warjx4/L+fPnE2WQWs09efKk5M2b1+7v2LFDUqRIISVLloyxbZkyZSys6j532zoAAABwP7pw4YIVKrUC6x1gvbkLjo4KsVpxXb9+vbUR5MyZ0+6fO3fOVhnQCVQLFy6UqlWrWivAiy++aCsDJIY0adJY/8XkyZPt2LrGq/bVRm8lUDqJbPDgwTJw4ECr5taqVctWMtCQre0Jehx324FWbbXn488//5S9e/d6QnDq1KkTZdwAAABOcvjwYau0Ri8U6upT2u6pdHWnV199VRwVYjUE6kz/qVOnWpjUKqxOrmrSpIkFyp49e0rlypVtwpMmeA2TiUEbiXUymVZ5L168aBPFtOobF10BQVci0IZk7evQng4d17BhwzzbPPvsszbpzM3dL/vzzz9L4cKFE2XcAAAAycHXX39txcGnn37aJtEHggSFWK24hoWFxfo7DYG6Vqw3nfTlLXoPhYbF6I/perGx9Vpoz63eYjN//vwYZW5d5UBvcdEJaAAAAPAtHGqO0snz3tzr6ifWWfbEcF9fsQsAAAD/Ry9kpSs8vfXWW3Yxq0BGiAUAAICHtmzqdQF0ntNHH30kBw8etMqsXkjqhx9+sNWqAkGQK1DWSXAY7QnWhX9/bPSwZEiV4EUegPtOi+bz/D0EwDFuVI173gf842ZEpByoM8Mmi7vXsU/OTp06ZfOb3Bc70JWfdPJ727ZtbQkux13sAAAAAMlfnjx57HoAegtUtBMAAADAcQixAAAAcBxCLAAAAByHEAsAAADHIcQCAADAcQixAAAAcBxCLAAAAByHdWL/ptZDK0twaIi/hwEEvNTh/h4BYrO+XFV/DwGxuDamvL+HgGiu/BUlJfw9CPigEgsAAADHoRILAACQhE63eCTJniv3is2SXFGJBQAAgAQFBd32NmbMGAkkVGIBAAAgp06d8vz80UcfyahRo+TQoUOex0JDQyWQEGIBAAAguXPn9vycKVMmq756PxZoaCcAAACA4xBiAQAA4DiEWAAAADgOIRYAAACOQ4gFAACA4xBiAQAA4DgssQUAAJCEkvNVtJISlVgAAAD46NKli/zxxx8SyFIklze6VatW/h4GAAAAkkiyCLEJ9e2330r79u2lQIECkjZtWildurRMmzbN38MCAABAPN2XPbHffPON5MyZU95//30Lstu2bZOePXtKcHCw9O3b19/DAwAAQGJXYpcsWSLly5e3Cma2bNmkfv36cvXqVdm1a5c0aNBAsmfPbtfbrV27tuzevdtnX70G76xZs6R58+aSLl06q4Bu375dDh8+LHXq1JH06dNLjRo15MiRI559xowZI5UqVbL9NHDqfu3atZNLly7FOcZbt27JxIkTpUiRIjbOihUr2rjdunXrZpVXHWPRokWlY8eO0rVrV/n0008T+nYAAAAg0EPsqVOn7DS8hsCDBw/Khg0bpE2bNuJyueTKlSvSuXNn2bJli+zYsUOKFy8uTZs2tce9jR8/Xjp16iR79+6VUqVKSYcOHaRXr14ydOhQCQ8Pt2NFr4ZqyF28eLGsWLFCwsLCZM+ePdK7d+84x6kBdsGCBfLOO+/IgQMHZODAgRZUN27cGOc+GoqzZs0a5+8jIyPl8uXLPjcAAAA4oJ1AQ2xUVJQF10KFCtljWpVVdevW9dl29uzZkjlzZguOWnl104qnVlLVkCFDpHr16jJy5Ehp1KiRPTZgwADbxtv169ctlObLl8/uT58+XZo1ayZTpkyR3LlzxwibEyZMkHXr1tmxlVZbNVxrNVerr9FpO8FHH30kK1euvG0wHjt2bELeLgAAAARCJVZPy9erV8+Ca9u2bWXOnDny+++/2+/OnDkjPXr0sAqsthNkzJhRIiIi5NixYz7HqFChgufnXLly+QRh92MaWr0rnQULFvQEWKXhVFsGDh06FGOMWrW9du2atTaEhoZ6bhqCvdsU3Pbv3y8tW7aU0aNHS8OGDeN87Vop1mqt+3b8+PEEvHMAAADwWyVWJz6tXbvWKpdr1qyxiujw4cNl586d8vzzz8uFCxes11SrtCEhIRY2b9y44XOMVKlS+fTIxvWYhtS7ocFZaVXVO/gqHZO377//3kK5TuoaMWLEbY+r+0bfHwAAAA5ZnUBDZs2aNe02atQoC6xLly6VrVu3yowZM6wPVmml8vz584kySK3mnjx5UvLmzWv3tec2RYoUUrJkyRjblilTxsKm7hNb64Cb9spqC4T28b788suJMk4AAAAEYIjViuv69evttLsuUaX3z507Z6sMaBvBwoULpWrVqtYK8OKLL9rKAIkhTZo0FjYnT55sx+7fv7/11Ubvh1UZMmSQwYMH22QurebWqlXLTv9ryNYWBz2OthBogNU+3EGDBsnp06c9leYcOXIkypgBAABi8+WWHUn2XHVrPSzJVYJ6YjUEbtq0yaqtJUqUsFPwOrmqSZMmMnfuXOuPrVy5sjzzzDMWNDXoJoZixYrZZDJ9Xg3Q2lerVd+46AoIOllMJ2NpwG7cuLG1F+iSW0qX29LwrevE5smTx3OrVq1aoowXAADAqVdBDQoKkldeecXn8WXLlnlaPu9E85ZOsL/Xgly6plUA03Vi9Y3TJbkCiVaEdQJb2Q29JTiUXlngTlKHx70sHvxnfbmq/h4CYnHtlf+b8IzAcOWvKCmxeoed3dWiXnKtxHbp0sVWbNKz4EePHpUsWbLY45rFWrdubUuh3s53330njz76qBULvec83Qv35WVnAQAAEDu9kJW2bOoZ7YRavny5nQG/1wFWEWIBAADgoXOEdM19XYXqxIkTkhCfffaZLV2aFAI+xGo7QaC1EgAAACRnrVu3lkqVKtk6+vH122+/WTuBzpVKCgEfYgEAAJD0Xn31VXnvvffk4MGD8a7C6qpQesXW2Dz33HM+F6L6uwixAAAAiEEnaOlypHrV0viG2MceeyzO348bN87OrrtvSX6xAwAAANwfXnnlFWsriO0CU9GvmPrVV1/JzJkz49xGl15NrOVXFZVYAAAAxKp8+fLy9NNPy5tvvim3ExYWZtcQKFy4sCQVKrEAAABJyGlX0Ro3bpytHXunpbVu10pwLxBiAQAAYObPny/RaXU1MjJS4hIVFSVffPGFrFq1SpIS7QQAAAC4axcvXpSBAwdKtWrVJClRiQUAAMBd08laI0aMkKRGiP2blk7cLRlS8TYCd9bV3wNALOpJuL+HgFjcGDPD30NANDcjIkVW7/D3MOCFdgIAAAA4DiEWAAAAjkOIBQAAgOMQYgEAAOA4hFgAAAA4DiEWAAAAjkOIBQAAgOMQYgEAAOA4hFgAAAA4DiEWAAAAjkOIBQAAgOMQYgEAAOA4hFgAAAA4DiEWAAAAjpMsQmyXLl2kVatW/h4GAAAAkkiyCLEJdeHCBWncuLHkzZtXQkJCpECBAtK3b1+5fPmyv4cGAACAeLgvQ2yKFCmkZcuW8tlnn8mPP/4o8+fPl3Xr1slzzz3n76EBAADgXoTYJUuWSPny5SVt2rSSLVs2qV+/vly9elV27dolDRo0kOzZs0umTJmkdu3asnv3bp99g4KCZNasWdK8eXNJly6dlC5dWrZv3y6HDx+WOnXqSPr06aVGjRpy5MgRzz5jxoyRSpUq2X5aMdX92rVrJ5cuXYpzjLdu3ZKJEydKkSJFbJwVK1a0cbtlyZJFnn/+ealataoUKlRI6tWrJ71795bNmzfHeczIyEir1HrfAAAA4IAQe+rUKWnfvr1069ZNDh48KBs2bJA2bdqIy+WSK1euSOfOnWXLli2yY8cOKV68uDRt2tQe9zZ+/Hjp1KmT7N27V0qVKiUdOnSQXr16ydChQyU8PNyOpaf2vWnIXbx4saxYsULCwsJkz549FjrjogF2wYIF8s4778iBAwdk4MCB0rFjR9m4cWOs2588eVI+/fRTC963O6aGc/dNAzUAAAD8I8ilqTGetLJapUoV+eWXX6yCeTtaDc2cObN8+OGHVnm1JwsKkhEjRliQVRp2q1evLnPnzrVgrBYtWiRdu3aVP//801OJfemll+TXX3+VfPny2WMaZJs1aya//fab5M6d2yZ2/fHHH7Js2TKrmGbNmtXaA/TYbs8++6xcu3bNxuOmgXz58uX2XC1atLCgnCZNmlhfjx5Xb25aidUg+2OjhyVDqpTxfQsBIKC0aD7P30NALG5UneHvISCamxGRcqDODDsTnDFjRn8PBwmtxOppeT31ru0Ebdu2lTlz5sjvv/9uvztz5oz06NHDKrBaqdQPOCIiQo4dO+ZzjAoVKnh+zpUrl/2px/N+7Pr16z6n6wsWLOgJsErDqYbkQ4cOxRijVm01rGprQ2hoqOemlVnvNgX1xhtvWDDXIKu/GzRoUJyvXSeA6WvyvgEAAMA/ElRCDA4OlrVr18q2bdtkzZo1Mn36dBk+fLjs3LnTekx11v+0adOsSquhT8PmjRs3fI6RKlUqz89amY3rMQ2pd0ODs1q5cqVP8FU6Jm9axdWbtjVo9faRRx6RkSNHSp48ee7quQEAAJA0EnweXENmzZo17TZq1CgLrEuXLpWtW7fKjBkzrA9WHT9+XM6fP58og9Rqrvat6pJY7jYEXWGgZMmSMbYtU6aMhVXd53Y9rtG5Q7N3ywAAAACSQYjViuv69eulYcOGkjNnTrt/7tw5W2VA2wgWLlxoM/61FeDFF1+0lQESg/ap6qSxyZMn27H79+9vKxRoFTW6DBkyyODBg20ylwbTWrVqWf+KhmxtAdDjfPHFF9b+UK1aNWs10MlfOl4N5oULF06UMQMAACBAQqyGwE2bNsnUqVMtTGoVdsqUKdKkSRMLlD179pTKlSvbhKcJEyZYmEwMxYoVs1UQtMp78eJFmyimVd+46MSxHDly2IoCR48etQlmOq5hw4bZ7zVcaz+vBl2tvOp49fj/+c9/EmW8AAAACKDVCfxBVyfQVQd0Sa5AoiFeJ7CxOgEAJ2N1gsDE6gSBh9UJAs99ecUuAAAAOBshFgAAAI4T8CFW2wkCrZUAAAAA/hXwIRYAAACIjhALAAAAxyHEAgAAwHEIsQAAAHAcQiwAAAAchxALAAAAxyHEAgAAwHEIsQAAAHCclP4eAJDYGo2p5u8hIBapw3v7ewiIxfpyVf09BMTi2pjy/h4CornyV5SU8Pcg4INKLAAAAByHEAsAAADHIcQCAADAcQixAAAAcBxCLAAAAByHEAsAAADHIcQCAADAcQixAAAAcBxCLAAAAByHEAsAAADHIcQCAADAcQixAAAAcBxCLAAAABwnWYTYLl26SKtWrfw9DAAAACSRZBFi/44LFy5I/vz5JSgoSP744w9/DwcAAADxcN+H2O7du0uFChX8PQwAAADcyxC7ZMkSKV++vKRNm1ayZcsm9evXl6tXr8quXbukQYMGkj17dsmUKZPUrl1bdu/e7bOvVjtnzZolzZs3l3Tp0knp0qVl+/btcvjwYalTp46kT59eatSoIUeOHPHsM2bMGKlUqZLtV6BAAduvXbt2cunSpTjHeOvWLZk4caIUKVLExlmxYkUbd3QzZ8606uvgwYMT+jYAAADAKSH21KlT0r59e+nWrZscPHhQNmzYIG3atBGXyyVXrlyRzp07y5YtW2THjh1SvHhxadq0qT3ubfz48dKpUyfZu3evlCpVSjp06CC9evWSoUOHSnh4uB2rb9++PvtoyF28eLGsWLFCwsLCZM+ePdK7d+84x6kBdsGCBfLOO+/IgQMHZODAgdKxY0fZuHGjZ5vvv/9exo0bZ9ulSHHntyEyMlIuX77scwMAAIB/pExoiI2KirLgWqhQIXtMq7Kqbt26PtvOnj1bMmfObMFRK69uXbt2tUqqGjJkiFSvXl1GjhwpjRo1sscGDBhg23i7fv26hc18+fLZ/enTp0uzZs1kypQpkjt37hhhc8KECbJu3To7tipatKiFa63maoVYt9EwPmnSJClYsKAcPXr0jq9dg/HYsWMT8nYBAAAgECqxelq+Xr16Flzbtm0rc+bMkd9//91+d+bMGenRo4dVYLWdIGPGjBIRESHHjh3zOYZ3/2muXLl8grD7MQ2t3pVODZruAKs0nGrLwKFDh2KMUau2165ds9aG0NBQz01DsLtNQau+2sqg1dn40n20hcF9O378eLz3BQAAgB8rscHBwbJ27VrZtm2brFmzxiqiw4cPl507d8rzzz9vM/2nTZtmVdqQkBALmzdu3PA5RqpUqXx6ZON6TEPq3dDgrFauXOkTfJWOSX355Zeyb98+T5+stjAo7efV1xNbxVX3de8PAAAAB4VYd8isWbOm3UaNGmWBdenSpbJ161aZMWOG9cEqrVSeP38+UQap1dyTJ09K3rx57b723Gofa8mSJWNsW6ZMGQubuo+2DsTmk08+kT///NNzXyelaZ/v5s2b5R//+EeijBkAAAABEmK14rp+/Xpp2LCh5MyZ0+6fO3fOTs1rG8HChQulatWq1grw4osv2soAiSFNmjQ2aWzy5Ml27P79+1tfbfR+WJUhQwZbbUAnc2k1t1atWnb6X0O2tjjocaIHVXfY1tehfbwAAABIRiFWQ+CmTZtk6tSpFia1CquTq5o0aWKBsmfPnlK5cmVbCksnVyXW0lXFihWzyWRa5b148aJNFNOqb1x0BYQcOXLYZCydtKXBVMc1bNiwRBkPAAAA/CvI5W4IDVC6TuyyZctsSa5AoiFeJ7D92OhhyZAqwV0ZuIcajanm7yEgFqnD414WD/6zvlxVfw8Bsbj2yv9NeEZguPJXlJRYvcPO7mpRD/5331+xCwAAAM5DiAUAAIDjBHyI1XaCQGslAAAAgH8FfIgFAAAAoiPEAgAAwHEIsQAAAHAcQiwAAAAchxALAAAAxyHEAgAAwHEIsQAAAHAcQiwAAAAcJ6W/B+B0rYdWluDQEH8PAwDuSr394f4eAmKxQrr6ewhAwKMSCwAAAMchxAIAAMBxCLEAAABwHEIsAAAAHIcQCwAAAMchxAIAAMBxCLEAAABwHEIsAAAAHIcQCwAAAMchxAIAAMBxCLEAAABwHEIsAAAAHIcQCwAAAMchxAIAAMBxkkWI7dKli7Rq1crfwwAAAEASSRYh9m4EBQXFuC1atMjfwwIAAEA8pJT72Lx586Rx48ae+5kzZ/breAAAAHCPKrFLliyR8uXLS9q0aSVbtmxSv359uXr1quzatUsaNGgg2bNnl0yZMknt2rVl9+7dPvtqtXPWrFnSvHlzSZcunZQuXVq2b98uhw8fljp16kj69OmlRo0acuTIEc8+Y8aMkUqVKtl+BQoUsP3atWsnly5dinOMt27dkokTJ0qRIkVsnBUrVrRxR6ehNXfu3J5bmjRp4jxmZGSkXL582ecGAAAAB4TYU6dOSfv27aVbt25y8OBB2bBhg7Rp00ZcLpdcuXJFOnfuLFu2bJEdO3ZI8eLFpWnTpva4t/Hjx0unTp1k7969UqpUKenQoYP06tVLhg4dKuHh4Xasvn37+uyjIXfx4sWyYsUKCQsLkz179kjv3r3jHKcG2AULFsg777wjBw4ckIEDB0rHjh1l48aNPtv16dPHQveDDz4o7777rj337Y6p4dx900ANAAAAB7QTaIiNioqy4FqoUCF7TKuyqm7duj7bzp492yqdGhy18urWtWtXq6SqIUOGSPXq1WXkyJHSqFEje2zAgAG2jbfr169bKM2XL5/dnz59ujRr1kymTJliFdToFdMJEybIunXr7NiqaNGiFq61mqsVYjVu3Dgbs1Z216xZY6E4IiJC+vfvH+tr15A9aNAgz32txBJkAQAAHBBi9bR8vXr1LLhq6GzYsKE88cQTkiVLFjlz5oyMGDHCqrNnz56VmzdvyrVr1+TYsWM+x6hQoYLn51y5cvkEYfdjGlo1JGbMmNEeK1iwoCfAKg2n2jJw6NChGCFWq7b6vNra4O3GjRvywAMPeO5rcHbTx7UlYtKkSXGG2JCQELsBAADAYSE2ODhY1q5dK9u2bbPqpVZEhw8fLjt37pTnn39eLly4INOmTbMqrQY+DZsaHr2lSpXKp0c2rsc0pN4NraaqlStX+gRfdbsQ+tBDD1mrg1ZyCasAAADJbHUCDZk1a9a026hRoyywLl26VLZu3SozZsywPlh1/PhxOX/+fKIMUqu5J0+elLx589p97blNkSKFlCxZMsa2ZcqUsRCq+7hbB+JDe3S1okyABQAASGYhViuu69evtzaCnDlz2v1z587ZKgM6kWvhwoVStWpVawV48cUXbWWAxKCrBuikscmTJ9ux9ZS/9tVGbyVQGTJkkMGDB9tkLq3m1qpVy1Yy0JCt7Ql6HJ0gpu0PDz/8sB1bq8vaR6v7AQAAIJmFWA2BmzZtkqlTp1qY1CqsTq5q0qSJBcqePXtK5cqVbcJTYobCYsWK2WQyrfJevHjRJopp1Tcu2haQI0cOW1Hg6NGjNsFMxzVs2DBP+8Lbb79tQVdXJNDjv/7669KjR49EGS8AAADurSDX7daVCgC6TuyyZcvsdH8g0RCvS22V3dBbgkNpQQDuJHV43MviAfC14nPfVXrgf1f+ipISq3fY2V33xHP413172VkAAAA4FyEWAAAAjhPwIVbbCQKtlQAAAAD+FfAhFgAAAIiOEAsAAADHIcQCAADAcQixAAAAcBxCLAAAAByHEAsAAADHIcQCAADAcVL6ewBO5b5a782rN/w9FMARbv55xd9DABx1iVMEloioKJ///4f/Bbn4NO7KiRMnpECBAv4eBgAASELHjx+X/Pnz+3sYIMTevVu3bsnJkyclQ4YMEhQUJE52+fJlC+T6DzNjxoz+Hg7+F59LYOJzCUx8LoEpOX0uGpeuXLkiefPmlRQp6MYMBLQT3CX9C5zcvonpf2Cc/h+Z5IjPJTDxuQQmPpfAlFw+l0yZMvl7CPDCVwkAAAA4DiEWAAAAjkOIhYSEhMjo0aPtTwQOPpfAxOcSmPhcAhOfC+4lJnYBAADAcajEAgAAwHEIsQAAAHAcQiwAAAAchxALAAAAxyHE3sc2bdokLVq0sKuP6FXHli1b5u8hQUQmTpwo1apVs6vB5cyZU1q1aiWHDh3y97DuezNnzpQKFSp4Fm2vXr26rFq1yt/DgpdXXnnF/lv2wgsv+Hso970xY8bYZ+F9K1WqlL+HhWSGEHsfu3r1qlSsWFHefvttfw8FXjZu3Ch9+vSRHTt2yNq1a+Wvv/6Shg0b2ucF/9Er9GlI+uabbyQ8PFzq1q0rLVu2lAMHDvh7aBCRXbt2yaxZs+yLBgJD2bJl5dSpU57bli1b/D0kJDNcdvY+1qRJE7shsISFhfncnz9/vlVkNTw9+uijfhvX/U7PWnh7+eWXrTqrXzb0/6zhPxEREfL000/LnDlz5KWXXvL3cPC/UqZMKblz5/b3MJCMUYkFAtylS5fsz6xZs/p7KPhfN2/elEWLFll1XNsK4F965qJZs2ZSv359fw8FXn766SdrVytatKh9yTh27Ji/h4RkhkosEMBu3bpl/X01a9aUcuXK+Xs49719+/ZZaL1+/bqEhobK0qVLpUyZMv4e1n1Nv0zs3r3b2gkQOB566CE7i1SyZElrJRg7dqw88sgjsn//fuv3BxIDIRYI8AqT/kefXrLAoP+HvHfvXquOL1myRDp37mw9zARZ/zh+/LgMGDDAesfTpEnj7+HAi3ermvYpa6gtVKiQLF68WLp37+7XsSH5IMQCAapv377y+eef2yoSOqkI/pc6dWopVqyY/VylShWr/k2bNs0mFCHpaZ/42bNnpXLlyj6tHvpv5q233pLIyEgJDg726xjxPzJnziwlSpSQw4cP+3soSEYIsUCAcblc0q9fPztVvWHDBilSpIi/h4TbtHtoUIJ/1KtXz1o8vHXt2tWWchoyZAgBNsAm3x05ckSeeeYZfw8FyQgh9j7/j4r3t+Kff/7ZTpXqBKKCBQv6dWz3ewvBhx9+KMuXL7fesdOnT9vjmTJlkrRp0/p7ePetoUOH2ilS/bdx5coV+4z0S8bq1av9PbT7lv77iN4rnj59esmWLRs95H42ePBgW9FDWwhOnjwpo0ePti8V7du39/fQkIwQYu9jutblP//5T8/9QYMG2Z/a56cN+fAPXbZJ1alTx+fxefPmSZcuXfw0Kuhp606dOtkkFf1CoX1+GmAbNGjg76EBAefEiRMWWC9cuCA5cuSQWrVq2XJ0+jOQWIJceu4SAAAAcBDWiQUAAIDjEGIBAADgOIRYAAAAOA4hFgAAAI5DiAUAAIDjEGIBAADgOIRYAAAAOA4hFgAAAI5DiAXgeL/88osEBQXZZZPjS69+1qpVq9tuo1dNe+GFF+ReKVy4sEydOvWeHR8AkjNCLADHK1CggF0Otly5cuIku3btkp49e3ruaxBftmxZoj7HrVu3JGPGjPLjjz/a/RIlSsimTZsS9TkAwB9S+uVZASCR3LhxQ1KnTi25c+cWp0mK68jv379f0qRJY+H1zJkz8uuvv0q1atXu+fMCwL1GJRZAkpg9e7bkzZvXKoPeWrZsKd26dbOfjxw5Yvdz5coloaGhFrbWrVsX4xT8+PHjpVOnTlZh1Epm9HaCmzdvSvfu3aVIkSKSNm1aKVmypEybNi3WcY0dO9bCpB7rueees1Acl8jISBk8eLDky5dP0qdPLw899JBs2LAhzu1dLpeMGTNGChYsKCEhIfb6+/fvH2s7gf6sWrduba/FfV8tX75cKleubGG0aNGiNuaoqCiJj23btkmNGjXs5y1btsgDDzxg7wkAOB2VWABJom3bttKvXz/56quvpF69evbYxYsXJSwsTL744gu7HxERIU2bNpWXX37ZQt+CBQukRYsWcujQIQuCbpMnT5ZRo0bJ6NGjY30uDcr58+eXjz/+WLJly2ZBTsNunjx5pF27dp7t1q9fb8FQg6gG4a5du9r2+vyx6du3r3z//feyaNEiC6RLly6Vxo0by759+6R48eIxtv/kk0/kjTfesO3Lli0rp0+flm+//TbO1oKcOXPKvHnz7JjBwcH2+ObNmy2wv/nmm/LII49Y0He3IMT1+lXmzJntz+vXr1uY1vsawjXg68+1atWSzz//PM79ASDguQAgibRs2dLVrVs3z/1Zs2a58ubN67p582ac+5QtW9Y1ffp0z/1ChQq5WrVq5bPNzz//7NL/nO3ZsyfO4/Tp08f1+OOPe+537tzZlTVrVtfVq1c9j82cOdMVGhrqGU/t2rVdAwYMsJ9//fVXV3BwsOu3337zOW69evVcQ4cOjfU5p0yZ4ipRooTrxo0bsf5eX8sbb7zhua+vYenSpTGOP2HCBJ/HFi5c6MqTJ4/rdvQ9OXr0qCtLliyuVatW2f3ixYu7PvjgA/v51KlTt90fAAId7QQAkszTTz9t1UmtCKoPPvhAnnrqKUmRIoWnEqun60uXLm3VQm0pOHjwoBw7dsznOFWrVr3jc7399ttSpUoVaxXQ42g7Q/TjVKxYUdKlS+e5X716dRvD8ePHYxxPq61axdTeUj2e+7Zx40arjsZVff7zzz+tBaBHjx5WuY1vG4CbVm7HjRvn85x6LJ3Idu3atTj303aEc+fO2evTym7KlCnl5MmT8vjjj9vvnNhDDADeaCcAkGS0NUALjitXrrR+Vz1Vrqfb3TTArl271toFihUrZr2bTzzxRIw+Ve1HvR09fa/HmjJligXTDBkyyKRJk2Tnzp13PXYNt3qK/5tvvvGc6nfTYBnXqgnaCqF9vfq6evfubePQ4JsqVap4P6/2wLZp0ybG77QVIjZNmjSx91YDs950fBrA9cuDtku4jwsATkaIBZBkNHRpGNMK7OHDh23ClU5Yctu6daut36qTm9xBS3tVE0qPo5OZNDS6xVYt1SqnVkrdE5127NhhgU/DZ3Q6IUqD4NmzZ603Nb702Bre9danTx8pVaqUVXW9X7ebBlt9Dm+6nQZhDfXx9d///tdeV+fOne391slyGur1uZ999tl4HwcAAhkhFkCStxQ0b95cDhw4IB07dvT5nU6O+vTTTy3w6Qz9kSNHxljNID70ODopbPXq1bZCwcKFC23ilP7sTSu8uorBiBEjLCzrRCmdvOVub/CmbQQ6dp1kpRVeDbV6ul4nh1WoUEGaNWsWY5/58+dbKNVVDPS0/vvvv2+htlChQrGOW0/z6/Fq1qxpE9uyZMliE9j0/dKJbVqV1rFp+Nals1566aVYj6OrJ2gF9rvvvrPn1NetPw8ZMiRBYRgAAhk9sQCSVN26dSVr1qxWXezQoYPP715//XULblpF1SDbqFGjWCuWd9KrVy+rQD755JMWIC9cuOBTlXXTVRI08D766KO27WOPPWZLYsVFVw7QEPuvf/3Lqsh6xS8Nx94rJ3jTvt45c+ZYKNWgq20FK1as8JzSj07DsbYdaCVYQ7LS90BXEVizZo21YDz88MPWghFXEHYLDw+359cAe+LECVsjNj69xADgFEE6u8vfgwAAAAASgkosAAAAHIcQCwAAAMchxAIAAMBxCLEAAABwHEIsAAAAHIcQCwAAAMchxAIAAMBxCLEAAABwHEIsAAAAHIcQCwAAAMchxAIAAMBx/j8/R9NJ6TbYRgAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "from matplotlib.colors import ListedColormap\n", + "from matplotlib.patches import Patch\n", + "\n", + "# Read the reduced FASTA\n", + "names, rows = [], []\n", + "for block in open(\"snps.fasta\").read().split(\">\")[1:]:\n", + " head, seq = block.split(chr(10), 1)\n", + " names.append(head.strip())\n", + " rows.append(seq.replace(chr(10), \"\").strip())\n", + "\n", + "code = {\"A\": 0, \"C\": 1, \"G\": 2, \"T\": 3, \"N\": 4, \"-\": 4}\n", + "mat = [[code.get(b, 4) for b in r] for r in rows]\n", + "\n", + "colors = [\"#2ecc71\", \"#3498db\", \"#f1c40f\", \"#e74c3c\", \"#bdc3c7\"]\n", + "cmap = ListedColormap(colors)\n", + "\n", + "fig, ax = plt.subplots(figsize=(7, 3.2))\n", + "ax.imshow(mat, cmap=cmap, vmin=0, vmax=4, aspect='auto')\n", + "ax.set_yticks(range(len(names)), names)\n", + "ax.set_xticks(range(len(rows[0])), range(1, len(rows[0]) + 1))\n", + "ax.set_xlabel(\"variable site #\")\n", + "ax.set_title(\"SNPick reduced alignment\")\n", + "legend = [Patch(facecolor=c, label=b) for b, c in zip(\"ACGT\", colors)]\n", + "legend.append(Patch(facecolor=colors[4], label=\"N / -\"))\n", + "ax.legend(handles=legend, bbox_to_anchor=(1.01, 1), loc=\"upper left\", frameon=False)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "fc968876", + "metadata": {}, + "source": "## 5. Include gaps\n\nBy default gaps (`-`) are ignored and never make a column variable. With `-g` a gap becomes a 5th allele β€” in the VCF an ALT gap is written `*`, while a gap in the reference is written `N` (`*` is not a legal REF allele)." + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "4339ae5d", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-18T16:42:28.487899Z", + "iopub.status.busy": "2026-07-18T16:42:28.487801Z", + "iopub.status.idle": "2026-07-18T16:42:28.893575Z", + "shell.execute_reply": "2026-07-18T16:42:28.892888Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Without -g:\n", + "[snpick] 0 variable, 8 constant (A:2 C:2 G:2 T:2), 0 ambiguous-only, 8 total.\r\n", + "[snpick] No variable positions β€” writing empty output.\r\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "With -g:\n", + "[snpick] 1 variable, 7 constant (A:2 C:2 G:1 T:2), 0 ambiguous-only, 8 total.\r\n" + ] + } + ], + "source": [ + "gapped = {\n", + " \"ref\": \"ACGTACGT\",\n", + " \"sample1\": \"AC-TACGT\",\n", + " \"sample2\": \"ACGTACGT\",\n", + "}\n", + "with open(\"gaps.fasta\", \"w\") as fh:\n", + " for name, seq in gapped.items():\n", + " fh.write(f\">{name}\\n{seq}\\n\")\n", + "\n", + "print(\"Without -g:\")\n", + "!snpick -f gaps.fasta -o nogap.fasta 2>&1 | grep variable\n", + "print(chr(10) + \"With -g:\")\n", + "!snpick -f gaps.fasta -o withgap.fasta -g 2>&1 | grep variable" + ] + }, + { + "cell_type": "markdown", + "id": "c86b864d", + "metadata": {}, + "source": [ + "## Next steps\n", + "\n", + "- [Usage reference](../../usage/) β€” every flag, with examples.\n", + "- [Output formats](../../output/) β€” the reduced FASTA, VCF v4.2 and `fconst` in detail.\n", + "- [Architecture](../../architecture/) β€” how the two-pass, memory-mapped scan works.\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/docs/usage.md b/docs/usage.md new file mode 100644 index 0000000..d32f0ca --- /dev/null +++ b/docs/usage.md @@ -0,0 +1,120 @@ +--- +tags: + - cli + - usage +--- + +# Usage + +``` +snpick [OPTIONS] --fasta +``` + +`-f/--fasta` is always required. `-o/--output` is required **unless** `--dry-run` or `--check` +is given. Use `-` in place of a path for `--fasta` (read from stdin) or for `--output`, +`--vcf-output`, `--sites-output` and `--stats-json` (write to stdout). **At most one** output may +use `-` at a time β€” snpick errors (exit 2) if two or more `-` sinks are given. Run +`snpick --help` for the authoritative, always-current list. + +## Core + +| Argument | Description | +|---|---| +| `-f, --fasta ` | Input FASTA alignment (`-` = stdin; gzip/bgzip auto-detected) | +| `-o, --output ` | Output FASTA of variable sites (`-` = stdout) | +| `-g, --include-gaps` | Treat gaps (`-`) as a 5th character | +| `--format ` | Reduced-alignment format: `fasta` (default), `phylip`, `nexus` | +| `-t, --threads ` | Threads for the parallel scan (default: all logical cores) | +| `-q, --quiet` | Silence progress logs (errors still shown) | +| `-v`, `-vv` | Increase detail (diagnostics; allele-class histogram) | + +## VCF & coordinates + +| Argument | Description | +|---|---| +| `--vcf` | Write a VCF named after the output (`snps.fasta` β†’ `snps.vcf`) | +| `--vcf-output ` | Write the VCF to a custom path (implies `--vcf`; `-` = stdout) | +| `--chrom ` | `CHROM` / `##contig` name (default `1`, e.g. `NC_000962.3`) | +| `--reference ` | Sequence ID used for REF/polarity (default: first) | +| `--ref-coords` | Write `POS` as the ungapped reference position, not the alignment column | +| `--sites-output ` | Map each site: `alignment_pos`, `ref_pos`, `ref`, `alt` (`-` = stdout) | + +## Filtering & masking + +| Argument | Description | +|---|---| +| `--max-missing ` | Drop sites with more than fraction `f` of missing genotypes | +| `--mac ` | Drop sites with minor-allele count below `n` (e.g. `2` drops singletons) | +| `--maf ` | Drop sites with minor-allele frequency below `f` | +| `--min-samples ` | Drop sites with fewer than `n` samples with data | +| `--max-alleles ` | Drop sites with more than `n` distinct alleles (`2` = biallelic) | +| `--keep-samples ` | Keep only these samples (comma-separated, or `@file`) | +| `--exclude-samples ` | Drop these samples (mutually exclusive with `--keep-samples`) | +| `--mask ` | Mask BED regions (excluded from output **and** `fconst`) | +| `--mask-ref` | Interpret `--mask` coordinates as reference positions | + +!!! info "Filtering keeps ASC valid" + Per-site filters remove **variable** sites from the output only; they are never reclassified + as constant, so the `fconst` counts are untouched. Sample selection filters *before* the + scan, so `fconst` is recomputed for exactly the retained samples. + +## QC, reporting & robustness + +| Argument | Description | +|---|---| +| `--stats-json ` | Machine-readable JSON summary (`-` = stdout) | +| `--dry-run` | Report statistics without writing any FASTA/VCF | +| `--check` | Audit the alignment composition and exit | +| `--on-invalid ` | `ignore` (default), `warn`, or `error` on out-of-alphabet bytes | +| `--iupac-mode ` | `missing` (default) or `resolve` IUPAC codes to bases | +| `--allow-dup-ids` | Permit duplicate sequence IDs instead of erroring | + +**Exit codes:** `0` success Β· `1` I/O error Β· `2` bad input/data. + +## Examples + +### Basic extraction + +```bash +snpick -f alignment.fasta -o snps.fasta +``` + +### VCF with reference coordinates and a stats sidecar + +```bash +snpick -f alignment.fasta.gz -o snps.fasta \ + --vcf --chrom NC_000962.3 --reference H37Rv --ref-coords \ + --stats-json stats.json # (1)! +``` + +1. `--ref-coords` sets `POS` (and the contig length) to ungapped reference positions; + `--reference` pins REF polarity to H37Rv; `--stats-json` gives IQ-TREE the `fconst` array + without scraping stderr. + +### Filter and mask + +```bash +# Drop singletons, keep biallelic sites, mask repetitive regions (reference coords) +snpick -f alignment.fasta -o snps.fasta \ + --mac 2 --max-alleles 2 --mask exclude.bed --mask-ref +``` + +### Subset samples + +```bash +snpick -f alignment.fasta -o clade.fasta --keep-samples @clade_ids.txt +``` + +### Alternative formats and piping + +```bash +snpick -f alignment.fasta -o snps.nex --format nexus +gzip -dc big.fasta.gz | snpick -f - -o - --vcf-output snps.vcf > snps.fasta +``` + +### Dry run and composition audit + +```bash +snpick -f alignment.fasta --dry-run --stats-json - # stats only, JSON to stdout +snpick -f alignment.fasta --check # composition breakdown, then exit +``` diff --git a/fuzz/Cargo.toml b/fuzz/Cargo.toml new file mode 100644 index 0000000..772114e --- /dev/null +++ b/fuzz/Cargo.toml @@ -0,0 +1,24 @@ +[package] +name = "snpick-fuzz" +version = "0.0.0" +publish = false +edition = "2021" + +[package.metadata] +cargo-fuzz = true + +# Own workspace so it does not interfere with the parent package's build. +[workspace] + +[dependencies] +libfuzzer-sys = "0.4" + +[dependencies.snpick] +path = ".." + +[[bin]] +name = "index_fasta" +path = "fuzz_targets/index_fasta.rs" +test = false +doc = false +bench = false diff --git a/fuzz/fuzz_targets/index_fasta.rs b/fuzz/fuzz_targets/index_fasta.rs new file mode 100644 index 0000000..3318789 --- /dev/null +++ b/fuzz/fuzz_targets/index_fasta.rs @@ -0,0 +1,15 @@ +#![no_main] +//! Fuzz the FASTA indexer: `index_fasta` does manual index arithmetic over untrusted, +//! machine-generated input, so it must never panic β€” only return `Ok`/`Err`. +//! +//! Run with a nightly toolchain and cargo-fuzz: +//! cargo +nightly fuzz run index_fasta + +use libfuzzer_sys::fuzz_target; + +fuzz_target!(|data: &[u8]| { + if let Ok((records, seq_len, layout)) = snpick::fasta::index_fasta(data) { + // Exercise the reference reader on the first record too. + let _ = snpick::fasta::get_ref_seq(data, &records[0], seq_len, layout); + } +}); diff --git a/mkdocs.yml b/mkdocs.yml new file mode 100644 index 0000000..e91f4c1 --- /dev/null +++ b/mkdocs.yml @@ -0,0 +1,163 @@ +site_name: SNPick +site_description: Fast, memory-efficient extraction of variable sites from FASTA alignments +site_url: https://pathogenomics-lab.github.io/snpick/ +site_author: Paula Ruiz-Rodriguez +repo_url: https://github.com/PathoGenOmics-Lab/snpick +repo_name: PathoGenOmics-Lab/snpick +edit_uri: edit/main/docs/ +copyright: Copyright © Paula Ruiz-Rodriguez & Mireia Coscolla β€” IΒ²SysBio (UV-CSIC) + +theme: + name: material + custom_dir: overrides + logo: assets/logo.svg + favicon: assets/logo.png + language: en + palette: + - media: "(prefers-color-scheme)" + toggle: + icon: material/brightness-auto + name: Switch to light mode + - media: "(prefers-color-scheme: light)" + scheme: default + primary: deep purple + accent: pink + toggle: + icon: material/brightness-7 + name: Switch to dark mode + - media: "(prefers-color-scheme: dark)" + scheme: slate + primary: deep purple + accent: pink + toggle: + icon: material/brightness-4 + name: Switch to system preference + font: + text: Inter + code: JetBrains Mono + features: + - announce.dismiss + - navigation.instant + - navigation.instant.prefetch + - navigation.instant.progress + - navigation.tabs + - navigation.tabs.sticky + - navigation.sections + - navigation.indexes + - navigation.top + - navigation.footer + - navigation.tracking + - navigation.path + - toc.follow + - search.suggest + - search.highlight + - search.share + - content.code.copy + - content.code.annotate + - content.code.select + - content.tabs.link + - content.tooltips + - content.action.edit + - content.action.view + icon: + repo: fontawesome/brands/github + edit: material/pencil + view: material/eye + +nav: + - Home: index.md + - Installation: installation.md + - Usage: usage.md + - Tutorial: tutorials/tutorial.ipynb + - Output formats: output.md + - Benchmarks: benchmarks.md + - Architecture: architecture.md + - Contributing: contributing.md + - Changelog: changelog.md + - Tags: tags.md + +markdown_extensions: + - abbr + - admonition + - attr_list + - def_list + - footnotes + - md_in_html + - tables + - toc: + permalink: true + title: On this page + - pymdownx.betterem + - pymdownx.caret + - pymdownx.critic + - pymdownx.details + - pymdownx.highlight: + anchor_linenums: true + line_spans: __span + pygments_lang_class: true + - pymdownx.inlinehilite + - pymdownx.keys + - pymdownx.mark + - pymdownx.smartsymbols + - pymdownx.snippets: + base_path: [".", "docs"] + auto_append: ["includes/abbreviations.md"] + check_paths: true + - pymdownx.superfences: + custom_fences: + - name: mermaid + class: mermaid + format: !!python/name:pymdownx.superfences.fence_code_format + - pymdownx.tabbed: + alternate_style: true + slugify: !!python/object/apply:pymdownx.slugs.slugify + kwds: + case: lower + - pymdownx.tasklist: + custom_checkbox: true + - pymdownx.tilde + - pymdownx.emoji: + emoji_index: !!python/name:material.extensions.emoji.twemoji + emoji_generator: !!python/name:material.extensions.emoji.to_svg + +plugins: + - search + - tags + - mkdocs-jupyter: + execute: false + include_source: true + ignore_h1_titles: true + - glightbox: + touchNavigation: true + loop: false + effect: zoom + draggable: true + - git-revision-date-localized: + type: timeago + enable_creation_date: true + enable_git_follow: false + fallback_to_build_date: true + - social: + # Social cards need Cairo/Pango; generate them only in CI (where the libs + # are installed), so local builds don't require the native dependencies. + enabled: !ENV [CI, false] + cards_layout_options: + background_color: "#5e35b1" + - minify: + minify_html: true + +extra: + social: + - icon: fontawesome/brands/github + link: https://github.com/PathoGenOmics-Lab + name: PathoGenOmics-Lab on GitHub + - icon: fontawesome/brands/python + link: https://anaconda.org/bioconda/snpick + name: SNPick on Bioconda + +extra_css: + - stylesheets/extra.css + +# The abbreviations glossary is auto-appended via snippets, not a standalone page. +exclude_docs: | + includes/ diff --git a/overrides/main.html b/overrides/main.html new file mode 100644 index 0000000..98d292a --- /dev/null +++ b/overrides/main.html @@ -0,0 +1,10 @@ +{% extends "base.html" %} + +{% block announce %} + + + {% include ".icons/material/party-popper.svg" %} + + SNPick 1.0.2 is here β€” see what changed in the changelog. + +{% endblock %} diff --git a/snpick.def b/snpick.def new file mode 100644 index 0000000..8c6906c --- /dev/null +++ b/snpick.def @@ -0,0 +1,17 @@ +Bootstrap: docker +From: ghcr.io/pathogenomics-lab/snpick:latest + +# Apptainer / Singularity definition. Build with: +# apptainer build snpick.sif snpick.def +# The published container image already ships a static snpick binary. + +%labels + Author Paula Ruiz-Rodriguez + Description Fast extraction of variable sites from FASTA alignments + +%runscript + exec /usr/local/bin/snpick "$@" + +%help + snpick β€” extract variable (SNP) sites from FASTA alignments. + Usage: apptainer run snpick.sif -f alignment.fasta -o snps.fasta diff --git a/src/audit.rs b/src/audit.rs new file mode 100644 index 0000000..16ea091 --- /dev/null +++ b/src/audit.rs @@ -0,0 +1,83 @@ +//! Alignment composition audit: a per-byte histogram over the sequence data, used by +//! `--on-invalid` and `--check`. Only run on demand (it is a full extra pass). + +use crate::fasta::FastaRecord; +use crate::types::SeqLayout; + +/// Count every sequence byte (newlines skipped) into a 256-entry histogram. +pub fn composition( + data: &[u8], records: &[FastaRecord], seq_length: usize, layout: SeqLayout, +) -> [u64; 256] { + let mut h = [0u64; 256]; + for rec in records { + if layout.single_line { + for &b in &data[rec.seq_offset..rec.seq_offset + seq_length] { + h[b as usize] += 1; + } + } else { + let mut pos = rec.seq_offset; + let end = data.len(); + let mut i = 0; + while i < seq_length && pos < end { + let b = data[pos]; + pos += 1; + if b == b'\n' || b == b'\r' { + continue; + } + h[b as usize] += 1; + i += 1; + } + } + } + h +} + +/// Whether a byte is part of the accepted nucleotide alphabet (ACGTU, N, IUPAC ambiguity +/// codes, and gap characters). Anything else (e.g. protein residues) is "invalid". +pub fn is_valid_nucleotide(b: u8) -> bool { + matches!( + b.to_ascii_uppercase(), + b'A' | b'C' | b'G' | b'T' | b'U' | b'N' + | b'R' | b'Y' | b'S' | b'W' | b'K' | b'M' | b'B' | b'D' | b'H' | b'V' + | b'-' | b'.' | b'?' + ) +} + +/// Total count of out-of-alphabet bytes. +pub fn invalid_count(h: &[u64; 256]) -> u64 { + (0..256u16).filter(|&i| !is_valid_nucleotide(i as u8)).map(|i| h[i as usize]).sum() +} + +/// A short human-readable composition breakdown (for `--check`). +pub fn summary(h: &[u64; 256]) -> String { + let sum = |bytes: &[u8]| -> u64 { bytes.iter().map(|&b| h[b as usize]).sum() }; + let total: u64 = h.iter().sum(); + let acgt = sum(b"ACGTacgt"); + let n = sum(b"Nn"); + let gap = sum(b"-."); + let iupac = sum(b"RYSWKMBDHVryswkmbdhvUu"); + let invalid = invalid_count(h); + let pct = |c: u64| if total == 0 { 0.0 } else { 100.0 * c as f64 / total as f64 }; + format!( + "total={} A/C/G/T={} ({:.2}%) N={} ({:.2}%) gap={} ({:.2}%) IUPAC={} ({:.2}%) invalid={} ({:.2}%)", + total, acgt, pct(acgt), n, pct(n), gap, pct(gap), iupac, pct(iupac), invalid, pct(invalid), + ) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn invalid_detection() { + let mut h = [0u64; 256]; + h[b'A' as usize] = 10; + h[b'N' as usize] = 2; + h[b'-' as usize] = 1; + h[b'E' as usize] = 3; // protein residue -> invalid + assert_eq!(invalid_count(&h), 3); + assert!(is_valid_nucleotide(b'R')); + assert!(!is_valid_nucleotide(b'E')); + assert!(summary(&h).contains("invalid=3")); + } +} diff --git a/src/coords.rs b/src/coords.rs new file mode 100644 index 0000000..40d7665 --- /dev/null +++ b/src/coords.rs @@ -0,0 +1,127 @@ +//! Reference-anchored coordinates and BED region masking. + +use std::io; + +/// Ungapped reference position (1-based) for each alignment column. A column where the +/// reference base is a gap inherits the position of the preceding non-gap base. +pub fn ref_positions(ref_seq: &[u8]) -> Vec { + let mut out = Vec::with_capacity(ref_seq.len()); + let mut rp: u32 = 0; + for &b in ref_seq { + if b != b'-' { + rp += 1; + } + out.push(rp); + } + out +} + +/// Parse a BED file into (start, end) 0-based half-open intervals (the chrom column is ignored; +/// snpick works over a single alignment). +pub fn parse_bed(path: &str) -> io::Result> { + let content = std::fs::read_to_string(path).map_err(|e| { + io::Error::new(e.kind(), format!("Cannot read BED '{}': {}", path, e)) + })?; + let mut ivs = Vec::new(); + for (n, line) in content.lines().enumerate() { + let line = line.trim(); + if line.is_empty() || line.starts_with('#') { + continue; + } + // Skip BED header lines, but only when the FIRST token is exactly `track`/`browser` + // (not a chrom that merely starts with those letters, e.g. `track_scaffold_1`). + let first = line.split_whitespace().next().unwrap_or(""); + if first == "track" || first == "browser" { + continue; + } + let mut f = line.split('\t'); + let _chrom = f.next(); + let start = f.next().and_then(|s| s.trim().parse::().ok()); + let end = f.next().and_then(|s| s.trim().parse::().ok()); + match (start, end) { + (Some(s), Some(e)) if s <= e => ivs.push((s, e)), + _ => { + return Err(io::Error::new(io::ErrorKind::InvalidData, format!( + "Malformed BED at line {}: expected 'chromstartend' with start <= end.", + n + 1 + ))) + } + } + } + Ok(ivs) +} + +/// Build a per-column boolean mask (true = masked). When `ref_coords` is `Some`, the intervals +/// are 0-based half-open reference positions; otherwise they are alignment columns. +pub fn build_mask( + intervals: &[(usize, usize)], seq_length: usize, ref_coords: Option<&[u32]>, +) -> Vec { + let mut mask = vec![false; seq_length]; + match ref_coords { + None => { + for &(s, e) in intervals { + let (s, e) = (s.min(seq_length), e.min(seq_length)); + for m in &mut mask[s..e] { + *m = true; + } + } + } + Some(rp) => { + for &(s, e) in intervals { + // 0-based half-open [s, e) over the reference == 1-based positions s+1 ..= e. + // Compare in usize (saturating on the +1) so an out-of-range or degenerate + // interval never wraps or truncates to u32 and masks the wrong columns β€” the + // None branch above is likewise robust via `.min(seq_length)`. Reference + // positions are always <= MAX_SEQ_LENGTH, so widening `p` loses nothing. + let lo = s.saturating_add(1); + for (col, &p) in rp.iter().enumerate() { + let p = p as usize; + if p >= lo && p <= e { + mask[col] = true; + } + } + } + } + } + mask +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn ref_positions_gaps_inherit() { + // A - C G -> positions 1 1 2 3 + assert_eq!(ref_positions(b"A-CG"), vec![1, 1, 2, 3]); + assert_eq!(ref_positions(b"ACGT"), vec![1, 2, 3, 4]); + } + + #[test] + fn mask_alignment_columns() { + // columns [1,3) masked + let m = build_mask(&[(1, 3)], 5, None); + assert_eq!(m, vec![false, true, true, false, false]); + } + + #[test] + fn mask_reference_coordinates() { + // ref A-CG -> positions [1,1,2,3]; mask ref [1,2) = ref position 2 -> column 2 + let rp = ref_positions(b"A-CG"); + let m = build_mask(&[(1, 2)], 4, Some(&rp)); + assert_eq!(m, vec![false, false, true, false]); + } + + #[test] + fn mask_reference_out_of_range_is_noop() { + // Reference positions are small; an interval far beyond them (or a degenerate empty one + // at usize::MAX) must mask nothing. The old `as u32` cast wrapped/truncated these into + // the valid range and masked the wrong columns (and panicked on the +1 in debug). + let rp = ref_positions(b"ACGT"); // positions [1,2,3,4] + let huge = 1usize << 40; // ≑ 0 mod 2^32 + assert_eq!(build_mask(&[(huge, huge + 4)], 4, Some(&rp)), vec![false; 4]); + assert_eq!(build_mask(&[(usize::MAX, usize::MAX)], 4, Some(&rp)), vec![false; 4]); + // An in-range interval still masks correctly. + assert_eq!(build_mask(&[(1, 3)], 4, Some(&rp)), vec![false, true, true, false]); + } +} diff --git a/src/extract.rs b/src/extract.rs index 7c907fd..6c0e82f 100644 --- a/src/extract.rs +++ b/src/extract.rs @@ -9,6 +9,36 @@ use std::io::{self, BufWriter, Write}; use crate::fasta::FastaRecord; use crate::types::*; +/// Write a NEXUS taxon label, single-quoting (and doubling embedded quotes) when it contains +/// characters NEXUS treats as punctuation/whitespace, so labels like `iso[London]` survive. +fn write_nexus_label(w: &mut impl Write, id: &[u8]) -> io::Result<()> { + let safe = !id.is_empty() + && id.iter().all(|&b| b.is_ascii_alphanumeric() || b == b'_' || b == b'.'); + if safe { + return w.write_all(id); + } + w.write_all(b"'")?; + for &b in id { + if b == b'\'' { + w.write_all(b"''")?; + } else { + w.write_all(&[b])?; + } + } + w.write_all(b"'") +} + +/// Open an output sink: a file, or stdout when the path is "-". +pub fn open_sink(path: &str) -> io::Result> { + if path == "-" { + Ok(Box::new(io::stdout().lock())) + } else { + let f = File::create(path).map_err(|e| io::Error::new(e.kind(), + format!("Cannot create output '{}': {}", path, e)))?; + Ok(Box::new(f)) + } +} + /// Parameters for variable site extraction (pass 2). pub struct ExtractParams<'a> { pub records: &'a [FastaRecord<'a>], @@ -17,6 +47,8 @@ pub struct ExtractParams<'a> { pub lookup: &'a [u8; 256], pub upper: &'a [u8; 256], pub layout: SeqLayout, + pub format: OutputFormat, + pub progress: bool, } /// Pass 2: extract variable sites from alignment and write output FASTA. @@ -28,16 +60,29 @@ pub struct ExtractParams<'a> { pub fn pass2_extract( data: &[u8], var_positions: &mut [VariablePosition], params: &ExtractParams<'_>, ) -> io::Result>> { - let ExtractParams { records, output, collect_vcf, lookup, upper, layout } = params; + let ExtractParams { records, output, collect_vcf, lookup, upper, layout, format, progress } = params; let collect_vcf = *collect_vcf; let layout = *layout; + let format = *format; + let progress = *progress; let num_var = var_positions.len(); let num_samples = records.len(); let pos_indices: Vec = var_positions.iter().map(|v| v.index).collect(); - let out_file = File::create(output).map_err(|e| io::Error::new(e.kind(), - format!("Cannot create output '{}': {}", output, e)))?; - let mut writer = BufWriter::with_capacity(IO_BUF, out_file); + let mut writer = BufWriter::with_capacity(IO_BUF, open_sink(output)?); + + // Format preamble. + match format { + OutputFormat::Fasta => {} + OutputFormat::Phylip => writeln!(writer, "{} {}", num_samples, num_var)?, + OutputFormat::Nexus => { + writeln!(writer, "#NEXUS")?; + writeln!(writer, "BEGIN DATA;")?; + writeln!(writer, " DIMENSIONS NTAX={} NCHAR={};", num_samples, num_var)?; + writeln!(writer, " FORMAT DATATYPE=DNA MISSING=N GAP=-;")?; + writeln!(writer, " MATRIX")?; + } + } let mut vcf_geno: Vec = if collect_vcf { vec![0u8; num_var * num_samples] } else { Vec::new() }; let mut ns_counts: Vec = if collect_vcf { vec![0usize; num_var] } else { Vec::new() }; @@ -66,15 +111,32 @@ pub fn pass2_extract( } } - writer.write_all(b">")?; - writer.write_all(rec.id)?; - if !rec.desc.is_empty() { - writer.write_all(b" ")?; - writer.write_all(rec.desc)?; + match format { + OutputFormat::Fasta => { + writer.write_all(b">")?; + writer.write_all(rec.id)?; + if !rec.desc.is_empty() { + writer.write_all(b" ")?; + writer.write_all(rec.desc)?; + } + writer.write_all(b"\n")?; + writer.write_all(&var_buf)?; + writer.write_all(b"\n")?; + } + OutputFormat::Phylip => { + writer.write_all(rec.id)?; + writer.write_all(b" ")?; + writer.write_all(&var_buf)?; + writer.write_all(b"\n")?; + } + OutputFormat::Nexus => { + writer.write_all(b" ")?; + write_nexus_label(&mut writer, rec.id)?; + writer.write_all(b" ")?; + writer.write_all(&var_buf)?; + writer.write_all(b"\n")?; + } } - writer.write_all(b"\n")?; - writer.write_all(&var_buf)?; - writer.write_all(b"\n")?; if collect_vcf { for (vi, &nuc) in var_buf.iter().enumerate() { @@ -84,6 +146,18 @@ pub fn pass2_extract( } } } + + if progress && (si % 256 == 0) { + eprint!("\r[snpick] Extracting sequence {}/{}", si + 1, num_samples); + } + } + if progress { + eprintln!("\r[snpick] Extracted {} sequences. ", num_samples); + } + + if let OutputFormat::Nexus = format { + writeln!(writer, " ;")?; + writeln!(writer, "END;")?; } writer.flush()?; @@ -94,6 +168,5 @@ pub fn pass2_extract( } } - eprintln!("[snpick] Pass 2: Wrote {} sequences to {}.", num_samples, output); if collect_vcf { Ok(Some(vcf_geno)) } else { Ok(None) } } diff --git a/src/fasta.rs b/src/fasta.rs index 08f21fd..8813b88 100644 --- a/src/fasta.rs +++ b/src/fasta.rs @@ -57,6 +57,7 @@ pub fn index_fasta(data: &[u8]) -> io::Result<(Vec>, usize, SeqL // Count bases let mut seq_len = 0usize; let mut line_count = 0usize; + let mut has_blank_line = false; while pos < len && data[pos] != b'>' { let line_start = pos; while pos < len && data[pos] != b'\n' && data[pos] != b'\r' { pos += 1; } @@ -64,12 +65,17 @@ pub fn index_fasta(data: &[u8]) -> io::Result<(Vec>, usize, SeqL if line_len > 0 { seq_len += line_len; line_count += 1; + } else { + has_blank_line = true; } if pos < len && data[pos] == b'\r' { pos += 1; } if pos < len && data[pos] == b'\n' { pos += 1; } } - if line_count > 1 { is_single_line = false; } + // A blank line inside the sequence shifts byte offsets, so the single-line + // fast path (data[seq_offset + pos]) would read misaligned bases. Route such + // records through the newline-skipping scanner instead. + if line_count > 1 || has_blank_line { is_single_line = false; } if records.is_empty() { if seq_len == 0 { diff --git a/src/filter.rs b/src/filter.rs new file mode 100644 index 0000000..e8b7954 --- /dev/null +++ b/src/filter.rs @@ -0,0 +1,198 @@ +//! Per-site filtering: recount alleles and missingness over the variable columns. +//! +//! Filtering removes *variable* sites from the output only β€” it never reclassifies a site as +//! constant, so the ASC `fconst` counts (constant sites) are untouched and stay valid. + +use crate::fasta::FastaRecord; +use crate::types::{SeqLayout, BIT_A, BIT_C, BIT_G, BIT_T, BIT_GAP}; + +/// Per-variable-site allele and missingness counts across all samples. +#[derive(Clone, Copy, Default)] +pub struct SiteStat { + /// Counts of A, C, G, T. + pub counts: [u32; 4], + /// Count of gaps ('-'). + pub gap: u32, + /// Count of ambiguous / uncalled bases. + pub missing: u32, +} + +/// Thresholds for per-site filtering. `None` = no constraint. +#[derive(Clone, Copy, Default)] +pub struct SiteFilters { + /// Max allowed fraction of missing genotypes (0.0–1.0). + pub max_missing: Option, + /// Min minor-allele count. + pub mac: Option, + /// Min minor-allele frequency (0.0–1.0). + pub maf: Option, + /// Min number of samples with data. + pub min_samples: Option, + /// Max number of distinct alleles (e.g. 2 = biallelic only). + pub max_alleles: Option, +} + +impl SiteFilters { + /// Whether any threshold is set (skip the counting pass entirely if not). + pub fn active(&self) -> bool { + self.max_missing.is_some() + || self.mac.is_some() + || self.maf.is_some() + || self.min_samples.is_some() + || self.max_alleles.is_some() + } +} + +impl SiteStat { + /// Number of samples with a called allele (i.e. any set bit; gaps count only under + /// `include_gaps`, since the lookup then maps '-' to a bit). + pub fn ns(&self, num_samples: u32) -> u32 { + num_samples - self.missing + } + + /// Counts of the alleles actually present (ACGT with count > 0, plus gap if `include_gaps`). + fn allele_counts(&self, include_gaps: bool) -> Vec { + let mut v: Vec = self.counts.iter().copied().filter(|&c| c > 0).collect(); + if include_gaps && self.gap > 0 { + v.push(self.gap); + } + v + } + + /// Whether this site passes every set threshold. + pub fn passes(&self, f: &SiteFilters, num_samples: u32, include_gaps: bool) -> bool { + let ns = self.ns(num_samples); + let missing = self.missing; + if let Some(ms) = f.min_samples { + if ns < ms { + return false; + } + } + if let Some(mm) = f.max_missing { + if num_samples > 0 && (missing as f64 / num_samples as f64) > mm { + return false; + } + } + let alleles = self.allele_counts(include_gaps); + if let Some(ma) = f.max_alleles { + if alleles.len() as u32 > ma { + return false; + } + } + let mac = alleles.iter().copied().min().unwrap_or(0); + if let Some(m) = f.mac { + if mac < m { + return false; + } + } + if let Some(maf) = f.maf { + if ns == 0 || (mac as f64 / ns as f64) < maf { + return false; + } + } + true + } +} + +#[inline] +fn tally(s: &mut SiteStat, bits: u8) { + // `bits` may carry several flags under --iupac-mode resolve (e.g. R = A|G), so credit every + // set base β€” matching the classifier's resolved allele set. Zero = uncalled/ambiguous. + if bits == 0 { + s.missing += 1; + return; + } + if bits & BIT_A != 0 { s.counts[0] += 1; } + if bits & BIT_C != 0 { s.counts[1] += 1; } + if bits & BIT_G != 0 { s.counts[2] += 1; } + if bits & BIT_T != 0 { s.counts[3] += 1; } + if bits & BIT_GAP != 0 { s.gap += 1; } +} + +/// Count alleles and missing calls at each variable position (sparse pass over just the +/// variable columns). Gap handling follows `lookup` (built with the same `include_gaps`). +pub fn count_sites( + data: &[u8], records: &[FastaRecord], pos_indices: &[usize], + layout: SeqLayout, lookup: &[u8; 256], +) -> Vec { + let num_var = pos_indices.len(); + let mut stats = vec![SiteStat::default(); num_var]; + for rec in records { + if layout.single_line { + let base = rec.seq_offset; + for (vi, &p) in pos_indices.iter().enumerate() { + tally(&mut stats[vi], lookup[data[base + p] as usize]); + } + } else { + let mut pos = rec.seq_offset; + let end = data.len(); + let mut base_idx = 0usize; + let mut var_idx = 0usize; + while var_idx < num_var && pos < end { + let b = data[pos]; + pos += 1; + if b == b'\n' || b == b'\r' { + continue; + } + if base_idx == pos_indices[var_idx] { + tally(&mut stats[var_idx], lookup[b as usize]); + var_idx += 1; + } + base_idx += 1; + } + } + } + stats +} + +#[cfg(test)] +mod tests { + use super::*; + + fn stat(a: u32, c: u32, g: u32, t: u32, gap: u32, missing: u32) -> SiteStat { + SiteStat { counts: [a, c, g, t], gap, missing } + } + + #[test] + fn mac_drops_singletons() { + let f = SiteFilters { mac: Some(2), ..Default::default() }; + // A:9 T:1 -> minor-allele count 1, dropped by --mac 2 + assert!(!stat(9, 0, 0, 1, 0, 0).passes(&f, 10, false)); + // A:8 T:2 -> minor-allele count 2, kept + assert!(stat(8, 0, 0, 2, 0, 0).passes(&f, 10, false)); + } + + #[test] + fn max_missing_and_min_samples() { + // 3 of 10 missing = 0.3 + let s = stat(5, 2, 0, 0, 0, 3); + assert!(s.passes(&SiteFilters { max_missing: Some(0.3), ..Default::default() }, 10, false)); + assert!(!s.passes(&SiteFilters { max_missing: Some(0.2), ..Default::default() }, 10, false)); + assert!(s.passes(&SiteFilters { min_samples: Some(7), ..Default::default() }, 10, false)); + assert!(!s.passes(&SiteFilters { min_samples: Some(8), ..Default::default() }, 10, false)); + } + + #[test] + fn max_alleles_biallelic() { + let f = SiteFilters { max_alleles: Some(2), ..Default::default() }; + assert!(stat(5, 5, 0, 0, 0, 0).passes(&f, 10, false)); // 2 alleles + assert!(!stat(4, 3, 3, 0, 0, 0).passes(&f, 10, false)); // 3 alleles + } + + #[test] + fn resolved_iupac_counts_both_bases() { + // Under --iupac-mode resolve, an R credits both A and G, so a site is 2-allelic and + // has no missing β€” it must fail --max-alleles 1 and pass --max-missing 0. + let s = stat(4, 0, 1, 0, 0, 0); // A=4 (incl. R), G=1 (from R) + assert!(!s.passes(&SiteFilters { max_alleles: Some(1), ..Default::default() }, 4, false)); + assert!(s.passes(&SiteFilters { max_missing: Some(0.0), ..Default::default() }, 4, false)); + assert_eq!(s.ns(4), 4); + } + + #[test] + fn maf_frequency() { + let f = SiteFilters { maf: Some(0.15), ..Default::default() }; + assert!(!stat(9, 0, 0, 1, 0, 0).passes(&f, 10, false)); // maf 0.1 + assert!(stat(8, 0, 0, 2, 0, 0).passes(&f, 10, false)); // maf 0.2 + } +} diff --git a/src/input.rs b/src/input.rs new file mode 100644 index 0000000..7f1ad25 --- /dev/null +++ b/src/input.rs @@ -0,0 +1,158 @@ +//! Input handling: transparent gzip/bgzip decompression and stdin, mapped to memory. +//! +//! The two-pass architecture needs random access, which a pipe or a compressed stream can't +//! provide, so compressed or piped input is spooled to a temp file that is then memory-mapped +//! and removed on drop. + +use std::fs::File; +use std::io::{self, BufReader, BufWriter, Read, Write}; +use std::path::PathBuf; + +use flate2::read::MultiGzDecoder; +use memmap2::Mmap; + +/// A memory-mapped input, owning any temp file created for decompressed/piped data. +pub struct MappedInput { + pub mmap: Mmap, + temp: Option, +} + +impl Drop for MappedInput { + fn drop(&mut self) { + if let Some(p) = &self.temp { + let _ = std::fs::remove_file(p); + } + } +} + +fn temp_path() -> PathBuf { + std::env::temp_dir().join(format!("snpick-{}.tmp", std::process::id())) +} + +/// Read up to 2 magic bytes, tolerating short reads (a pipe's first read may return < 2 bytes). +fn read_magic(mut r: impl Read) -> io::Result<(usize, [u8; 2])> { + let mut buf = [0u8; 2]; + let mut n = 0; + while n < 2 { + match r.read(&mut buf[n..]) { + Ok(0) => break, + Ok(k) => n += k, + Err(ref e) if e.kind() == io::ErrorKind::Interrupted => continue, + Err(e) => return Err(e), + } + } + Ok((n, buf)) +} + +fn spool_to_temp(mut reader: impl Read) -> io::Result<(Mmap, PathBuf)> { + let tp = temp_path(); + let build = (|| -> io::Result { + { + let mut out = BufWriter::new(File::create(&tp)?); + io::copy(&mut reader, &mut out)?; + out.flush()?; + } + let f = File::open(&tp)?; + unsafe { Mmap::map(&f) } + })(); + match build { + Ok(mmap) => Ok((mmap, tp)), + Err(e) => { + // Don't leak the spool file if decompression / mmap failed. + let _ = std::fs::remove_file(&tp); + Err(e) + } + } +} + +/// Map the input for reading. `path == "-"` reads stdin; gzip/bgzip input (detected by magic +/// bytes) is decompressed transparently. +pub fn map_input(path: &str) -> io::Result { + if path == "-" { + // Peek the first two bytes, then chain them back so gzip over stdin is detected. + let mut reader = BufReader::new(io::stdin().lock()); + let (n, magic) = read_magic(&mut reader)?; + let head = std::io::Cursor::new(magic[..n].to_vec()); + let stream = head.chain(reader); + let (mmap, tp) = if n == 2 && magic == [0x1f, 0x8b] { + spool_to_temp(MultiGzDecoder::new(stream))? + } else { + spool_to_temp(stream)? + }; + return Ok(MappedInput { mmap, temp: Some(tp) }); + } + + let f = File::open(path) + .map_err(|e| io::Error::new(e.kind(), format!("Cannot open '{}': {}", path, e)))?; + + // Peek the first two bytes for the gzip magic (0x1f 0x8b); bgzip is a valid gzip stream. + let (n, magic) = read_magic(&f)?; + let is_gzip = n == 2 && magic == [0x1f, 0x8b]; + + if is_gzip { + // Chain the peeked magic back onto the same handle rather than reopening the path. + // Reopening rewinds a regular file but NOT a non-seekable path (process substitution + // `<(gzip -c ...)` or a FIFO), where the magic bytes are already consumed β€” the decoder + // would then start mid-stream and fail with "invalid gzip header". This mirrors the + // stdin branch and works for both seekable and non-seekable inputs. + let head = std::io::Cursor::new(magic[..n].to_vec()); + let stream = head.chain(BufReader::new(f)); + let (mmap, tp) = spool_to_temp(MultiGzDecoder::new(stream))?; + Ok(MappedInput { mmap, temp: Some(tp) }) + } else { + // Mmap maps from offset 0 regardless of the read position left by read_magic. + let mmap = unsafe { Mmap::map(&f)? }; + Ok(MappedInput { mmap, temp: None }) + } +} + +#[cfg(test)] +mod tests { + use super::map_input; + use flate2::write::GzEncoder; + use flate2::Compression; + use std::fs::File; + use std::io::Write; + + fn gzip(bytes: &[u8]) -> Vec { + let mut e = GzEncoder::new(Vec::new(), Compression::fast()); + e.write_all(bytes).unwrap(); + e.finish().unwrap() + } + + #[test] + fn gzip_regular_file_decodes() { + let raw = b">s1\nATGC\n>s2\nATCC\n"; + let p = std::env::temp_dir().join(format!("snpick-gz-{}.fa.gz", std::process::id())); + std::fs::write(&p, gzip(raw)).unwrap(); + let m = map_input(p.to_str().unwrap()).unwrap(); + assert_eq!(&m.mmap[..], raw); + std::fs::remove_file(&p).ok(); + } + + // A gzip stream delivered through a non-seekable FIFO must still decode. On the old code the + // gzip branch reopened the path, which on a pipe yields a fresh handle whose magic bytes were + // already consumed, so decoding failed with "invalid gzip header". + #[cfg(unix)] + #[test] + fn gzip_over_fifo_is_rewind_safe() { + use std::process::Command; + let raw = b">s1\nATGC\n>s2\nATCC\n"; + let gz = gzip(raw); + let fifo = std::env::temp_dir().join(format!("snpick-fifo-{}.gz", std::process::id())); + let _ = std::fs::remove_file(&fifo); + assert!(Command::new("mkfifo").arg(&fifo).status().unwrap().success()); + + let wpath = fifo.clone(); + let writer = std::thread::spawn(move || { + // File::create on a FIFO blocks until the reader opens; then stream and close. + if let Ok(mut f) = File::create(&wpath) { + let _ = f.write_all(&gz); + } + }); + let m = map_input(fifo.to_str().unwrap()).unwrap(); + assert_eq!(&m.mmap[..], raw); + writer.join().unwrap(); + std::fs::remove_file(&fifo).ok(); + } +} diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..ab108ae --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,43 @@ +//! # SNPick +//! +//! Fast, memory-efficient extraction of variable (SNP) sites from FASTA alignments. +//! +//! This library exposes the core of the `snpick` command-line tool as reusable modules: +//! +//! - [`fasta`] β€” zero-copy FASTA indexing over memory-mapped data. +//! - [`scan`] β€” pass 1: the parallel bitmask scan and site classification ([`scan::analyze`]). +//! - [`extract`] β€” pass 2: reduced FASTA / genotype-matrix extraction. +//! - [`vcf`] β€” VCF v4.2 writer. +//! - [`filter`] β€” per-site allele/missingness counting and filtering. +//! - [`coords`] β€” reference-anchored coordinates and BED masking. +//! - [`audit`] β€” alignment composition audit. +//! - [`input`] β€” transparent gzip/stdin input mapping. +//! - [`types`] β€” shared constants, lookup tables and data structures. +//! +//! # Example +//! +//! ```no_run +//! use snpick::fasta::index_fasta; +//! use snpick::scan::{analyze, pass1_scan}; +//! use snpick::types::build_lookup; +//! +//! let data: &[u8] = b">a\nACGT\n>b\nACGA\n"; +//! let lookup = build_lookup(false); +//! let (records, seq_len, layout) = index_fasta(data).unwrap(); +//! let bitmask = pass1_scan(data, &records, seq_len, layout, &lookup); +//! let ref_seq: Vec = data[3..7].to_vec(); +//! let (variable_sites, counts) = analyze(&bitmask, &ref_seq, &lookup, false); +//! assert_eq!(variable_sites.len(), 1); +//! assert_eq!(counts.constant.total(), 3); +//! ``` + +pub mod audit; +pub mod coords; +pub mod extract; +pub mod fasta; +pub mod filter; +#[cfg(feature = "mmap")] +pub mod input; +pub mod scan; +pub mod types; +pub mod vcf; diff --git a/src/main.rs b/src/main.rs index f1bb6cd..ce08bd0 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,39 +1,152 @@ -mod extract; -mod fasta; -mod scan; -mod types; -mod vcf; use clap::Parser; -use memmap2::Mmap; -use std::fs::File; -use std::io::{self, BufWriter, Write}; +use std::collections::HashSet; +use std::io::{self, BufWriter, IsTerminal, Write}; use std::path::Path; use std::time::Instant; -use crate::extract::{pass2_extract, ExtractParams}; -use crate::fasta::{get_ref_seq, index_fasta}; -use crate::scan::{analyze, pass1_scan}; -use crate::types::*; -use crate::vcf::write_vcf; +use snpick::extract::{pass2_extract, ExtractParams}; +use snpick::fasta::{get_ref_seq, index_fasta, FastaRecord}; +use snpick::filter::{count_sites, SiteFilters}; +use snpick::scan::{analyze, pass1_scan}; +use snpick::types::*; +use snpick::vcf::write_vcf; +use snpick::{audit, coords, input}; + +/// Emit a `[snpick]` progress line to stderr unless `--quiet` was passed. +/// Errors are always printed; only progress chatter is gated. +macro_rules! progress { + ($quiet:expr, $($arg:tt)*) => { + if !$quiet { eprintln!($($arg)*); } + }; +} // ============================================================================= // CLI // ============================================================================= +/// How to treat IUPAC ambiguity codes (R, Y, S, ...). +#[derive(clap::ValueEnum, Clone, Copy, Debug, PartialEq, Eq, Default)] +enum IupacMode { + /// Treat ambiguity codes as missing data (default). + #[default] + Missing, + /// Resolve ambiguity codes to their constituent bases when classifying sites. + Resolve, +} + +/// Policy for out-of-alphabet bytes in the input. +#[derive(clap::ValueEnum, Clone, Copy, Debug, PartialEq, Eq, Default)] +enum OnInvalid { + /// Skip the check entirely (no cost). + #[default] + Ignore, + /// Warn about out-of-alphabet bytes but continue. + Warn, + /// Error out if any out-of-alphabet byte is present. + Error, +} + #[derive(Parser, Debug)] #[command( name = "snpick", version = env!("CARGO_PKG_VERSION"), author = "Paula Ruiz-Rodriguez ", - about = "A fast, memory-efficient tool for extracting variable sites from FASTA alignments." + about = "Fast, memory-efficient extraction of variable sites from FASTA alignments.", + long_about = "snpick extracts variable (SNP) sites from a whole-genome FASTA alignment, \ +producing a reduced alignment ready for phylogenetic inference (optionally with a VCF and the \ +ASC fconst constant-site counts for IQ-TREE / RAxML). It uses a zero-copy, memory-mapped, \ +parallel two-pass scan that scales to thousands of genomes with minimal RAM.", + after_help = "NOTES:\n \ +- All input sequences must have the same length (an alignment). The first sequence is the\n \ +reference for REF/ALT polarity.\n \ +- fconst is printed as A,C,G,T (the order IQ-TREE's -fconst expects).\n \ +- In the VCF, POS is the 1-based ALIGNMENT column, not an ungapped reference coordinate.\n \ +- IUPAC ambiguous bases (N, R, Y, ...) are treated as missing data, never as alleles.\n \ +- Gaps ('-') are ignored unless -g is given, where they become a 5th allele ('*' in the VCF).\n\n\ +EXAMPLES:\n \ +snpick -f aln.fasta -o snps.fasta\n \ +snpick -f aln.fasta -o snps.fasta --vcf --chrom NC_000962.3\n \ +snpick -f aln.fasta -o snps.fasta -g -t 8 -q" )] struct Args { + /// Input FASTA alignment (all sequences must have equal length). #[arg(short, long)] fasta: String, - #[arg(short, long)] output: String, + /// Output FASTA containing only the variable sites (not needed with --dry-run/--check). + #[arg(short, long, required_unless_present_any = ["dry_run", "check"])] output: Option, + /// Treat gaps ('-') as a 5th character instead of ignoring them. #[arg(short = 'g', long)] include_gaps: bool, + /// Also write a VCF, named after the output (snps.fasta -> snps.vcf). #[arg(long)] vcf: bool, + /// Write the VCF to a custom path (implies --vcf). #[arg(long)] vcf_output: Option, + /// Silence progress logs on stderr (errors are still reported). + #[arg(short = 'q', long)] quiet: bool, + /// Increase detail: -v adds diagnostics, -vv adds an allele-class histogram. + #[arg(short, long, action = clap::ArgAction::Count)] verbose: u8, + /// Number of threads for the parallel scan (default: all logical cores). + #[arg(short = 't', long, value_parser = parse_threads)] + threads: Option, + /// CHROM / contig name written to the VCF (e.g. NC_000962.3). + #[arg(long, default_value = "1")] + chrom: String, + /// Sequence ID to use as the REF/polarity reference (default: first sequence). + #[arg(long)] reference: Option, + /// Permit duplicate sequence IDs instead of erroring. + #[arg(long)] allow_dup_ids: bool, + /// Write a machine-readable JSON run summary to this path ('-' for stdout). + #[arg(long)] stats_json: Option, + /// Report site statistics without writing any FASTA/VCF output. + #[arg(long)] dry_run: bool, + /// Keep only these sample IDs (comma-separated, or @file with one ID per line). + #[arg(long)] keep_samples: Option, + /// Drop these sample IDs (comma-separated, or @file). Excludes --keep-samples. + #[arg(long, conflicts_with = "keep_samples")] exclude_samples: Option, + /// Drop sites whose fraction of missing genotypes exceeds this (0.0-1.0). + #[arg(long, value_parser = parse_fraction)] max_missing: Option, + /// Drop sites whose minor-allele count is below this. + #[arg(long)] mac: Option, + /// Drop sites whose minor-allele frequency is below this (0.0-1.0). + #[arg(long, value_parser = parse_fraction)] maf: Option, + /// Drop sites with fewer than this many samples with data. + #[arg(long)] min_samples: Option, + /// Drop sites with more than this many distinct alleles (e.g. 2 = biallelic). + #[arg(long)] max_alleles: Option, + /// BED file of regions to mask out (excluded from output AND fconst). + #[arg(long)] mask: Option, + /// Interpret --mask coordinates as reference positions, not alignment columns. + #[arg(long, requires = "mask")] mask_ref: bool, + /// Use ungapped reference positions for VCF POS instead of the alignment column. + #[arg(long)] ref_coords: bool, + /// Write a TSV mapping each variable site (alignment_pos, ref_pos, ref, alt) to this path. + #[arg(long)] sites_output: Option, + /// What to do about out-of-alphabet bytes: ignore (default), warn, or error. + #[arg(long, value_enum, default_value_t = OnInvalid::Ignore)] on_invalid: OnInvalid, + /// Audit the alignment composition and exit without writing output. + #[arg(long)] check: bool, + /// How to treat IUPAC ambiguity codes: missing (default) or resolve to bases. + #[arg(long, value_enum, default_value_t = IupacMode::Missing)] iupac_mode: IupacMode, + /// Output format for the reduced alignment: fasta (default), phylip, or nexus. + #[arg(long, value_enum, default_value_t = OutputFormat::Fasta)] format: OutputFormat, +} + +/// Parse a positive thread count (Rayon treats 0 as "use default", which would +/// be a confusing silent no-op, so reject it explicitly). +fn parse_threads(s: &str) -> Result { + let n: usize = s.parse().map_err(|_| format!("'{}' is not a valid thread count", s))?; + if n == 0 { + return Err("thread count must be at least 1".to_string()); + } + Ok(n) +} + +/// Parse a fraction in [0.0, 1.0] (rejects NaN, infinities and out-of-range values). +fn parse_fraction(s: &str) -> Result { + let v: f64 = s.parse().map_err(|_| format!("'{}' is not a number", s))?; + if !v.is_finite() || !(0.0..=1.0).contains(&v) { + return Err(format!("'{}' must be a fraction between 0.0 and 1.0", s)); + } + Ok(v) } // ============================================================================= @@ -45,7 +158,12 @@ fn resolve_path(p: &str) -> io::Result { if path.exists() { return std::fs::canonicalize(path); } - let parent = path.parent().unwrap_or(Path::new(".")); + // A bare filename has an *empty* parent (Some("")), not None, so canonicalize + // would fail with ENOENT. Treat that as the current directory. + let parent = match path.parent() { + Some(p) if !p.as_os_str().is_empty() => p, + _ => Path::new("."), + }; let parent_abs = std::fs::canonicalize(parent).map_err(|e| { io::Error::new(e.kind(), format!("Cannot resolve parent of '{}': {}", p, e)) })?; @@ -53,6 +171,9 @@ fn resolve_path(p: &str) -> io::Result { } fn check_paths_differ(a: &str, b: &str) -> io::Result<()> { + if a == "-" || b == "-" { + return Ok(()); // stdin/stdout can't collide with a file + } let pa = resolve_path(a)?; let pb = resolve_path(b)?; if pa == pb { @@ -62,86 +183,384 @@ fn check_paths_differ(a: &str, b: &str) -> io::Result<()> { Ok(()) } +/// Reject empty sequence IDs, and duplicate IDs unless `allow_dups`. +/// Duplicate sample names yield spec-invalid VCFs and trees downstream tools reject. +fn validate_ids(records: &[FastaRecord], allow_dups: bool) -> io::Result<()> { + let mut seen: HashSet<&[u8]> = HashSet::with_capacity(records.len()); + for (i, rec) in records.iter().enumerate() { + if rec.id.is_empty() { + return Err(io::Error::new(io::ErrorKind::InvalidData, + format!("Record #{} has an empty sequence ID.", i + 1))); + } + if !allow_dups && !seen.insert(rec.id) { + let id = std::str::from_utf8(rec.id).unwrap_or("?"); + return Err(io::Error::new(io::ErrorKind::InvalidData, + format!("Duplicate sequence ID '{}' (record #{}). Use --allow-dup-ids to permit.", + id, i + 1))); + } + } + Ok(()) +} + +/// Emit a flat JSON run summary to `path` ('-' = stdout). Hand-written (no serde); +/// only the string fields need escaping. +#[allow(clippy::too_many_arguments)] +fn write_stats_json( + path: &str, input: &str, reference: &str, seq_length: usize, num_samples: usize, + sc: &SiteCounts, written: usize, include_gaps: bool, threads: usize, +) -> io::Result<()> { + let esc = |s: &str| s.replace('\\', "\\\\").replace('"', "\\\""); + let c = &sc.constant; + let json = format!( + "{{\n \"snpick_version\": \"{}\",\n \"input\": \"{}\",\n \"reference\": \"{}\",\n \ +\"sequences\": {},\n \"alignment_length\": {},\n \"variable_sites\": {},\n \ +\"written_sites\": {},\n \ +\"constant_sites\": {},\n \"constant_by_base\": {{ \"A\": {}, \"C\": {}, \"G\": {}, \"T\": {} }},\n \ +\"ambiguous_sites\": {},\n \"fconst\": [{}, {}, {}, {}],\n \ +\"include_gaps\": {},\n \"threads\": {}\n}}\n", + env!("CARGO_PKG_VERSION"), esc(input), esc(reference), + num_samples, seq_length, sc.variable, written, + c.total(), c.a, c.c, c.g, c.t, sc.ambiguous, c.a, c.c, c.g, c.t, + include_gaps, threads, + ); + if path == "-" { + io::stdout().write_all(json.as_bytes()) + } else { + std::fs::write(path, json) + } +} + +/// Parse a sample-ID selector: a comma-separated list, or `@path` to read one ID per line +/// (blank lines ignored). Every non-blank line is a literal ID β€” sample IDs come from FASTA +/// headers and may legitimately start with '#', so no line is treated as a comment (that would +/// silently drop such IDs and keep the wrong sample set, unlike the comma-separated form). +fn parse_id_set(spec: &str) -> io::Result>> { + let mut set = HashSet::new(); + if let Some(path) = spec.strip_prefix('@') { + let content = std::fs::read_to_string(path).map_err(|e| io::Error::new(e.kind(), + format!("Cannot read sample list '{}': {}", path, e)))?; + for line in content.lines() { + let line = line.trim(); + if !line.is_empty() { + set.insert(line.as_bytes().to_vec()); + } + } + } else { + for id in spec.split(',') { + let id = id.trim(); + if !id.is_empty() { set.insert(id.as_bytes().to_vec()); } + } + } + Ok(set) +} + +/// Apply --keep-samples / --exclude-samples in place. Filtering happens before the scan, so +/// site classification and fconst are computed for exactly the retained samples. Unknown IDs +/// are an error (catches typos). +fn apply_sample_filter( + records: &mut Vec, keep: &Option, exclude: &Option, +) -> io::Result<()> { + let (spec, keeping) = match (keep, exclude) { + (Some(s), _) => (s, true), + (_, Some(s)) => (s, false), + _ => return Ok(()), + }; + let set = parse_id_set(spec)?; + let present: HashSet<&[u8]> = records.iter().map(|r| r.id).collect(); + for id in &set { + if !present.contains(id.as_slice()) { + return Err(io::Error::new(io::ErrorKind::InvalidInput, + format!("Sample '{}' not found in the alignment.", String::from_utf8_lossy(id)))); + } + } + records.retain(|r| set.contains(r.id) == keeping); + if records.is_empty() { + return Err(io::Error::new(io::ErrorKind::InvalidInput, + "No sequences remain after sample filtering.")); + } + Ok(()) +} + +/// Write a TSV mapping each variable site to its alignment column, reference position, +/// REF and ALT alleles (gaps rendered as '*'). +fn write_sites_tsv(path: &str, var: &[VariablePosition], ref_pos: Option<&[u32]>) -> io::Result<()> { + let mut w = BufWriter::new(snpick::extract::open_sink(path)?); + writeln!(w, "alignment_pos\tref_pos\tref\talt")?; + for vp in var { + // Clamp to 1 to match the VCF POS (a leading-gap reference column has ungapped + // position 0, which is not a valid 1-based coordinate). + let rp = ref_pos.map(|r| r[vp.index].max(1)).unwrap_or((vp.index + 1) as u32); + let refc = if vp.ref_base == b'-' { '*' } else { vp.ref_base as char }; + let alt: String = vp.alt_bases.iter() + .map(|&b| if b == b'-' { "*".to_string() } else { (b as char).to_string() }) + .collect::>().join(","); + writeln!(w, "{}\t{}\t{}\t{}", vp.index + 1, rp, refc, alt)?; + } + w.flush() +} + +/// Index of the record to use as the REF/polarity reference. +fn reference_index(records: &[FastaRecord], reference: &Option) -> io::Result { + match reference { + Some(id) => { + let target = id.as_bytes(); + records.iter().position(|r| r.id == target).ok_or_else(|| { + io::Error::new(io::ErrorKind::InvalidInput, + format!("--reference '{}' not found among sequence IDs.", id)) + }) + } + None => Ok(0), + } +} + // ============================================================================= // Pipeline // ============================================================================= fn run() -> io::Result<()> { let args = Args::parse(); + let quiet = args.quiet; + + // Configure the Rayon pool up front. Absent, Rayon uses all logical cores; + // pinning it keeps wall-clock deterministic on shared HPC/SLURM nodes. + if let Some(n) = args.threads { + rayon::ThreadPoolBuilder::new().num_threads(n).build_global().map_err(|e| { + io::Error::other(format!("Cannot configure {} thread(s): {}", n, e)) + })?; + } + let start = Instant::now(); let lookup = build_lookup(args.include_gaps); let upper = build_upper(); - let do_vcf = args.vcf || args.vcf_output.is_some(); + let dry_run = args.dry_run; + // --dry-run reports statistics only; --check audits and exits. Neither writes the reduced + // alignment or VCF, and --check additionally writes no sidecars, so it needs no output paths. + let writes_align = !dry_run && !args.check; + let do_vcf = (args.vcf || args.vcf_output.is_some()) && writes_align; - // Validate paths - check_paths_differ(&args.fasta, &args.output)?; + // A CHROM with whitespace would break the tab-delimited VCF columns. + if do_vcf && (args.chrom.is_empty() || args.chrom.bytes().any(|b| b.is_ascii_whitespace())) { + return Err(io::Error::new(io::ErrorKind::InvalidInput, + "--chrom must be non-empty and contain no whitespace.")); + } + + // Output path (clap guarantees it is present unless --dry-run/--check). + let out_path: Option = if writes_align { args.output.clone() } else { None }; + + // Validate paths (skip when writing nothing β€” e.g. `--check -o input.fa` must not be + // rejected for a collision with a file it never touches). + if let Some(ref out) = out_path { + check_paths_differ(&args.fasta, out)?; + } let vcf_path = if do_vcf { - let vp = args.vcf_output.unwrap_or_else(|| { - let out = Path::new(&args.output); - let stem = out.file_stem().and_then(|s| s.to_str()).unwrap_or("output"); - let parent = out.parent().unwrap_or(Path::new(".")); - parent.join(format!("{}.vcf", stem)).to_string_lossy().into_owned() - }); + let out = out_path.as_deref().unwrap_or("output"); + let vp = match args.vcf_output.clone() { + Some(vp) => vp, + None => { + // Can't derive `.vcf` when the alignment streams to stdout (`-o -`); + // that used to fabricate a file literally named "-.vcf". + if out == "-" { + return Err(io::Error::new(io::ErrorKind::InvalidInput, + "Cannot derive a VCF path when the alignment is written to stdout (-o -); \ + pass --vcf-output (or --vcf-output - to stream the VCF to stdout).")); + } + let o = Path::new(out); + let stem = o.file_stem().and_then(|s| s.to_str()).unwrap_or("output"); + let parent = o.parent().unwrap_or(Path::new(".")); + parent.join(format!("{}.vcf", stem)).to_string_lossy().into_owned() + } + }; check_paths_differ(&args.fasta, &vp)?; - check_paths_differ(&args.output, &vp)?; + check_paths_differ(out, &vp)?; Some(vp) } else { None }; - // Memory-map input - let file = File::open(&args.fasta).map_err(|e| io::Error::new(e.kind(), - format!("Cannot open '{}': {}", args.fasta, e)))?; - let file_len = file.metadata()?.len(); - if file_len == 0 { + // Sidecars are emitted for normal and --dry-run runs, but not --check (which returns before + // they are written), so only validate their paths when they will actually be produced. + if !args.check { + // Sidecar outputs must not collide with the input (which is mmapped in place β€” + // overwriting it truncates the live mapping and corrupts the run) or the other outputs. + for sidecar in [args.stats_json.as_deref(), args.sites_output.as_deref()].into_iter().flatten() { + check_paths_differ(&args.fasta, sidecar)?; + if let Some(ref out) = out_path { + check_paths_differ(out, sidecar)?; + } + if let Some(ref vp) = vcf_path { + check_paths_differ(vp, sidecar)?; + } + } + // ...and the two sidecars must not clobber each other. + if let (Some(a), Some(b)) = (args.stats_json.as_deref(), args.sites_output.as_deref()) { + check_paths_differ(a, b)?; + } + // At most one output may stream to stdout; check_paths_differ exempts "-" (a file can't + // collide with stdout), so several "-" sinks would otherwise interleave β€” e.g. the JSON + // summary and the sites TSV concatenated onto one unparseable stream. + let to_stdout = [ + out_path.as_deref(), vcf_path.as_deref(), + args.stats_json.as_deref(), args.sites_output.as_deref(), + ].into_iter().flatten().filter(|p| *p == "-").count(); + if to_stdout > 1 { + return Err(io::Error::new(io::ErrorKind::InvalidInput, + "Multiple outputs are directed to stdout ('-'); at most one may use '-'.")); + } + } + + // Map the input, transparently decompressing gzip/bgzip and reading stdin ("-") as needed. + let input = input::map_input(&args.fasta)?; + if input.mmap.is_empty() { return Err(io::Error::new(io::ErrorKind::InvalidData, - format!("Input file '{}' is empty (0 bytes).", args.fasta))); + format!("Input '{}' is empty (0 bytes).", args.fasta))); } - let mmap = unsafe { Mmap::map(&file).map_err(|e| io::Error::new(e.kind(), - format!("Cannot memory-map '{}': {}", args.fasta, e)))? }; // Hint: pass 1 reads sequentially; OS can prefetch and release pages eagerly - mmap.advise(memmap2::Advice::Sequential).ok(); - let data = &mmap[..]; + input.mmap.advise(memmap2::Advice::Sequential).ok(); + let data = &input.mmap[..]; // Index records - let (records, seq_length, layout) = index_fasta(data)?; + let (mut records, seq_length, layout) = index_fasta(data)?; + + validate_ids(&records, args.allow_dup_ids)?; + apply_sample_filter(&mut records, &args.keep_samples, &args.exclude_samples)?; let num_samples = records.len(); - eprintln!("[snpick] Mapped {} bytes. {} sequences Γ— {} positions.{}", + let ref_idx = reference_index(&records, &args.reference)?; + let ref_name = std::str::from_utf8(records[ref_idx].id).unwrap_or("reference").to_string(); + + progress!(quiet, "[snpick] Mapped {} bytes. {} sequences Γ— {} positions.{}", data.len(), num_samples, seq_length, if layout.single_line { "" } else { " (multi-line FASTA)" }); - // Pass 1: bitmask scan - let bitmask = pass1_scan(data, &records, seq_length, layout, &lookup); - let ref_seq = get_ref_seq(data, &records[0], seq_length, layout); + // Composition audit for --on-invalid / --check (a full extra pass, only on demand). + if args.check || args.on_invalid != OnInvalid::Ignore { + let hist = audit::composition(data, &records, seq_length, layout); + let invalid = audit::invalid_count(&hist); + if args.on_invalid == OnInvalid::Error && invalid > 0 { + return Err(io::Error::new(io::ErrorKind::InvalidData, format!( + "{} out-of-alphabet byte(s) found β€” is this a nucleotide alignment? \ + (use --on-invalid ignore to allow).", invalid))); + } + if args.on_invalid == OnInvalid::Warn && invalid > 0 { + eprintln!("[snpick] Warning: {} out-of-alphabet byte(s) in the alignment.", invalid); + } + if args.check { + progress!(quiet, "[snpick] Composition: {}", audit::summary(&hist)); + return Ok(()); + } + } + + // Pass 1: bitmask scan. Under --iupac-mode resolve, ambiguity codes are expanded to their + // bases *only here* (for classification); NS counting and REF selection keep the strict table. + let scan_lookup = if args.iupac_mode == IupacMode::Resolve { + build_iupac_lookup(args.include_gaps) + } else { + lookup + }; + let mut bitmask = pass1_scan(data, &records, seq_length, layout, &scan_lookup); + let ref_seq = get_ref_seq(data, &records[ref_idx], seq_length, layout); let t1 = start.elapsed().as_secs_f64(); + // Ungapped reference positions (only computed when a feature needs them). + let ref_pos: Option> = + if args.ref_coords || args.sites_output.is_some() || args.mask_ref { + Some(coords::ref_positions(&ref_seq)) + } else { + None + }; + + // Mask out BED regions before classification: masked columns are zeroed, so they classify + // as ambiguous and never enter the output or the fconst constant counts. + if let Some(bed) = &args.mask { + let ivs = coords::parse_bed(bed)?; + let rc = if args.mask_ref { ref_pos.as_deref() } else { None }; + let mask = coords::build_mask(&ivs, seq_length, rc); + let masked = mask.iter().filter(|&&m| m).count(); + for (col, &m) in mask.iter().enumerate() { + if m { bitmask[col] = 0; } + } + progress!(quiet, "[snpick] Masked {} columns from {} BED region(s).", masked, ivs.len()); + } + let (mut var_positions, site_counts) = analyze(&bitmask, &ref_seq, &lookup, args.include_gaps); - let num_var = var_positions.len(); drop(bitmask); drop(ref_seq); - eprintln!("[snpick] {} variable, {} constant ({}), {} ambiguous-only, {} total.", + progress!(quiet, "[snpick] {} variable, {} constant ({}), {} ambiguous-only, {} total.", site_counts.variable, site_counts.constant.total(), site_counts.constant, site_counts.ambiguous, seq_length); - eprintln!("[snpick] ASC fconst: {}", site_counts.constant.fconst()); - eprintln!("[snpick] Pass 1 took {:.2}s.", t1); + progress!(quiet, "[snpick] ASC fconst: {}", site_counts.constant.fconst()); + progress!(quiet, "[snpick] Pass 1 took {:.2}s.", t1); - // Handle zero-variant case - if num_var == 0 { - eprintln!("[snpick] No variable positions β€” writing empty output."); - let out = File::create(&args.output)?; - let mut w = BufWriter::new(out); - for rec in &records { - w.write_all(b">")?; - w.write_all(rec.id)?; - if !rec.desc.is_empty() { w.write_all(b" ")?; w.write_all(rec.desc)?; } - writeln!(w)?; writeln!(w)?; + // Per-site filtering (opt-in). Removes variable sites from the OUTPUT only; the fconst + // constant-site counts are unchanged, so ASC stays valid. + let filters = SiteFilters { + max_missing: args.max_missing, mac: args.mac, maf: args.maf, + min_samples: args.min_samples, max_alleles: args.max_alleles, + }; + if filters.active() && !var_positions.is_empty() { + let pos_indices: Vec = var_positions.iter().map(|v| v.index).collect(); + // Count over the SAME table pass 1 used to classify (resolves IUPAC codes under + // --iupac-mode resolve), so the filter judges the same allele set. + let stats = count_sites(data, &records, &pos_indices, layout, &scan_lookup); + let ns_total = num_samples as u32; + let before = var_positions.len(); + let mut i = 0; + var_positions.retain(|_| { + let keep = stats[i].passes(&filters, ns_total, args.include_gaps); + i += 1; + keep + }); + progress!(quiet, "[snpick] Filtered {} of {} variable sites ({} remain).", + before - var_positions.len(), before, var_positions.len()); + } + let num_var = var_positions.len(); + + if args.verbose >= 1 && !quiet { + let mbps = data.len() as f64 / 1e6 / t1.max(1e-9); + eprintln!("[snpick] Threads: {}. Layout: {}. Scan throughput: {:.0} MB/s.", + rayon::current_num_threads(), + if layout.single_line { "single-line" } else { "multi-line" }, mbps); + } + if args.verbose >= 2 && !quiet { + let (mut bi, mut tri, mut tetra) = (0usize, 0usize, 0usize); + for vp in &var_positions { + match 1 + vp.alt_bases.len() { + 2 => bi += 1, + 3 => tri += 1, + _ => tetra += 1, + } } - w.flush()?; + eprintln!("[snpick] Allele classes: biallelic={} triallelic={} 4+-allelic={}", bi, tri, tetra); + } + + // Machine-readable stats sidecar (also emitted for dry-run and zero-variant). + if let Some(ref sp) = args.stats_json { + write_stats_json(sp, &args.fasta, &ref_name, seq_length, num_samples, + &site_counts, num_var, args.include_gaps, rayon::current_num_threads())?; + progress!(quiet, "[snpick] Stats JSON written to {}.", + if sp == "-" { "stdout" } else { sp }); + } + + // Variable-site coordinate map β€” a reporting sidecar, so emitted for --dry-run too. + if let Some(sp) = &args.sites_output { + write_sites_tsv(sp, &var_positions, ref_pos.as_deref())?; + progress!(quiet, "[snpick] Sites TSV written to {}.", sp); + } + + if dry_run { + progress!(quiet, "[snpick] Dry run β€” no alignment/VCF written."); return Ok(()); } + // Past the dry-run gate, an output path is guaranteed (clap-enforced). + let out = out_path.as_deref().expect("output is required unless --dry-run"); + let pos_map = if args.ref_coords { ref_pos.as_deref() } else { None }; + + if num_var == 0 { + progress!(quiet, "[snpick] No variable positions found β€” writing empty alignment."); + } + // VCF size guard if do_vcf { let geno_bytes = num_var.saturating_mul(num_samples); @@ -153,20 +572,23 @@ fn run() -> io::Result<()> { } } - // Pass 2: extract variable sites + // Pass 2: extract variable sites (show an in-place progress line on a TTY for big inputs). + let show_progress = !quiet && num_samples >= 500 && io::stderr().is_terminal() && out != "-"; let ep = ExtractParams { - records: &records, output: &args.output, - collect_vcf: do_vcf, lookup: &lookup, upper: &upper, layout, + records: &records, output: out, + collect_vcf: do_vcf, lookup: &lookup, upper: &upper, layout, format: args.format, + progress: show_progress, }; let vcf_geno = pass2_extract(data, &mut var_positions, &ep)?; + progress!(quiet, "[snpick] Pass 2: Wrote {} sequences to {}.", num_samples, out); // Write VCF if let (Some(ref geno), Some(ref vp)) = (&vcf_geno, &vcf_path) { - write_vcf(geno, num_samples, &var_positions, vp, &records, seq_length)?; - eprintln!("[snpick] VCF written to {}.", vp); + write_vcf(geno, num_samples, &var_positions, vp, &records, seq_length, &args.chrom, &ref_name, pos_map)?; + progress!(quiet, "[snpick] VCF written to {}.", if vp == "-" { "stdout" } else { vp }); } - eprintln!("[snpick] Done in {:.2}s. {} vars from {} seqs Γ— {} pos.", + progress!(quiet, "[snpick] Done in {:.2}s. {} vars from {} seqs Γ— {} pos.", start.elapsed().as_secs_f64(), num_var, num_samples, seq_length); Ok(()) } @@ -174,7 +596,12 @@ fn run() -> io::Result<()> { fn main() { if let Err(e) = run() { eprintln!("[snpick] Error: {}", e); - std::process::exit(1); + // Exit codes: 2 for bad input/data, 1 for I/O and everything else. + let code = match e.kind() { + io::ErrorKind::InvalidData | io::ErrorKind::InvalidInput => 2, + _ => 1, + }; + std::process::exit(code); } } @@ -185,10 +612,12 @@ fn main() { #[cfg(test)] mod tests { use super::*; - use crate::extract::ExtractParams; - use crate::fasta::{get_ref_seq, index_fasta}; - use crate::scan::{analyze, pass1_scan}; - use crate::vcf::write_vcf; + use memmap2::Mmap; + use std::fs::File; + use snpick::extract::ExtractParams; + use snpick::fasta::{get_ref_seq, index_fasta}; + use snpick::scan::{analyze, pass1_scan}; + use snpick::vcf::write_vcf; fn tmp(name: &str, c: &str) -> String { let p = format!("/tmp/snpick_t_{}.fa", name); @@ -206,6 +635,19 @@ mod tests { assert_eq!(build_lookup(true)[b'-' as usize], BIT_GAP); } + #[test] fn parse_id_set_at_file_keeps_hash_ids() { + // @file must treat every non-blank line as a literal ID β€” including IDs that start + // with '#' (valid FASTA headers) β€” and agree with the comma-separated form. + let p = format!("/tmp/snpick_idset_{}.txt", std::process::id()); + std::fs::write(&p, "#iso1\n\n normal \n").unwrap(); + let from_file = parse_id_set(&format!("@{}", p)).unwrap(); + assert!(from_file.contains(&b"#iso1"[..])); + assert!(from_file.contains(&b"normal"[..])); + assert_eq!(from_file.len(), 2); + assert_eq!(from_file, parse_id_set("#iso1,normal").unwrap()); + std::fs::remove_file(&p).ok(); + } + #[test] fn test_index() { let p = tmp("idx3", ">s1\nATGC\n>s2\nATCC\n"); let m = setup(&p); @@ -238,7 +680,7 @@ mod tests { let bm = pass1_scan(&m, &recs, sl, layout, &lk); let rs = get_ref_seq(&m, &recs[0], sl, layout); let (mut v, _) = analyze(&bm, &rs, &lk, false); - let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; pass2_extract(&m, &mut v, &ep).unwrap(); let c = std::fs::read_to_string(o).unwrap(); let l: Vec<&str> = c.lines().collect(); @@ -256,9 +698,9 @@ mod tests { let bm = pass1_scan(&m, &recs, sl, layout, &lk); let rs = get_ref_seq(&m, &recs[0], sl, layout); let (mut v, _) = analyze(&bm, &rs, &lk, false); - let ep = ExtractParams { records: &recs, output: fo, collect_vcf: true, lookup: &lk, upper: &up, layout }; + let ep = ExtractParams { records: &recs, output: fo, collect_vcf: true, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; let g = pass2_extract(&m, &mut v, &ep).unwrap().unwrap(); - write_vcf(&g, recs.len(), &v, vo, &recs, sl).unwrap(); + write_vcf(&g, recs.len(), &v, vo, &recs, sl, "1", "ref", None).unwrap(); let c = std::fs::read_to_string(vo).unwrap(); let dl: Vec<&str> = c.lines().filter(|l| !l.starts_with('#')).collect(); let f: Vec<&str> = dl[0].split('\t').collect(); @@ -277,9 +719,9 @@ mod tests { let bm = pass1_scan(&m, &recs, sl, layout, &lk); let rs = get_ref_seq(&m, &recs[0], sl, layout); let (mut v, _) = analyze(&bm, &rs, &lk, false); - let ep = ExtractParams { records: &recs, output: fo, collect_vcf: true, lookup: &lk, upper: &up, layout }; + let ep = ExtractParams { records: &recs, output: fo, collect_vcf: true, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; let g = pass2_extract(&m, &mut v, &ep).unwrap().unwrap(); - write_vcf(&g, recs.len(), &v, vo, &recs, sl).unwrap(); + write_vcf(&g, recs.len(), &v, vo, &recs, sl, "1", "ref", None).unwrap(); let c = std::fs::read_to_string(vo).unwrap(); let dl: Vec<&str> = c.lines().filter(|l| !l.starts_with('#')).collect(); let f: Vec<&str> = dl[0].split('\t').collect(); @@ -287,6 +729,97 @@ mod tests { std::fs::remove_file(&p).ok(); std::fs::remove_file(fo).ok(); std::fs::remove_file(vo).ok(); } + #[test] fn test_vcf_gap_ref() { + // Gap in the reference at a variable site: REF must be '*' (valid VCF v4.2), + // never a bare '-'. + let p = tmp("vgref", ">ref\nA-GC\n>s1\nATGC\n>s2\nA-GT\n"); + let fo = "/tmp/snpick_t_vgref_out.fa"; let vo = "/tmp/snpick_t_vgref.vcf"; + let m = setup(&p); + let lk = build_lookup(true); + let up = build_upper(); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + let bm = pass1_scan(&m, &recs, sl, layout, &lk); + let rs = get_ref_seq(&m, &recs[0], sl, layout); + let (mut v, _) = analyze(&bm, &rs, &lk, true); + let ep = ExtractParams { records: &recs, output: fo, collect_vcf: true, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; + let g = pass2_extract(&m, &mut v, &ep).unwrap().unwrap(); + write_vcf(&g, recs.len(), &v, vo, &recs, sl, "1", "ref", None).unwrap(); + let c = std::fs::read_to_string(vo).unwrap(); + let dl: Vec<&str> = c.lines().filter(|l| !l.starts_with('#')).collect(); + let f: Vec<&str> = dl[0].split('\t').collect(); + assert_eq!(f[1], "2"); // POS (gap-in-ref site) + assert_eq!(f[3], "N"); // REF: gap β†’ 'N' (valid VCF), not '-' or '*' + assert_eq!(f[4], "T"); // ALT + assert!(!c.contains("\t-\t")); // no bare hyphen field anywhere + std::fs::remove_file(&p).ok(); std::fs::remove_file(fo).ok(); std::fs::remove_file(vo).ok(); + } + + #[test] fn test_vcf_gap_alt() { + // Gap sample at a multi-allelic site renders as the '*' ALT allele. + let p = tmp("vgalt", ">ref\nATGC\n>s1\nTTGC\n>s2\n-TGC\n"); + let fo = "/tmp/snpick_t_vgalt_out.fa"; let vo = "/tmp/snpick_t_vgalt.vcf"; + let m = setup(&p); + let lk = build_lookup(true); + let up = build_upper(); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + let bm = pass1_scan(&m, &recs, sl, layout, &lk); + let rs = get_ref_seq(&m, &recs[0], sl, layout); + let (mut v, _) = analyze(&bm, &rs, &lk, true); + let ep = ExtractParams { records: &recs, output: fo, collect_vcf: true, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; + let g = pass2_extract(&m, &mut v, &ep).unwrap().unwrap(); + write_vcf(&g, recs.len(), &v, vo, &recs, sl, "1", "ref", None).unwrap(); + let c = std::fs::read_to_string(vo).unwrap(); + let dl: Vec<&str> = c.lines().filter(|l| !l.starts_with('#')).collect(); + let f: Vec<&str> = dl[0].split('\t').collect(); + assert_eq!(f[1], "1"); // POS + assert_eq!(f[3], "A"); // REF + assert_eq!(f[4], "T,*"); // ALT: gap sample β†’ '*' + std::fs::remove_file(&p).ok(); std::fs::remove_file(fo).ok(); std::fs::remove_file(vo).ok(); + } + + #[test] fn test_vcf_allgap_ref_contig_len() { + // All-gap reference under --ref-coords: ungapped length is 0, but POS clamps to 1, so the + // declared contig length must also be >= 1 β€” an emitted POS must never exceed it. + let p = tmp("agref", ">ref\n--\n>s1\nAC\n>s2\nGT\n"); + let fo = "/tmp/snpick_t_agref_out.fa"; let vo = "/tmp/snpick_t_agref.vcf"; + let m = setup(&p); + let lk = build_lookup(false); + let up = build_upper(); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + let bm = pass1_scan(&m, &recs, sl, layout, &lk); + let rs = get_ref_seq(&m, &recs[0], sl, layout); + let (mut v, _) = analyze(&bm, &rs, &lk, false); + let ep = ExtractParams { records: &recs, output: fo, collect_vcf: true, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; + let g = pass2_extract(&m, &mut v, &ep).unwrap().unwrap(); + let rp = snpick::coords::ref_positions(&rs); + write_vcf(&g, recs.len(), &v, vo, &recs, sl, "1", "ref", Some(&rp)).unwrap(); + let c = std::fs::read_to_string(vo).unwrap(); + let contig_len: usize = c.lines() + .find_map(|l| l.strip_prefix("##contig=').parse().ok()).unwrap(); + assert!(contig_len >= 1, "contig length was {}", contig_len); + for l in c.lines().filter(|l| !l.starts_with('#')) { + let pos: usize = l.split('\t').nth(1).unwrap().parse().unwrap(); + assert!(pos <= contig_len, "POS {} exceeds contig length {}", pos, contig_len); + } + std::fs::remove_file(&p).ok(); std::fs::remove_file(fo).ok(); std::fs::remove_file(vo).ok(); + } + + #[test] fn test_vcf_header_only() { + // Zero variable sites but VCF requested: a valid header-only VCF that + // still lists every sample, so downstream pipelines find the file. + let p = tmp("hdr", ">s1\nATGC\n>s2\nATGC\n"); + let vo = "/tmp/snpick_t_hdr.vcf"; + let m = setup(&p); + let (recs, sl, _layout) = index_fasta(&m).unwrap(); + write_vcf(&[], recs.len(), &[], vo, &recs, sl, "1", "ref", None).unwrap(); + let c = std::fs::read_to_string(vo).unwrap(); + let data: Vec<&str> = c.lines().filter(|l| !l.starts_with('#')).collect(); + assert_eq!(data.len(), 0); + assert!(c.contains("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\ts1\ts2")); + std::fs::remove_file(&p).ok(); std::fs::remove_file(vo).ok(); + } + #[test] fn test_desc_preserved() { let p = tmp("descg", ">s1 some description\nATGC\n>s2 another desc\nATCC\n"); let o = "/tmp/snpick_t_descg_out.fa"; @@ -297,7 +830,7 @@ mod tests { let bm = pass1_scan(&m, &recs, sl, layout, &lk); let rs = get_ref_seq(&m, &recs[0], sl, layout); let (mut v, _) = analyze(&bm, &rs, &lk, false); - let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; pass2_extract(&m, &mut v, &ep).unwrap(); let c = std::fs::read_to_string(o).unwrap(); assert!(c.contains(">s1 some description")); @@ -321,10 +854,157 @@ mod tests { std::fs::remove_file(&p).ok(); } + #[test] fn test_validate_ids() { + let p = tmp("ids", ">a\nAC\n>a\nAT\n>b\nAG\n"); + let m = setup(&p); + let (recs, _, _) = index_fasta(&m).unwrap(); + assert!(validate_ids(&recs, false).is_err()); // duplicate 'a' + assert!(validate_ids(&recs, true).is_ok()); // allowed + std::fs::remove_file(&p).ok(); + } + + #[test] fn test_stats_json() { + let sc = SiteCounts { + constant: ConstantSiteCounts { a: 7, c: 7, g: 7, t: 4 }, + variable: 5, + ambiguous: 1, + }; + let p = "/tmp/snpick_t_stats.json"; + write_stats_json(p, "aln.fasta", "H37Rv", 30, 6, &sc, 3, false, 8).unwrap(); + let c = std::fs::read_to_string(p).unwrap(); + assert!(c.contains("\"variable_sites\": 5")); + assert!(c.contains("\"written_sites\": 3")); + assert!(c.contains("\"fconst\": [7, 7, 7, 4]")); + assert!(c.contains("\"reference\": \"H37Rv\"")); + assert!(c.contains("\"sequences\": 6")); + std::fs::remove_file(p).ok(); + } + + #[test] fn test_phylip_output() { + let p = tmp("phy", ">s1\nATGC\n>s2\nATGT\n"); + let o = "/tmp/snpick_t_phy_out.phy"; + let m = setup(&p); + let lk = build_lookup(false); + let up = build_upper(); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + let bm = pass1_scan(&m, &recs, sl, layout, &lk); + let rs = get_ref_seq(&m, &recs[0], sl, layout); + let (mut v, _) = analyze(&bm, &rs, &lk, false); + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout, format: OutputFormat::Phylip, progress: false }; + pass2_extract(&m, &mut v, &ep).unwrap(); + let l: Vec = std::fs::read_to_string(o).unwrap().lines().map(|s| s.to_string()).collect(); + assert_eq!(l[0], "2 1"); + assert_eq!(l[1], "s1 C"); + assert_eq!(l[2], "s2 T"); + std::fs::remove_file(&p).ok(); std::fs::remove_file(o).ok(); + } + + #[test] fn test_count_sites() { + // pos0 singleton (A5 T1), pos1 common (A3 C3), pos2 triallelic (A2 C2 G2). + let p = tmp("cnt", ">s1\nAAA\n>s2\nAAA\n>s3\nAAC\n>s4\nACC\n>s5\nACG\n>s6\nTCG\n"); + let m = setup(&p); + let lk = build_lookup(false); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + let bm = pass1_scan(&m, &recs, sl, layout, &lk); + let rs = get_ref_seq(&m, &recs[0], sl, layout); + let (v, _) = analyze(&bm, &rs, &lk, false); + let idx: Vec = v.iter().map(|x| x.index).collect(); + let stats = snpick::filter::count_sites(&m, &recs, &idx, layout, &lk); + assert_eq!(stats.len(), 3); + assert_eq!(stats[0].counts, [5, 0, 0, 1]); + assert_eq!(stats[1].counts, [3, 3, 0, 0]); + assert_eq!(stats[2].counts, [2, 2, 2, 0]); + std::fs::remove_file(&p).ok(); + } + + #[test] fn test_exclude_samples_reclassifies() { + // A site variable only because of one sample becomes constant (and enters + // fconst) once that sample is excluded β€” the correctness snp-sites misses. + let p = tmp("excl", ">a\nAA\n>b\nAA\n>c\nAT\n"); + let m = setup(&p); + let lk = build_lookup(false); + let (mut recs, sl, layout) = index_fasta(&m).unwrap(); + let bm = pass1_scan(&m, &recs, sl, layout, &lk); + let rs = get_ref_seq(&m, &recs[0], sl, layout); + let (v, sc) = analyze(&bm, &rs, &lk, false); + assert_eq!(v.len(), 1); + assert_eq!(sc.constant.total(), 1); + apply_sample_filter(&mut recs, &None, &Some("c".to_string())).unwrap(); + assert_eq!(recs.len(), 2); + let bm2 = pass1_scan(&m, &recs, sl, layout, &lk); + let rs2 = get_ref_seq(&m, &recs[0], sl, layout); + let (v2, sc2) = analyze(&bm2, &rs2, &lk, false); + assert_eq!(v2.len(), 0); + assert_eq!(sc2.constant.total(), 2); + std::fs::remove_file(&p).ok(); + } + + #[test] fn test_keep_samples_and_unknown() { + let p = tmp("keep", ">a\nAT\n>b\nCT\n>c\nGT\n"); + let m = setup(&p); + let (mut recs, _, _) = index_fasta(&m).unwrap(); + apply_sample_filter(&mut recs, &Some("a,c".to_string()), &None).unwrap(); + assert_eq!(recs.len(), 2); + assert_eq!(recs[0].id, b"a"); + assert_eq!(recs[1].id, b"c"); + let (mut recs2, _, _) = index_fasta(&m).unwrap(); + assert!(apply_sample_filter(&mut recs2, &Some("a,zzz".to_string()), &None).is_err()); + std::fs::remove_file(&p).ok(); + } + + #[test] fn test_reference_index() { + let p = tmp("refi", ">a\nAC\n>b\nAT\n>c\nAG\n"); + let m = setup(&p); + let (recs, _, _) = index_fasta(&m).unwrap(); + assert_eq!(reference_index(&recs, &None).unwrap(), 0); + assert_eq!(reference_index(&recs, &Some("b".to_string())).unwrap(), 1); + assert!(reference_index(&recs, &Some("zzz".to_string())).is_err()); + std::fs::remove_file(&p).ok(); + } + + #[test] fn test_reference_polarity() { + // Choosing a different reference flips which base is REF vs ALT. + let p = tmp("refp", ">a\nAT\n>b\nCT\n"); + let m = setup(&p); + let lk = build_lookup(false); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + let bm = pass1_scan(&m, &recs, sl, layout, &lk); + // default reference = record 0 ('a') -> REF A + let rs0 = get_ref_seq(&m, &recs[0], sl, layout); + let (v0, _) = analyze(&bm, &rs0, &lk, false); + assert_eq!(v0[0].ref_base, b'A'); + assert_eq!(v0[0].alt_bases, vec![b'C']); + // reference = record 1 ('b') -> REF C + let rs1 = get_ref_seq(&m, &recs[1], sl, layout); + let (v1, _) = analyze(&bm, &rs1, &lk, false); + assert_eq!(v1[0].ref_base, b'C'); + assert_eq!(v1[0].alt_bases, vec![b'A']); + std::fs::remove_file(&p).ok(); + } + + #[test] fn test_all_gap_column_counted() { + // Under --include-gaps an all-gap column is tallied as ambiguous, so + // variable + constant + ambiguous still equals seq_length. + let p = tmp("allgap", ">s1\nA-\n>s2\nA-\n"); + let m = setup(&p); + let lk = build_lookup(true); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + let bm = pass1_scan(&m, &recs, sl, layout, &lk); + let rs = get_ref_seq(&m, &recs[0], sl, layout); + let (v, sc) = analyze(&bm, &rs, &lk, true); + assert_eq!(v.len(), 0); + assert_eq!(sc.constant.total(), 1); + assert_eq!(sc.ambiguous, 1); + assert_eq!(sc.variable + sc.constant.total() + sc.ambiguous, sl); + std::fs::remove_file(&p).ok(); + } + #[test] fn test_paths() { let p = tmp("pdg", ">x\nA\n"); assert!(check_paths_differ(&p, "/tmp/snpick_t_pdg2.fa").is_ok()); assert!(check_paths_differ(&p, &p).is_err()); + // A bare (directory-less) output name must resolve to the cwd, not error. + assert!(check_paths_differ(&p, "snpick_t_pdg_bare_out.fa").is_ok()); std::fs::remove_file(&p).ok(); } @@ -353,7 +1033,7 @@ mod tests { let bm = pass1_scan(&m, &recs, sl, layout, &lk); let rs = get_ref_seq(&m, &recs[0], sl, layout); let (mut v, _) = analyze(&bm, &rs, &lk, false); - let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; pass2_extract(&m, &mut v, &ep).unwrap(); let c = std::fs::read_to_string(o).unwrap(); let l: Vec<&str> = c.lines().collect(); @@ -371,9 +1051,9 @@ mod tests { let bm = pass1_scan(&m, &recs, sl, layout, &lk); let rs = get_ref_seq(&m, &recs[0], sl, layout); let (mut v, _) = analyze(&bm, &rs, &lk, false); - let ep = ExtractParams { records: &recs, output: fo, collect_vcf: true, lookup: &lk, upper: &up, layout }; + let ep = ExtractParams { records: &recs, output: fo, collect_vcf: true, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; let g = pass2_extract(&m, &mut v, &ep).unwrap().unwrap(); - write_vcf(&g, recs.len(), &v, vo, &recs, sl).unwrap(); + write_vcf(&g, recs.len(), &v, vo, &recs, sl, "1", "ref", None).unwrap(); let c = std::fs::read_to_string(fo).unwrap(); let l: Vec<&str> = c.lines().collect(); assert_eq!(l[1], "AG"); assert_eq!(l[3], "AC"); assert_eq!(l[5], "CG"); @@ -395,7 +1075,7 @@ mod tests { let rs = get_ref_seq(&m, &recs[0], sl, layout); let (mut v, _) = analyze(&bm, &rs, &lk, false); assert_eq!(v.len(), 1); assert_eq!(v[0].index, 4); - let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; pass2_extract(&m, &mut v, &ep).unwrap(); let c = std::fs::read_to_string(o).unwrap(); let l: Vec<&str> = c.lines().collect(); @@ -415,7 +1095,7 @@ mod tests { let rs = get_ref_seq(&m, &recs[0], sl, layout); let (mut v, _) = analyze(&bm, &rs, &lk, false); assert_eq!(v.len(), 1); assert_eq!(v[0].index, 6); - let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; pass2_extract(&m, &mut v, &ep).unwrap(); let c = std::fs::read_to_string(o).unwrap(); let l: Vec<&str> = c.lines().collect(); @@ -439,7 +1119,7 @@ mod tests { assert_eq!(v.len(), 1); assert_eq!(v[0].index, 2); let o = "/tmp/snpick_t_crlfml_out.fa"; - let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; pass2_extract(&m, &mut v, &ep).unwrap(); let c = std::fs::read_to_string(o).unwrap(); let l: Vec<&str> = c.lines().collect(); @@ -447,6 +1127,30 @@ mod tests { std::fs::remove_file(&p).ok(); std::fs::remove_file(o).ok(); } + #[test] fn test_single_line_blank_line() { + // Blank line after the header in an otherwise single-line file: must not + // drop the SNP. The record is forced onto the newline-skipping path. + let p = tmp("blank", ">s1\n\nAAAA\n>s2\n\nAAAT\n"); + let o = "/tmp/snpick_t_blank_out.fa"; + let m = setup(&p); + let lk = build_lookup(false); + let up = build_upper(); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + assert_eq!(sl, 4); + assert!(!layout.single_line); + let bm = pass1_scan(&m, &recs, sl, layout, &lk); + let rs = get_ref_seq(&m, &recs[0], sl, layout); + let (mut v, _) = analyze(&bm, &rs, &lk, false); + assert_eq!(v.len(), 1); + assert_eq!(v[0].index, 3); + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; + pass2_extract(&m, &mut v, &ep).unwrap(); + let c = std::fs::read_to_string(o).unwrap(); + let l: Vec<&str> = c.lines().collect(); + assert_eq!(l[1], "A"); assert_eq!(l[3], "T"); + std::fs::remove_file(&p).ok(); std::fs::remove_file(o).ok(); + } + #[test] fn test_single_sequence() { // Single sequence should produce 0 variable sites let p = tmp("sing", ">s1\nATGC\n"); @@ -476,7 +1180,7 @@ mod tests { let (mut v, _) = analyze(&bm, &rs, &lk, false); assert_eq!(v.len(), 1); let o = "/tmp/snpick_t_noeof_out.fa"; - let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; pass2_extract(&m, &mut v, &ep).unwrap(); let c = std::fs::read_to_string(o).unwrap(); let l: Vec<&str> = c.lines().collect(); @@ -484,6 +1188,83 @@ mod tests { std::fs::remove_file(&p).ok(); std::fs::remove_file(o).ok(); } + #[test] fn test_ambiguous_ref_base() { + // Reference (first sequence) is N at a variable site: REF falls back to + // the first observed base in A,C,G,T order and the ref genotype is '.'. + let p = tmp("nref", ">ref\nNTGC\n>s1\nATGC\n>s2\nCTGC\n"); + let fo = "/tmp/snpick_t_nref_out.fa"; let vo = "/tmp/snpick_t_nref.vcf"; + let m = setup(&p); + let lk = build_lookup(false); + let up = build_upper(); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + let bm = pass1_scan(&m, &recs, sl, layout, &lk); + let rs = get_ref_seq(&m, &recs[0], sl, layout); + let (mut v, _) = analyze(&bm, &rs, &lk, false); + assert_eq!(v.len(), 1); + assert_eq!(v[0].ref_base, b'A'); + assert_eq!(v[0].alt_bases, vec![b'C']); + let ep = ExtractParams { records: &recs, output: fo, collect_vcf: true, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; + let g = pass2_extract(&m, &mut v, &ep).unwrap().unwrap(); + write_vcf(&g, recs.len(), &v, vo, &recs, sl, "1", "ref", None).unwrap(); + let c = std::fs::read_to_string(vo).unwrap(); + let dl: Vec<&str> = c.lines().filter(|l| !l.starts_with('#')).collect(); + let f: Vec<&str> = dl[0].split('\t').collect(); + assert_eq!(f[3], "A"); // REF fallback + assert_eq!(f[4], "C"); // ALT + assert_eq!(f[7], "NS=2"); // ref N excluded + assert_eq!(f[9], "."); // ref sample genotype missing + assert_eq!(f[10], "0"); // s1 = A = ref + assert_eq!(f[11], "1"); // s2 = C = alt + std::fs::remove_file(&p).ok(); std::fs::remove_file(fo).ok(); std::fs::remove_file(vo).ok(); + } + + #[test] fn test_lowercase_normalization() { + // Soft-masked (lowercase) input classifies case-insensitively and is + // written uppercase in the reduced FASTA. + let p = tmp("lower", ">ref\natgc\n>s1\natgt\n"); + let o = "/tmp/snpick_t_lower_out.fa"; + let m = setup(&p); + let lk = build_lookup(false); + let up = build_upper(); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + let bm = pass1_scan(&m, &recs, sl, layout, &lk); + let rs = get_ref_seq(&m, &recs[0], sl, layout); + let (mut v, _) = analyze(&bm, &rs, &lk, false); + assert_eq!(v.len(), 1); + assert_eq!(v[0].index, 3); + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout, format: OutputFormat::Fasta, progress: false }; + pass2_extract(&m, &mut v, &ep).unwrap(); + let c = std::fs::read_to_string(o).unwrap(); + let l: Vec<&str> = c.lines().collect(); + assert_eq!(l[1], "C"); assert_eq!(l[3], "T"); + std::fs::remove_file(&p).ok(); std::fs::remove_file(o).ok(); + } + + #[test] fn test_inconsistent_length_error() { + // Sequences of differing lengths are a malformed alignment β†’ error. + let p = tmp("badlen", ">s1\nATGC\n>s2\nATG\n"); + let m = setup(&p); + assert!(index_fasta(&m).is_err()); + std::fs::remove_file(&p).ok(); + } + + #[test] fn test_crlf_single_line() { + // CRLF line endings on single-line sequences use the fast path. + let p = tmp("crlf1", ""); + std::fs::write(&p, b">s1\r\nATGC\r\n>s2\r\nATCC\r\n").unwrap(); + let m = setup(&p); + let lk = build_lookup(false); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + assert_eq!(sl, 4); + assert!(layout.single_line); + let bm = pass1_scan(&m, &recs, sl, layout, &lk); + let rs = get_ref_seq(&m, &recs[0], sl, layout); + let (v, _) = analyze(&bm, &rs, &lk, false); + assert_eq!(v.len(), 1); + assert_eq!(v[0].index, 2); + std::fs::remove_file(&p).ok(); + } + #[test] fn test_empty() { let p = tmp("empg", ""); let m = setup(&p); diff --git a/src/scan.rs b/src/scan.rs index 2e18df9..8c611a2 100644 --- a/src/scan.rs +++ b/src/scan.rs @@ -3,22 +3,36 @@ //! Builds a per-position bitmask by OR-ing each sequence's nucleotide flags, //! then classifies positions as variable (>1 allele), constant, or ambiguous. +#[cfg(feature = "parallel")] use rayon::prelude::*; use crate::fasta::FastaRecord; use crate::types::*; -/// Prefault mmap pages by touching one byte per OS page. -/// Eliminates soft page faults during the scan loop (~0.5s on 1 GB files). +/// Minimum total scan work (records Γ— positions) before the parallel path is +/// worth its overhead β€” roughly 50 seqs Γ— 4M bp. +#[cfg(feature = "parallel")] +const PARALLEL_MIN_WORK: usize = 200_000_000; + +/// Prefault mmap pages by touching one byte per OS page. Eliminates soft page faults during +/// the scan loop. For large files this is parallelised so it is no longer a single-core stall +/// (~10 s on a 20 GB alignment) ahead of the parallel scan. #[inline(never)] fn prefault(data: &[u8]) { const PAGE: usize = 4096; - let mut sum = 0u8; - let mut off = 0; - while off < data.len() { - sum = sum.wrapping_add(data[off]); - off += PAGE; - } + let n = data.len(); + let num_pages = n.div_ceil(PAGE); + #[cfg(feature = "parallel")] + let sum: u8 = if n >= 256 * 1024 * 1024 { + (0..num_pages) + .into_par_iter() + .map(|p| data[p * PAGE]) + .reduce(|| 0u8, |a, b| a.wrapping_add(b)) + } else { + (0..num_pages).fold(0u8, |s, p| s.wrapping_add(data[p * PAGE])) + }; + #[cfg(not(feature = "parallel"))] + let sum: u8 = (0..num_pages).fold(0u8, |s, p| s.wrapping_add(data[p * PAGE])); std::hint::black_box(sum); } @@ -31,54 +45,85 @@ pub fn pass1_scan( data: &[u8], records: &[FastaRecord], seq_length: usize, layout: SeqLayout, lookup: &[u8; 256], ) -> Vec { - let mut bitmask = vec![0u8; seq_length]; - // Prefault all pages into RAM before the hot loop prefault(data); - // Parallel: each thread scans a chunk of sequences into its own bitmask, - // then merge all partial bitmasks with OR. Threads share the mmap read-only. - let num_threads = rayon::current_num_threads().min(records.len()); - - // Parallelism only pays off when there's enough work per thread. - // Threshold: total scan work > ~200M bases (e.g., 50 seqs Γ— 4M bp). - let total_work = records.len() * seq_length; - if num_threads <= 1 || total_work < 200_000_000 { - // Sequential fallback for small inputs - scan_sequential(data, records, seq_length, layout, lookup, &mut bitmask); - } else { - // Split records into chunks, one per thread - let chunk_size = records.len().div_ceil(num_threads); - let partial_bitmasks: Vec> = records - .par_chunks(chunk_size) - .map(|chunk| { - let mut local_bm = vec![0u8; seq_length]; - scan_sequential(data, chunk, seq_length, layout, lookup, &mut local_bm); - local_bm - }) - .collect(); - - // Merge: OR all partial bitmasks into the final one - for partial in &partial_bitmasks { - for (bm, &p) in bitmask.iter_mut().zip(partial.iter()) { - *bm |= p; - } + // Parallelism only pays off when there's enough work per thread (~200M bases). Small + // inputs (and builds without the `parallel` feature, e.g. wasm) scan sequentially. + #[cfg(feature = "parallel")] + { + let num_threads = rayon::current_num_threads().min(records.len()); + let total_work = records.len() * seq_length; + if num_threads > 1 && total_work >= PARALLEL_MIN_WORK { + return scan_parallel(data, records, seq_length, layout, lookup, num_threads); } } + let mut bitmask = vec![0u8; seq_length]; + scan_sequential(data, records, seq_length, layout, lookup, &mut bitmask); + bitmask +} + +/// Parallel scan: each thread scans a disjoint chunk of records into its own +/// bitmask, then all partials are merged with OR. OR is commutative and +/// associative over the disjoint chunks, so the result is byte-for-byte +/// identical to a sequential scan. +#[cfg(feature = "parallel")] +fn scan_parallel( + data: &[u8], records: &[FastaRecord], seq_length: usize, + layout: SeqLayout, lookup: &[u8; 256], num_threads: usize, +) -> Vec { + let chunk_size = records.len().div_ceil(num_threads); + let partial_bitmasks: Vec> = records + .par_chunks(chunk_size) + .map(|chunk| { + let mut local_bm = vec![0u8; seq_length]; + scan_sequential(data, chunk, seq_length, layout, lookup, &mut local_bm); + local_bm + }) + .collect(); + let mut bitmask = vec![0u8; seq_length]; + for partial in &partial_bitmasks { + for (bm, &p) in bitmask.iter_mut().zip(partial.iter()) { + *bm |= p; + } + } bitmask } +/// Branchless nucleotide β†’ bitmask for the standard A/C/G/T(+gap) alphabet. Unlike a table +/// gather, this auto-vectorises (SSE2/AVX2/NEON). `gap_mask` is 1 when gaps are counted, else 0. +/// Byte-identical to the strict lookup table for the ACGT(+gap) case. +#[inline(always)] +fn acgt_bits(b: u8, gap_mask: u8) -> u8 { + let u = b & !0x20; // fold letter case (A/a -> A) + (u == b'A') as u8 + | (((u == b'C') as u8) << 1) + | (((u == b'G') as u8) << 2) + | (((u == b'T') as u8) << 3) + | ((((b == b'-') as u8) & gap_mask) << 4) +} + /// Sequential scan of a set of records into a bitmask. fn scan_sequential( data: &[u8], records: &[FastaRecord], seq_length: usize, layout: SeqLayout, lookup: &[u8; 256], bitmask: &mut [u8], ) { + // The branchless kernel is only valid for the standard alphabet; IUPAC-resolve tables + // (where e.g. 'R' maps to a nonzero bitmask) fall back to the general gather. + let fast = lookup[b'R' as usize] == 0; + let gap_mask = if lookup[b'-' as usize] != 0 { 1u8 } else { 0u8 }; if layout.single_line { for rec in records { let seq = &data[rec.seq_offset..rec.seq_offset + seq_length]; - for (bm_byte, &seq_byte) in bitmask.iter_mut().zip(seq.iter()) { - *bm_byte |= lookup[seq_byte as usize]; + if fast { + for (bm_byte, &seq_byte) in bitmask.iter_mut().zip(seq.iter()) { + *bm_byte |= acgt_bits(seq_byte, gap_mask); + } + } else { + for (bm_byte, &seq_byte) in bitmask.iter_mut().zip(seq.iter()) { + *bm_byte |= lookup[seq_byte as usize]; + } } } } else { @@ -118,6 +163,7 @@ pub fn analyze( else if bits & BIT_C != 0 { cs.c += 1; } else if bits & BIT_G != 0 { cs.g += 1; } else if bits & BIT_T != 0 { cs.t += 1; } + else { ambiguous += 1; } // all-gap column (BIT_GAP only) under --include-gaps } else { ambiguous += 1; } @@ -126,3 +172,54 @@ pub fn analyze( let num_variable = vars.len(); (vars, SiteCounts { constant: cs, variable: num_variable, ambiguous }) } + +// ============================================================================= +// Tests +// ============================================================================= + +#[cfg(test)] +mod tests { + use super::*; + use crate::fasta::index_fasta; + + #[test] + fn acgt_bits_matches_lookup() { + // The branchless kernel must be byte-identical to the strict lookup table. + for gap in [false, true] { + let lk = build_lookup(gap); + let gm = if gap { 1u8 } else { 0u8 }; + for b in 0u16..256 { + assert_eq!(acgt_bits(b as u8, gm), lk[b as usize], "byte {}", b); + } + } + } + + #[cfg(feature = "parallel")] + #[test] + fn parallel_matches_sequential() { + // The parallel chunk/merge path is the production path for real + // (whole-genome) inputs but is gated behind a size threshold that no + // small test reaches. Drive it directly and assert it equals a + // sequential scan byte-for-byte. + let mut fa = Vec::new(); + for i in 0..97usize { + fa.extend_from_slice(format!(">s{}\n", i).as_bytes()); + let mut seq = vec![b'A'; 40]; + seq[i % 40] = b'C'; + seq[(i * 7) % 40] = b'G'; + seq[(i * 13) % 40] = b'T'; + fa.extend_from_slice(&seq); + fa.push(b'\n'); + } + let lk = build_lookup(false); + let (recs, sl, layout) = index_fasta(&fa).unwrap(); + + let mut seq_bm = vec![0u8; sl]; + scan_sequential(&fa, &recs, sl, layout, &lk, &mut seq_bm); + + let threads = rayon::current_num_threads().min(recs.len()).max(2); + let par_bm = scan_parallel(&fa, &recs, sl, layout, &lk, threads); + + assert_eq!(seq_bm, par_bm); + } +} diff --git a/src/types.rs b/src/types.rs index a8dab5a..8049313 100644 --- a/src/types.rs +++ b/src/types.rs @@ -10,8 +10,10 @@ pub const BIT_G: u8 = 0b00100; pub const BIT_T: u8 = 0b01000; pub const BIT_GAP: u8 = 0b10000; -/// Maximum alignment length (prevents OOM on malicious input). -pub const MAX_SEQ_LENGTH: usize = 10_000_000_000; +/// Maximum alignment length. Kept below `u32::MAX` so alignment columns and reference +/// positions (VCF `POS`) always fit in a `u32` without truncation; also bounds memory on +/// malicious input (the bitmask alone is one byte per column). +pub const MAX_SEQ_LENGTH: usize = 4_000_000_000; /// Maximum VCF genotype matrix size in bytes. pub const MAX_VCF_GENO_BYTES: usize = 4_000_000_000; @@ -19,6 +21,18 @@ pub const MAX_VCF_GENO_BYTES: usize = 4_000_000_000; /// I/O buffer size for BufWriter (16 MB). pub const IO_BUF: usize = 16 * 1024 * 1024; +/// Output format for the reduced alignment. +#[derive(clap::ValueEnum, Clone, Copy, Debug, PartialEq, Eq, Default)] +pub enum OutputFormat { + /// FASTA (default). + #[default] + Fasta, + /// Relaxed PHYLIP (IQ-TREE / RAxML). + Phylip, + /// NEXUS DATA block (MrBayes / PAUP* / SplitsTree). + Nexus, +} + /// Whether all sequences are single-line (no embedded newlines). /// When true, `data[seq_offset + pos]` gives the base at `pos` directly. /// When false, we must scan skipping newlines for each record. @@ -72,6 +86,29 @@ pub fn build_lookup(include_gaps: bool) -> [u8; 256] { t } +/// Build a nucleotide β†’ bitmask lookup that also resolves IUPAC ambiguity codes to their +/// constituent bases (R = A|G, etc.). Used only for pass-1 classification under +/// `--iupac-mode resolve`; NS counting and REF selection keep the strict table. +pub fn build_iupac_lookup(include_gaps: bool) -> [u8; 256] { + let mut t = build_lookup(include_gaps); + let mut set = |c: u8, bits: u8| { + t[c as usize] = bits; + t[(c + 32) as usize] = bits; // lowercase + }; + set(b'R', BIT_A | BIT_G); + set(b'Y', BIT_C | BIT_T); + set(b'S', BIT_C | BIT_G); + set(b'W', BIT_A | BIT_T); + set(b'K', BIT_G | BIT_T); + set(b'M', BIT_A | BIT_C); + set(b'B', BIT_C | BIT_G | BIT_T); + set(b'D', BIT_A | BIT_G | BIT_T); + set(b'H', BIT_A | BIT_C | BIT_T); + set(b'V', BIT_A | BIT_C | BIT_G); + // N and other fully-ambiguous codes stay 0. + t +} + /// Build lowercase β†’ uppercase lookup table. pub fn build_upper() -> [u8; 256] { let mut t = [0u8; 256]; diff --git a/src/vcf.rs b/src/vcf.rs index e69b877..2962004 100644 --- a/src/vcf.rs +++ b/src/vcf.rs @@ -3,26 +3,36 @@ //! Generates VCF v4.2 from the genotype matrix built during pass 2. //! Uses a per-position lookup table for O(1) allele β†’ index mapping. -use std::fs::File; use std::io::{self, BufWriter, Write}; +use crate::extract::open_sink; use crate::fasta::FastaRecord; -use crate::types::VariablePosition; +use crate::types::{VariablePosition, IO_BUF}; /// Write VCF output from genotype matrix and variable positions. +#[allow(clippy::too_many_arguments)] pub fn write_vcf( vcf_geno: &[u8], num_samples: usize, var_positions: &[VariablePosition], - vcf_path: &str, records: &[FastaRecord], seq_length: usize, + vcf_path: &str, records: &[FastaRecord], seq_length: usize, chrom: &str, reference: &str, + pos_map: Option<&[u32]>, ) -> io::Result<()> { - let out = File::create(vcf_path).map_err(|e| io::Error::new(e.kind(), - format!("Cannot create VCF '{}': {}", vcf_path, e)))?; - let mut w = BufWriter::with_capacity(4 * 1024 * 1024, out); + // `-` streams to stdout, like -o / --sites-output / --stats-json; otherwise a file. + let out = open_sink(vcf_path)?; + let mut w = BufWriter::with_capacity(IO_BUF, out); // Header writeln!(w, "##fileformat=VCFv4.2")?; writeln!(w, "##source=snpick v{}", env!("CARGO_PKG_VERSION"))?; - writeln!(w, "##reference=first_sequence")?; - writeln!(w, "##contig=", seq_length)?; + writeln!(w, "##reference={}", reference)?; + // With reference coordinates the contig length is the ungapped reference length. POS values + // are clamped to >= 1 below, so clamp the declared length the same way: an all-gap reference + // has ungapped length 0, which would otherwise declare `length=0` while emitting POS=1 + // records β€” a self-contradictory VCF that htslib rejects. + let contig_len = match pos_map { + Some(rp) => (*rp.last().unwrap_or(&(seq_length as u32))).max(1) as usize, + None => seq_length, + }; + writeln!(w, "##contig=", chrom, contig_len)?; writeln!(w, "##INFO=")?; writeln!(w, "##FORMAT=")?; write!(w, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT")?; @@ -32,30 +42,53 @@ pub fn write_vcf( } writeln!(w)?; - // Data rows + // Data rows. Each row is assembled in a reused byte buffer, so the + // num_var Γ— num_samples genotypes are single byte pushes rather than one + // formatted write!() per field. Output bytes are identical to the naive form. let mut lut = [255u8; 256]; + let mut row: Vec = Vec::with_capacity(64 + num_samples * 2); + let mut alt: Vec = Vec::with_capacity(8); for (vi, vp) in var_positions.iter().enumerate() { - let alt: String = vp.alt_bases.iter() - .map(|&b| if b == b'-' { "*".to_string() } else { (b as char).to_string() }) - .collect::>().join(","); + // ALT alleles, gap β†’ '*'. + alt.clear(); + for (i, &b) in vp.alt_bases.iter().enumerate() { + if i > 0 { alt.push(b','); } + alt.push(if b == b'-' { b'*' } else { b }); + } - // Build allele β†’ index LUT for this position + // Build allele β†’ index LUT for this position. lut[vp.ref_base as usize] = 0; for (i, &ab) in vp.alt_bases.iter().enumerate() { lut[ab as usize] = (i + 1) as u8; } - write!(w, "1\t{}\t.\t{}\t{}\t.\tPASS\tNS={}\tGT", - vp.index + 1, vp.ref_base as char, alt, vp.ns)?; + // VCF v4.2 REF must be A/C/G/T/N β€” never '-' and never '*' (which is an ALT-only + // spanning-deletion allele that bcftools/htslib reject in REF). Render a gap reference + // as 'N'. (ALT gaps stay '*', the snp-sites convention.) + let ref_byte = if vp.ref_base == b'-' { b'N' } else { vp.ref_base }; + + row.clear(); + let pos = match pos_map { + Some(rp) => rp[vp.index].max(1), + None => (vp.index + 1) as u32, + }; + write!(row, "{}\t{}\t.\t", chrom, pos)?; + row.push(ref_byte); + row.push(b'\t'); + row.extend_from_slice(&alt); + write!(row, "\t.\tPASS\tNS={}\tGT", vp.ns)?; - let row = vi * num_samples; + let base = vi * num_samples; for si in 0..num_samples { - let idx = lut[vcf_geno[row + si] as usize]; - if idx == 255 { write!(w, "\t.")?; } else { write!(w, "\t{}", idx)?; } + row.push(b'\t'); + let idx = lut[vcf_geno[base + si] as usize]; + // idx is 0..=4 (ref + at most 4 alts), so a single ASCII digit. + if idx == 255 { row.push(b'.'); } else { row.push(b'0' + idx); } } - writeln!(w)?; + row.push(b'\n'); + w.write_all(&row)?; - // Reset LUT entries + // Reset LUT entries. lut[vp.ref_base as usize] = 255; for &ab in &vp.alt_bases { lut[ab as usize] = 255; } } diff --git a/tests/cli.rs b/tests/cli.rs new file mode 100644 index 0000000..89052a7 --- /dev/null +++ b/tests/cli.rs @@ -0,0 +1,106 @@ +//! End-to-end CLI tests that drive the built `snpick` binary. Cargo sets +//! `CARGO_BIN_EXE_snpick` for integration tests, so these exercise argument +//! parsing and output routing exactly as a user would. + +use std::fs; +use std::path::PathBuf; +use std::process::Command; + +const BIN: &str = env!("CARGO_BIN_EXE_snpick"); + +/// A unique temp directory per test, removed on drop. +struct TempDir(PathBuf); +impl TempDir { + fn new(tag: &str) -> Self { + let p = std::env::temp_dir().join(format!("snpick-cli-{}-{}", tag, std::process::id())); + let _ = fs::remove_dir_all(&p); + fs::create_dir_all(&p).unwrap(); + TempDir(p) + } + fn path(&self, name: &str) -> PathBuf { + self.0.join(name) + } + fn write(&self, name: &str, contents: &str) -> PathBuf { + let p = self.path(name); + fs::write(&p, contents).unwrap(); + p + } +} +impl Drop for TempDir { + fn drop(&mut self) { + let _ = fs::remove_dir_all(&self.0); + } +} + +const ALN: &str = ">ref\nATGCAT\n>s1\nATGTAT\n>s2\nACGCAT\n>s3\nATGCAG\n"; + +#[test] +fn vcf_output_dash_streams_to_stdout() { + let d = TempDir::new("vcfdash"); + let fa = d.write("aln.fa", ALN); + let out = d.path("out.fa"); + let o = Command::new(BIN) + .current_dir(&d.0) + .args(["-f", fa.to_str().unwrap(), "-o", out.to_str().unwrap(), "--vcf-output", "-", "-q"]) + .output() + .unwrap(); + assert!(o.status.success(), "stderr: {}", String::from_utf8_lossy(&o.stderr)); + // The VCF must arrive on stdout, not in a file literally named "-". + let stdout = String::from_utf8_lossy(&o.stdout); + assert!(stdout.starts_with("##fileformat=VCFv4.2"), "stdout was: {:?}", stdout); + assert!(!d.path("-").exists(), "a file named '-' should NOT have been created"); + // The reduced FASTA still went to its file. + assert!(out.exists()); +} + +#[test] +fn check_does_not_reject_output_path() { + // --check writes nothing, so pointing -o at the input must not trip the collision guard, + // matching --dry-run's behavior. + let d = TempDir::new("checkpath"); + let fa = d.write("aln.fa", ALN); + for mode in ["--check", "--dry-run"] { + let o = Command::new(BIN) + .args(["-f", fa.to_str().unwrap(), mode, "-o", fa.to_str().unwrap(), "-q"]) + .output() + .unwrap(); + assert!(o.status.success(), "{mode} -o should succeed; stderr: {}", + String::from_utf8_lossy(&o.stderr)); + } +} + +#[test] +fn multiple_stdout_outputs_rejected() { + // Two sidecars both aimed at stdout would concatenate JSON+TSV into one unparseable stream. + let d = TempDir::new("multistdout"); + let fa = d.write("aln.fa", ALN); + let out = d.path("out.fa"); + let o = Command::new(BIN) + .args(["-f", fa.to_str().unwrap(), "-o", out.to_str().unwrap(), + "--stats-json", "-", "--sites-output", "-", "-q"]) + .output() + .unwrap(); + assert_eq!(o.status.code(), Some(2), "stdout: {}", String::from_utf8_lossy(&o.stdout)); + assert!(String::from_utf8_lossy(&o.stderr).contains("stdout")); + // Exactly one sidecar to stdout is still fine. + let ok = Command::new(BIN) + .args(["-f", fa.to_str().unwrap(), "-o", out.to_str().unwrap(), "--stats-json", "-", "-q"]) + .output() + .unwrap(); + assert!(ok.status.success()); + assert!(String::from_utf8_lossy(&ok.stdout).starts_with('{')); +} + +#[test] +fn vcf_derive_rejects_stdout_alignment() { + // `-o - --vcf` cannot derive a sensible VCF filename; it must error, not write "-.vcf". + let d = TempDir::new("vcfderive"); + let fa = d.write("aln.fa", ALN); + let o = Command::new(BIN) + .current_dir(&d.0) + .args(["-f", fa.to_str().unwrap(), "-o", "-", "--vcf", "-q"]) + .output() + .unwrap(); + assert_eq!(o.status.code(), Some(2), "stderr: {}", String::from_utf8_lossy(&o.stderr)); + assert!(!d.path("-.vcf").exists(), "no '-.vcf' file should be created"); +} diff --git a/tests/proptest.rs b/tests/proptest.rs new file mode 100644 index 0000000..8637074 --- /dev/null +++ b/tests/proptest.rs @@ -0,0 +1,63 @@ +//! Property-based tests over randomized alignments. + +use proptest::prelude::*; +use snpick::fasta::{get_ref_seq, index_fasta}; +use snpick::scan::{analyze, pass1_scan}; +use snpick::types::build_lookup; + +/// Build an equal-length FASTA of `nseq` sequences of `len` bases, drawn deterministically +/// from `seed` over the alphabet A/C/G/T/N/-. +fn make_alignment(nseq: usize, len: usize, seed: u64) -> Vec { + const ALPHABET: &[u8] = b"ACGTN-"; + let mut fa = Vec::new(); + let mut x = seed | 1; + for i in 0..nseq { + fa.extend_from_slice(format!(">s{}\n", i).as_bytes()); + for _ in 0..len { + x = x.wrapping_mul(6364136223846793005).wrapping_add(1); + fa.push(ALPHABET[(x >> 33) as usize % ALPHABET.len()]); + } + fa.push(b'\n'); + } + fa +} + +proptest! { + /// The core invariant: every alignment column is classified exactly once. + #[test] + fn classification_covers_every_column( + nseq in 2usize..8, + len in 1usize..40, + seed in any::(), + ) { + let fa = make_alignment(nseq, len, seed); + for gaps in [false, true] { + let lookup = build_lookup(gaps); + let (recs, sl, layout) = index_fasta(&fa).unwrap(); + let bm = pass1_scan(&fa, &recs, sl, layout, &lookup); + let rs = get_ref_seq(&fa, &recs[0], sl, layout); + let (v, sc) = analyze(&bm, &rs, &lookup, gaps); + prop_assert_eq!(v.len() + sc.constant.total() + sc.ambiguous, sl); + } + } + + /// The parallel and sequential scans always agree (production path vs. fallback). + #[test] + fn variable_count_stable_across_gaps_flag( + nseq in 2usize..8, + len in 1usize..40, + seed in any::(), + ) { + let fa = make_alignment(nseq, len, seed); + let lookup = build_lookup(false); + let (recs, sl, layout) = index_fasta(&fa).unwrap(); + let bm = pass1_scan(&fa, &recs, sl, layout, &lookup); + let rs = get_ref_seq(&fa, &recs[0], sl, layout); + let (v, _) = analyze(&bm, &rs, &lookup, false); + // Every reported variable position must be within bounds and have >= 1 ALT allele. + for vp in &v { + prop_assert!(vp.index < sl); + prop_assert!(!vp.alt_bases.is_empty()); + } + } +}