diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml new file mode 100644 index 0000000..c7bf005 --- /dev/null +++ b/.github/workflows/release.yml @@ -0,0 +1,127 @@ +name: Build & Release + +on: + pull_request: + branches: [ main ] + types: [ closed ] + +permissions: + contents: write + +env: + CARGO_TERM_COLOR: always + +jobs: + # Only run if the PR was actually merged (not just closed) + release: + if: github.event.pull_request.merged == true + runs-on: ubuntu-latest + outputs: + version: ${{ steps.version.outputs.version }} + tag: ${{ steps.version.outputs.tag }} + steps: + - uses: actions/checkout@v4 + + - name: Get version from Cargo.toml + id: version + run: | + VERSION=$(grep '^version' Cargo.toml | head -1 | sed 's/.*"\(.*\)"/\1/') + echo "version=$VERSION" >> "$GITHUB_OUTPUT" + echo "tag=v$VERSION" >> "$GITHUB_OUTPUT" + echo "πŸ“¦ Version: $VERSION" + + - name: Check if tag already exists + 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" + fi + + # Build matrix for cross-platform binaries + build: + needs: release + if: needs.release.outputs.version != '' + strategy: + fail-fast: false + matrix: + include: + - target: x86_64-unknown-linux-gnu + os: ubuntu-latest + name: snpick-linux-x86_64 + - target: aarch64-unknown-linux-gnu + os: ubuntu-latest + name: snpick-linux-aarch64 + - target: x86_64-apple-darwin + os: macos-13 + name: snpick-macos-x86_64 + - target: aarch64-apple-darwin + os: macos-14 + name: snpick-macos-aarch64 + + runs-on: ${{ matrix.os }} + steps: + - uses: actions/checkout@v4 + + - name: Install Rust toolchain + uses: dtolnay/rust-toolchain@stable + with: + targets: ${{ matrix.target }} + + - name: Install cross-compilation tools (Linux aarch64) + if: matrix.target == 'aarch64-unknown-linux-gnu' + run: | + sudo apt-get update + sudo apt-get install -y gcc-aarch64-linux-gnu + + - name: Build release binary + run: cargo build --release --target ${{ matrix.target }} + env: + CARGO_TARGET_AARCH64_UNKNOWN_LINUX_GNU_LINKER: aarch64-linux-gnu-gcc + + - name: Package binary + run: | + mkdir -p dist + cp target/${{ matrix.target }}/release/snpick dist/${{ matrix.name }} + chmod +x dist/${{ matrix.name }} + + - name: Upload artifact + uses: actions/upload-artifact@v4 + with: + name: ${{ matrix.name }} + path: dist/${{ matrix.name }} + + # Create GitHub release with all binaries + publish: + needs: [release, build] + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Download all artifacts + uses: actions/download-artifact@v4 + with: + path: binaries + + - name: List binaries + run: find binaries -type f | sort + + - name: Generate checksums + run: | + cd binaries + find . -type f -name 'snpick-*' -exec sh -c ' + for f; do sha256sum "$f"; done + ' _ {} + | sort > ../SHA256SUMS.txt + cat ../SHA256SUMS.txt + + - name: Create release + uses: softprops/action-gh-release@v2 + with: + tag_name: ${{ needs.release.outputs.tag }} + name: SNPick ${{ needs.release.outputs.tag }} + generate_release_notes: true + files: | + binaries/snpick-linux-x86_64/snpick-linux-x86_64 + binaries/snpick-linux-aarch64/snpick-linux-aarch64 + binaries/snpick-macos-x86_64/snpick-macos-x86_64 + binaries/snpick-macos-aarch64/snpick-macos-aarch64 + SHA256SUMS.txt diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index 9fd45e0..9fc0c18 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -2,7 +2,7 @@ name: Rust on: push: - branches: [ "main" ] + branches: [ "main", "1.0.1" ] pull_request: branches: [ "main" ] @@ -11,12 +11,18 @@ env: jobs: build: - runs-on: ubuntu-latest - steps: - uses: actions/checkout@v4 + + - name: Clippy + run: cargo clippy -- -D warnings + - name: Build run: cargo build --verbose - - name: Run tests + + - name: Test run: cargo test --verbose + + - name: Build (release) + run: cargo build --release diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e7035df --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +target/ +test_*.fasta +test_*.vcf +*.vcf +!src/ diff --git a/Cargo.lock b/Cargo.lock index 412cf73..b413335 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2,30 +2,6 @@ # It is not intended for manual editing. version = 4 -[[package]] -name = "aho-corasick" -version = "1.1.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8e60d3430d3a69478ad0993f19238d2df97c507009a52b3c10addcd7f6bcb916" -dependencies = [ - "memchr", -] - -[[package]] -name = "android-tzdata" -version = "0.1.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e999941b234f3131b00bc13c22d06e8c5ff726d1b6318ac7eb276997bbb4fef0" - -[[package]] -name = "android_system_properties" -version = "0.1.5" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "819e7219dbd41043ac279b19830f2efc897156490d7fd6ea916720117ee66311" -dependencies = [ - "libc", -] - [[package]] name = "anstream" version = "0.6.15" @@ -75,156 +51,11 @@ dependencies = [ "windows-sys", ] -[[package]] -name = "anyhow" -version = "1.0.89" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "86fdf8605db99b54d3cd748a44c6d04df638eb5dafb219b135d0149bd0db01f6" - -[[package]] -name = "approx" -version = "0.5.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" -dependencies = [ - "num-traits", -] - -[[package]] -name = "autocfg" -version = "1.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ace50bade8e6234aa140d9a2f552bbee1db4d353f69b8217bc503490fc1a9f26" - -[[package]] -name = "bio" -version = "0.41.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4a69423e30444738eccc5e54eccee75779dd3f15ecc0469b95d8529d4b6b7586" -dependencies = [ - "anyhow", - "approx", - "bio-types", - "bit-set", - "bv", - "bytecount", - "csv", - "custom_derive", - "enum-map", - "fxhash", - "getset", - "itertools", - "itertools-num", - "lazy_static", - "multimap", - "ndarray", - "newtype_derive", - "num-integer", - "num-traits", - "ordered-float", - "petgraph", - "rand", - "regex", - "serde", - "serde_derive", - "statrs", - "strum", - "strum_macros 0.23.1", - "thiserror", - "triple_accel", - "vec_map", -] - -[[package]] -name = "bio-types" -version = "1.0.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f4dcf54f8b7f51450207d54780bab09c05f30b8b0caa991545082842e466ad7e" -dependencies = [ - "derive-new", - "lazy_static", - "regex", - "strum_macros 0.26.4", - "thiserror", -] - -[[package]] -name = "bit-set" -version = "0.5.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0700ddab506f33b20a03b13996eccd309a48e5ff77d0d95926aa0210fb4e95f1" -dependencies = [ - "bit-vec", -] - -[[package]] -name = "bit-vec" -version = "0.6.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "349f9b6a179ed607305526ca489b34ad0a41aed5f7980fa90eb03160b69598fb" - -[[package]] -name = "bumpalo" -version = "3.16.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "79296716171880943b8470b5f8d03aa55eb2e645a4874bdbb28adb49162e012c" - -[[package]] -name = "bv" -version = "0.11.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8834bb1d8ee5dc048ee3124f2c7c1afcc6bc9aed03f11e9dfd8c69470a5db340" -dependencies = [ - "feature-probe", - "serde", -] - -[[package]] -name = "bytecount" -version = "0.6.8" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5ce89b21cab1437276d2650d57e971f9d548a2d9037cc231abdc0562b97498ce" - -[[package]] -name = "byteorder" -version = "1.5.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" - -[[package]] -name = "cc" -version = "1.1.30" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b16803a61b81d9eabb7eae2588776c4c1e584b738ede45fdbb4c972cec1e9945" -dependencies = [ - "shlex", -] - -[[package]] -name = "cfg-if" -version = "1.0.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" - -[[package]] -name = "chrono" -version = "0.4.38" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a21f936df1771bf62b77f047b726c4625ff2e8aa607c01ec06e5a05bd8463401" -dependencies = [ - "android-tzdata", - "iana-time-zone", - "js-sys", - "num-traits", - "wasm-bindgen", - "windows-targets", -] - [[package]] name = "clap" -version = "4.5.20" +version = "4.5.60" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b97f376d85a664d5837dbae44bf546e6477a679ff6610010f17276f686d867e8" +checksum = "2797f34da339ce31042b27d23607e051786132987f595b02ba4f6a6dffb7030a" dependencies = [ "clap_builder", "clap_derive", @@ -232,9 +63,9 @@ dependencies = [ [[package]] name = "clap_builder" -version = "4.5.20" +version = "4.5.60" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "19bc80abd44e4bed93ca373a0704ccbd1b710dc5749406201bb018272808dc54" +checksum = "24a241312cea5059b13574bb9b3861cabf758b879c15190b37b6d6fd63ab6876" dependencies = [ "anstream", "anstyle", @@ -244,21 +75,21 @@ dependencies = [ [[package]] name = "clap_derive" -version = "4.5.18" +version = "4.5.55" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4ac6a0c7b1a9e9a5186361f67dfa1b88213572f427fb9ab038efb2bd8c582dab" +checksum = "a92793da1a46a5f2a02a6f4c46c6496b28c43638adea8306fcb0caa1634f24e5" dependencies = [ - "heck 0.5.0", + "heck", "proc-macro2", "quote", - "syn 2.0.79", + "syn", ] [[package]] name = "clap_lex" -version = "0.7.2" +version = "1.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1462739cb27611015575c0c11df5df7601141071f07518d56fcc1be504cbec97" +checksum = "c8d4a3bb8b1e0c1050499d1815f5ab16d04f0959b233085fb31653fbfc9d98f9" [[package]] name = "colorchoice" @@ -266,30 +97,11 @@ version = "1.0.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d3fd119d74b830634cea2a0f58bbd0d54540518a14397557951e79340abc28c0" -[[package]] -name = "console" -version = "0.15.8" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0e1f83fc076bd6dd27517eacdf25fef6c4dfe5f1d7448bafaaf3a26f13b5e4eb" -dependencies = [ - "encode_unicode", - "lazy_static", - "libc", - "unicode-width", - "windows-sys", -] - -[[package]] -name = "core-foundation-sys" -version = "0.8.7" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "773648b94d0e5d620f64f280777445740e61fe701025087ec8b57f45c791888b" - [[package]] name = "crossbeam-deque" -version = "0.8.5" +version = "0.8.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "613f8cc01fe9cf1a3eb3d7f488fd2fa8388403e97039e2f73692932e291a770d" +checksum = "9dd111b7b7f7d55b72c0a6ae361660ee5853c9af73f70c3c2ef6858b950e2e51" dependencies = [ "crossbeam-epoch", "crossbeam-utils", @@ -306,144 +118,15 @@ dependencies = [ [[package]] name = "crossbeam-utils" -version = "0.8.20" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "22ec99545bb0ed0ea7bb9b8e1e9122ea386ff8a48c0922e43f36d45ab09e0e80" - -[[package]] -name = "csv" -version = "1.3.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ac574ff4d437a7b5ad237ef331c17ccca63c46479e5b5453eb8e10bb99a759fe" -dependencies = [ - "csv-core", - "itoa", - "ryu", - "serde", -] - -[[package]] -name = "csv-core" -version = "0.1.11" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5efa2b3d7902f4b634a20cae3c9c4e6209dc4779feb6863329607560143efa70" -dependencies = [ - "memchr", -] - -[[package]] -name = "custom_derive" -version = "0.1.7" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ef8ae57c4978a2acd8b869ce6b9ca1dfe817bff704c220209fdef2c0b75a01b9" - -[[package]] -name = "derive-new" -version = "0.6.0" +version = "0.8.21" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d150dea618e920167e5973d70ae6ece4385b7164e0d799fe7c122dd0a5d912ad" -dependencies = [ - "proc-macro2", - "quote", - "syn 2.0.79", -] +checksum = "d0a5c400df2834b80a4c3327b3aad3a4c4cd4de0629063962b03235697506a28" [[package]] name = "either" -version = "1.13.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0" - -[[package]] -name = "encode_unicode" -version = "0.3.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a357d28ed41a50f9c765dbfe56cbc04a64e53e5fc58ba79fbc34c10ef3df831f" - -[[package]] -name = "enum-map" -version = "1.1.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e893a7ba6116821058dec84a6fb14fb2a97cd8ce5fd0f85d5a4e760ecd7329d9" -dependencies = [ - "enum-map-derive", -] - -[[package]] -name = "enum-map-derive" -version = "0.6.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "84278eae0af6e34ff6c1db44c11634a694aafac559ff3080e4db4e4ac35907aa" -dependencies = [ - "proc-macro2", - "quote", - "syn 1.0.109", -] - -[[package]] -name = "equivalent" -version = "1.0.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5443807d6dff69373d433ab9ef5378ad8df50ca6298caf15de6e52e24aaf54d5" - -[[package]] -name = "feature-probe" -version = "0.1.1" +version = "1.15.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "835a3dc7d1ec9e75e2b5fb4ba75396837112d2060b03f7d43bc1897c7f7211da" - -[[package]] -name = "fixedbitset" -version = "0.4.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0ce7134b9999ecaf8bcd65542e436736ef32ddca1b3e06094cb6ec5755203b80" - -[[package]] -name = "fxhash" -version = "0.2.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c31b6d751ae2c7f11320402d34e41349dd1016f8d5d45e48c4312bc8625af50c" -dependencies = [ - "byteorder", -] - -[[package]] -name = "getrandom" -version = "0.2.15" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c4567c8db10ae91089c99af84c68c38da3ec2f087c3f82960bcdbf3656b6f4d7" -dependencies = [ - "cfg-if", - "libc", - "wasi", -] - -[[package]] -name = "getset" -version = "0.1.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f636605b743120a8d32ed92fc27b6cde1a769f8f936c065151eb66f88ded513c" -dependencies = [ - "proc-macro-error2", - "proc-macro2", - "quote", - "syn 2.0.79", -] - -[[package]] -name = "hashbrown" -version = "0.15.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1e087f84d4f86bf4b218b927129862374b72199ae7d8657835f1e89000eea4fb" - -[[package]] -name = "heck" -version = "0.3.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6d621efb26863f0e9924c6ac577e8275e5e6b77455db64ffa6c65c904e9e132c" -dependencies = [ - "unicode-segmentation", -] +checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719" [[package]] name = "heck" @@ -451,319 +134,25 @@ version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" -[[package]] -name = "iana-time-zone" -version = "0.1.61" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "235e081f3925a06703c2d0117ea8b91f042756fd6e7a6e5d901e8ca1a996b220" -dependencies = [ - "android_system_properties", - "core-foundation-sys", - "iana-time-zone-haiku", - "js-sys", - "wasm-bindgen", - "windows-core", -] - -[[package]] -name = "iana-time-zone-haiku" -version = "0.1.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f31827a206f56af32e590ba56d5d2d085f558508192593743f16b2306495269f" -dependencies = [ - "cc", -] - -[[package]] -name = "indexmap" -version = "2.6.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "707907fe3c25f5424cce2cb7e1cbcafee6bdbe735ca90ef77c29e84591e5b9da" -dependencies = [ - "equivalent", - "hashbrown", -] - -[[package]] -name = "indicatif" -version = "0.17.8" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "763a5a8f45087d6bcea4222e7b72c291a054edf80e4ef6efd2a4979878c7bea3" -dependencies = [ - "console", - "instant", - "number_prefix", - "portable-atomic", - "unicode-width", -] - -[[package]] -name = "instant" -version = "0.1.13" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e0242819d153cba4b4b05a5a8f2a7e9bbf97b6055b2a002b395c96b5ff3c0222" -dependencies = [ - "cfg-if", -] - [[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 = "itertools-num" -version = "0.1.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a872a22f9e6f7521ca557660adb96dd830e54f0f490fa115bb55dd69d38b27e7" -dependencies = [ - "num-traits", -] - -[[package]] -name = "itoa" -version = "1.0.11" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "49f1f14873335454500d59611f1cf4a4b0f786f9ac11f4312a78e4cf2566695b" - -[[package]] -name = "js-sys" -version = "0.3.72" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6a88f1bda2bd75b0452a14784937d796722fdebfe50df998aeb3f0b7603019a9" -dependencies = [ - "wasm-bindgen", -] - -[[package]] -name = "lazy_static" -version = "1.5.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe" - [[package]] name = "libc" -version = "0.2.159" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "561d97a539a36e26a9a5fad1ea11a3039a67714694aaa379433e580854bc3dc5" - -[[package]] -name = "libm" -version = "0.2.8" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4ec2a862134d2a7d32d7983ddcdd1c4923530833c9f2ea1a44fc5fa473989058" - -[[package]] -name = "log" -version = "0.4.22" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a7a70ba024b9dc04c27ea2f0c0548feb474ec5c54bba33a7f72f873a39d07b24" - -[[package]] -name = "matrixmultiply" -version = "0.3.9" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9380b911e3e96d10c1f415da0876389aaf1b56759054eeb0de7df940c456ba1a" -dependencies = [ - "autocfg", - "rawpointer", -] - -[[package]] -name = "memchr" -version = "2.7.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "78ca9ab1a0babb1e7d5695e3530886289c18cf2f87ec19a575a0abdce112e3a3" - -[[package]] -name = "multimap" -version = "0.8.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e5ce46fe64a9d73be07dcbe690a38ce1b293be448fd8ce1e6c1b8062c9f72c6a" -dependencies = [ - "serde", -] - -[[package]] -name = "nalgebra" -version = "0.27.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "462fffe4002f4f2e1f6a9dcf12cc1a6fc0e15989014efc02a941d3e0f5dc2120" -dependencies = [ - "approx", - "matrixmultiply", - "nalgebra-macros", - "num-complex", - "num-rational", - "num-traits", - "rand", - "rand_distr", - "simba", - "typenum", -] - -[[package]] -name = "nalgebra-macros" -version = "0.1.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "01fcc0b8149b4632adc89ac3b7b31a12fb6099a0317a4eb2ebff574ef7de7218" -dependencies = [ - "proc-macro2", - "quote", - "syn 1.0.109", -] - -[[package]] -name = "ndarray" -version = "0.15.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "adb12d4e967ec485a5f71c6311fe28158e9d6f4bc4a447b474184d0f91a8fa32" -dependencies = [ - "matrixmultiply", - "num-complex", - "num-integer", - "num-traits", - "rawpointer", -] - -[[package]] -name = "newtype_derive" -version = "0.1.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ac8cd24d9f185bb7223958d8c1ff7a961b74b1953fd05dba7cc568a63b3861ec" -dependencies = [ - "rustc_version", -] - -[[package]] -name = "ntapi" -version = "0.4.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e8a3895c6391c39d7fe7ebc444a87eb2991b2a0bc718fdabd071eec617fc68e4" -dependencies = [ - "winapi", -] - -[[package]] -name = "num-complex" -version = "0.4.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" -dependencies = [ - "num-traits", -] - -[[package]] -name = "num-integer" -version = "0.1.46" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f" -dependencies = [ - "num-traits", -] - -[[package]] -name = "num-rational" -version = "0.4.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" -dependencies = [ - "num-integer", - "num-traits", -] - -[[package]] -name = "num-traits" -version = "0.2.19" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" -dependencies = [ - "autocfg", - "libm", -] - -[[package]] -name = "number_prefix" -version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "830b246a0e5f20af87141b25c173cd1b609bd7779a4617d6ec582abaf90870f3" - -[[package]] -name = "once_cell" -version = "1.20.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1261fe7e33c73b354eab43b1273a57c8f967d0391e80353e51f764ac02cf6775" - -[[package]] -name = "ordered-float" -version = "1.1.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3305af35278dd29f46fcdd139e0b1fbfae2153f0e5928b39b035542dd31e37b7" -dependencies = [ - "num-traits", -] - -[[package]] -name = "paste" -version = "1.0.15" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" - -[[package]] -name = "petgraph" -version = "0.6.5" +version = "0.2.183" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b4c5cc86750666a3ed20bdaf5ca2a0344f9c67674cae0515bec2da16fbaa47db" -dependencies = [ - "fixedbitset", - "indexmap", -] - -[[package]] -name = "portable-atomic" -version = "1.9.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cc9c68a3f6da06753e9335d63e27f6b9754dd1920d941135b7ea8224f141adb2" - -[[package]] -name = "ppv-lite86" -version = "0.2.20" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "77957b295656769bb8ad2b6a6b09d897d94f05c41b069aede1fcdaa675eaea04" -dependencies = [ - "zerocopy", -] - -[[package]] -name = "proc-macro-error-attr2" -version = "2.0.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "96de42df36bb9bba5542fe9f1a054b8cc87e172759a1868aa05c1f3acc89dfc5" -dependencies = [ - "proc-macro2", - "quote", -] +checksum = "b5b646652bf6661599e1da8901b3b9522896f01e736bad5f723fe7a3a27f899d" [[package]] -name = "proc-macro-error2" -version = "2.0.1" +name = "memmap2" +version = "0.9.10" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "11ec05c52be0a07b08061f7dd003e7d7092e0472bc731b4af7bb1ef876109802" +checksum = "714098028fe011992e1c3962653c96b2d578c4b4bce9036e15ff220319b1e0e3" dependencies = [ - "proc-macro-error-attr2", - "proc-macro2", - "quote", - "syn 2.0.79", + "libc", ] [[package]] @@ -784,57 +173,11 @@ dependencies = [ "proc-macro2", ] -[[package]] -name = "rand" -version = "0.8.5" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "34af8d1a0e25924bc5b7c43c079c942339d8f0a8b57c39049bef581b46327404" -dependencies = [ - "libc", - "rand_chacha", - "rand_core", -] - -[[package]] -name = "rand_chacha" -version = "0.3.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88" -dependencies = [ - "ppv-lite86", - "rand_core", -] - -[[package]] -name = "rand_core" -version = "0.6.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c" -dependencies = [ - "getrandom", -] - -[[package]] -name = "rand_distr" -version = "0.4.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "32cb0b9bc82b0a0876c2dd994a7e7a2683d3e7390ca40e6886785ef0c7e3ee31" -dependencies = [ - "num-traits", - "rand", -] - -[[package]] -name = "rawpointer" -version = "0.2.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" - [[package]] name = "rayon" -version = "1.10.0" +version = "1.11.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b418a60154510ca1a002a752ca9714984e21e4241e804d32555251faf8b78ffa" +checksum = "368f01d005bf8fd9b1206fb6fa653e6c4a81ceb1466406b81792d87c5677a58f" dependencies = [ "either", "rayon-core", @@ -842,131 +185,21 @@ dependencies = [ [[package]] name = "rayon-core" -version = "1.12.1" +version = "1.13.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1465873a3dfdaa8ae7cb14b4383657caab0b3e8a0aa9ae8e04b044854c8dfce2" +checksum = "22e18b0f0062d30d4230b2e85ff77fdfe4326feb054b9783a3460d8435c8ab91" dependencies = [ "crossbeam-deque", "crossbeam-utils", ] -[[package]] -name = "regex" -version = "1.11.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "38200e5ee88914975b69f657f0801b6f6dccafd44fd9326302a4aaeecfacb1d8" -dependencies = [ - "aho-corasick", - "memchr", - "regex-automata", - "regex-syntax", -] - -[[package]] -name = "regex-automata" -version = "0.4.8" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "368758f23274712b504848e9d5a6f010445cc8b87a7cdb4d7cbee666c1288da3" -dependencies = [ - "aho-corasick", - "memchr", - "regex-syntax", -] - -[[package]] -name = "regex-syntax" -version = "0.8.5" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2b15c43186be67a4fd63bee50d0303afffcef381492ebe2c5d87f324e1b8815c" - -[[package]] -name = "rustc_version" -version = "0.1.7" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c5f5376ea5e30ce23c03eb77cbe4962b988deead10910c372b226388b594c084" -dependencies = [ - "semver", -] - -[[package]] -name = "rustversion" -version = "1.0.17" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "955d28af4278de8121b7ebeb796b6a45735dc01436d898801014aced2773a3d6" - -[[package]] -name = "ryu" -version = "1.0.18" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f3cb5ba0dc43242ce17de99c180e96db90b235b8a9fdc9543c96d2209116bd9f" - -[[package]] -name = "semver" -version = "0.1.20" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d4f410fedcf71af0345d7607d246e7ad15faaadd49d240ee3b24e5dc21a820ac" - -[[package]] -name = "serde" -version = "1.0.210" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c8e3592472072e6e22e0a54d5904d9febf8508f65fb8552499a1abc7d1078c3a" -dependencies = [ - "serde_derive", -] - -[[package]] -name = "serde_derive" -version = "1.0.210" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "243902eda00fad750862fc144cea25caca5e20d615af0a81bee94ca738f1df1f" -dependencies = [ - "proc-macro2", - "quote", - "syn 2.0.79", -] - -[[package]] -name = "shlex" -version = "1.3.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" - -[[package]] -name = "simba" -version = "0.5.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8e82063457853d00243beda9952e910b82593e4b07ae9f721b9278a99a0d3d5c" -dependencies = [ - "approx", - "num-complex", - "num-traits", - "paste", -] - [[package]] name = "snpick" -version = "0.1.0" +version = "1.0.1" dependencies = [ - "bio", - "chrono", "clap", - "indicatif", + "memmap2", "rayon", - "sysinfo", -] - -[[package]] -name = "statrs" -version = "0.15.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "05bdbb8e4e78216a85785a85d3ec3183144f98d0097b9281802c019bb07a6f05" -dependencies = [ - "approx", - "lazy_static", - "nalgebra", - "num-traits", - "rand", ] [[package]] @@ -975,232 +208,29 @@ version = "0.11.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7da8b5736845d9f2fcb837ea5d9e2628564b3b043a70948a3f0b778838c5fb4f" -[[package]] -name = "strum" -version = "0.23.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cae14b91c7d11c9a851d3fbc80a963198998c2a64eec840477fa92d8ce9b70bb" - -[[package]] -name = "strum_macros" -version = "0.23.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5bb0dc7ee9c15cea6199cde9a127fa16a4c5819af85395457ad72d68edc85a38" -dependencies = [ - "heck 0.3.3", - "proc-macro2", - "quote", - "rustversion", - "syn 1.0.109", -] - -[[package]] -name = "strum_macros" -version = "0.26.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4c6bee85a5a24955dc440386795aa378cd9cf82acd5f764469152d2270e581be" -dependencies = [ - "heck 0.5.0", - "proc-macro2", - "quote", - "rustversion", - "syn 2.0.79", -] - -[[package]] -name = "syn" -version = "1.0.109" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237" -dependencies = [ - "proc-macro2", - "quote", - "unicode-ident", -] - [[package]] name = "syn" -version = "2.0.79" +version = "2.0.87" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "89132cd0bf050864e1d38dc3bbc07a0eb8e7530af26344d3d2bbbef83499f590" +checksum = "25aa4ce346d03a6dcd68dd8b4010bcb74e54e62c90c573f394c46eae99aba32d" dependencies = [ "proc-macro2", "quote", "unicode-ident", ] -[[package]] -name = "sysinfo" -version = "0.29.11" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cd727fc423c2060f6c92d9534cef765c65a6ed3f428a03d7def74a8c4348e666" -dependencies = [ - "cfg-if", - "core-foundation-sys", - "libc", - "ntapi", - "once_cell", - "rayon", - "winapi", -] - -[[package]] -name = "thiserror" -version = "1.0.64" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d50af8abc119fb8bb6dbabcfa89656f46f84aa0ac7688088608076ad2b459a84" -dependencies = [ - "thiserror-impl", -] - -[[package]] -name = "thiserror-impl" -version = "1.0.64" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "08904e7672f5eb876eaaf87e0ce17857500934f4981c4a0ab2b4aa98baac7fc3" -dependencies = [ - "proc-macro2", - "quote", - "syn 2.0.79", -] - -[[package]] -name = "triple_accel" -version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "22048bc95dfb2ffd05b1ff9a756290a009224b60b2f0e7525faeee7603851e63" - -[[package]] -name = "typenum" -version = "1.17.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "42ff0bf0c66b8238c6f3b578df37d0b7848e55df8577b3f74f92a69acceeb825" - [[package]] name = "unicode-ident" version = "1.0.13" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e91b56cd4cadaeb79bbf1a5645f6b4f8dc5bde8834ad5894a8db35fda9efa1fe" -[[package]] -name = "unicode-segmentation" -version = "1.12.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f6ccf251212114b54433ec949fd6a7841275f9ada20dddd2f29e9ceea4501493" - -[[package]] -name = "unicode-width" -version = "0.1.14" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7dd6e30e90baa6f72411720665d41d89b9a3d039dc45b8faea1ddd07f617f6af" - [[package]] name = "utf8parse" version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "06abde3611657adf66d383f00b093d7faecc7fa57071cce2578660c9f1010821" -[[package]] -name = "vec_map" -version = "0.8.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f1bddf1187be692e79c5ffeab891132dfb0f236ed36a43c7ed39f1165ee20191" -dependencies = [ - "serde", -] - -[[package]] -name = "wasi" -version = "0.11.0+wasi-snapshot-preview1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423" - -[[package]] -name = "wasm-bindgen" -version = "0.2.95" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "128d1e363af62632b8eb57219c8fd7877144af57558fb2ef0368d0087bddeb2e" -dependencies = [ - "cfg-if", - "once_cell", - "wasm-bindgen-macro", -] - -[[package]] -name = "wasm-bindgen-backend" -version = "0.2.95" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cb6dd4d3ca0ddffd1dd1c9c04f94b868c37ff5fac97c30b97cff2d74fce3a358" -dependencies = [ - "bumpalo", - "log", - "once_cell", - "proc-macro2", - "quote", - "syn 2.0.79", - "wasm-bindgen-shared", -] - -[[package]] -name = "wasm-bindgen-macro" -version = "0.2.95" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e79384be7f8f5a9dd5d7167216f022090cf1f9ec128e6e6a482a2cb5c5422c56" -dependencies = [ - "quote", - "wasm-bindgen-macro-support", -] - -[[package]] -name = "wasm-bindgen-macro-support" -version = "0.2.95" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "26c6ab57572f7a24a4985830b120de1594465e5d500f24afe89e16b4e833ef68" -dependencies = [ - "proc-macro2", - "quote", - "syn 2.0.79", - "wasm-bindgen-backend", - "wasm-bindgen-shared", -] - -[[package]] -name = "wasm-bindgen-shared" -version = "0.2.95" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "65fc09f10666a9f147042251e0dda9c18f166ff7de300607007e96bdebc1068d" - -[[package]] -name = "winapi" -version = "0.3.9" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5c839a674fcd7a98952e593242ea400abe93992746761e38641405d28b00f419" -dependencies = [ - "winapi-i686-pc-windows-gnu", - "winapi-x86_64-pc-windows-gnu", -] - -[[package]] -name = "winapi-i686-pc-windows-gnu" -version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ac3b87c63620426dd9b991e5ce0329eff545bccbbb34f3be09ff6fb6ab51b7b6" - -[[package]] -name = "winapi-x86_64-pc-windows-gnu" -version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" - -[[package]] -name = "windows-core" -version = "0.52.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "33ab640c8d7e35bf8ba19b884ba838ceb4fba93a4e8c65a9059d08afcfc683d9" -dependencies = [ - "windows-targets", -] - [[package]] name = "windows-sys" version = "0.52.0" @@ -1273,24 +303,3 @@ name = "windows_x86_64_msvc" version = "0.52.6" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "589f6da84c646204747d1270a2a5661ea66ed1cced2631d546fdfb155959f9ec" - -[[package]] -name = "zerocopy" -version = "0.7.35" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1b9b4fd18abc82b8136838da5d50bae7bdea537c574d8dc1a34ed098d6c166f0" -dependencies = [ - "byteorder", - "zerocopy-derive", -] - -[[package]] -name = "zerocopy-derive" -version = "0.7.35" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "fa4f8080344d4671fb4e831a13ad1e68092748387dfc4f55e356242fae12ce3e" -dependencies = [ - "proc-macro2", - "quote", - "syn 2.0.79", -] diff --git a/Cargo.toml b/Cargo.toml index e0c11dc..257d500 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,16 +1,14 @@ [package] name = "snpick" -version = "1.0.0" +version = "1.0.1" edition = "2021" [dependencies] -bio = "2.0.3" -nalgebra = "0.33.2" -rayon = "1.10.0" clap = { version = "4.5.21", features = ["derive"] } -indicatif = "0.17.9" -chrono = "0.4.38" -sysinfo = "0.32.0" -twox-hash = "2.0.1" - +memmap2 = "0.9.10" +rayon = "1.11.0" +[profile.release] +opt-level = 3 +lto = "fat" +codegen-units = 1 diff --git a/README.md b/README.md index 76867d0..905337a 100644 --- a/README.md +++ b/README.md @@ -1,124 +1,147 @@

- + SNPick logo

- -[![License: GPL v3](https://img.shields.io/badge/license-GPL%20v3-%23af64d1?style=flat-square)](https://github.com/PathoGenOmics-Lab/snpick/blob/main/LICENSE) + +[![License: GPL v3](https://img.shields.io/badge/license-GPL%20v3-%23af64d1?style=flat-square)](LICENSE) +[![DOI](https://img.shields.io/badge/doi-10.5281%2Fzenodo.14191809-%23ff0077?style=flat-square)](https://doi.org/10.5281/zenodo.14191809) [![PGO](https://img.shields.io/badge/PathoGenOmics-lab-red?style=flat-square)](https://github.com/PathoGenOmics-Lab) -[![Anaconda-Server Badge](https://img.shields.io/conda/dn/bioconda/snpick.svg?style=flat-square)](https://anaconda.org/bioconda/snpick) [![Anaconda-Version Badge](https://anaconda.org/bioconda/snpick/badges/version.svg)](https://anaconda.org/bioconda/snpick) -[![DOI](https://img.shields.io/badge/doi-10.5281%2Fzenodo.14191809-%23ff0077?style=flat-square)](https://doi.org/10.5281/zenodo.14191809) +[![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) +**Fast, memory-efficient extraction of variable sites from FASTA alignments.** + +[Quick Start](#-quick-start) Β· [Features](#-features) Β· [Usage](#-usage) Β· [Benchmarks](#-benchmarks) Β· [Citation](#-citation) +
-__Paula Ruiz-Rodriguez1__ +__Paula Ruiz-Rodriguez1__ __and Mireia Coscolla1__
- 1. Institute for Integrative Systems Biology, I2SysBio, University of Valencia-CSIC, Valencia, Spain + 1. Institute for Integrative Systems Biology, I2SysBio, University of Valencia-CSIC, Valencia, Spain -# SNPick +--- -SNPick is a fast and memory-efficient tool designed to extract SNP (Single Nucleotide Polymorphism) sites from large-scale FASTA alignments, taking into consideration IUPAC ambiguous nucleotides. SNP analysis is critical for understanding genetic variation, phylogenetic relationships, and evolutionary biology. +## What is SNPick? -## Problematic +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. -Tools like **snp-sites** have been widely used to extract SNP positions from FASTA alignments. **snp-sites** is a valuable tool for smaller datasets, providing an effective solution for extracting SNPs in many research scenarios. However, it faces significant challenges when processing large alignments or a high number of sequences. **snp-sites** is not suitable for generating SNP datasets specifically designed for phylogenetic analyses of SNPs using ASC (Ascertainment Bias Correction) in tools like **RAxML** or **IQ-TREE**. +**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. -When dealing with thousands of sequences or large sequence lengths, **snp-sites** tends to become inefficient both in terms of runtime and memory usage. For instance, the tool may struggle or even fail entirely when processing datasets with hundreds of thousands of sequences or sequences that are several megabases in length. This limitation is particularly relevant for researchers who need to work with whole-genome alignments from many individuals, especially in epidemiological studies or population genomics. +### SNPick vs snp-sites -## πŸͺSNPick Advantage +| | **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 | +| ASC fconst output | βœ… Built-in | ❌ Not supported | +| VCF output | βœ… Optional | βœ… Default | +| Gap handling | βœ… Optional (`-g`) | βœ… Default | +| IUPAC ambiguous | βœ… Tracked as ambiguous | ⚠️ Treated as variant | -**SNPick** was developed to overcome these limitations, providing a highly scalable and optimized approach to SNP extraction. Unlike **snp-sites**, SNPick employs parallel processing and an efficient memory management strategy, making it suitable for very large datasets. This scalability allows users to work with massive alignments, extracting SNPs in a fraction of the time and with significantly lower memory requirements compared to traditional tools. +--- -The key benefits of SNPick include: -- **Scalability**: Able to handle datasets with tens of thousands of sequences and large genome alignments without a significant increase in memory consumption or processing time. -- **Efficiency**: Optimized for speed, taking advantage of multi-threading and efficient memory allocation. -- **Versatility**: Supports the extraction of SNPs while handling IUPAC ambiguous nucleotides, making it highly adaptable to different types of genomic data. +## πŸš€ Quick Start -## πŸ’Ύ Functionalities +```bash +# Install +conda install -c bioconda snpick -SNPick provides input that is compatible with phylogenetic software like **RAxML** and **IQ-TREE** when performing ascertainment bias correction (asc bias). This ensures that the output is suitable for accurate downstream phylogenetic inference and evolutionary analyses. +# Extract variable sites +snpick -f alignment.fasta -o snps.fasta -SNPick offers multiple functionalities to facilitate comprehensive SNP analysis and ensure compatibility with downstream applications: +# With VCF output +snpick -f alignment.fasta -o snps.fasta --vcf -### 1. **FASTA SNP Extraction** +# Include gaps as informative +snpick -f alignment.fasta -o snps.fasta -g +``` -SNPick allows users to extract SNP sites from a given FASTA alignment file, providing an output file that contains only the SNP positions. This feature makes downstream analysis more efficient by significantly reducing the size of the dataset while retaining only informative sites. +--- -### 2. **Multi-threading Support** +## ✨ Features -The tool supports multi-threading to expedite processing. Users can specify the number of threads using the `--threads` option, allowing SNPick to efficiently utilize available computational resources. +### Variable site extraction -### 3. **Handling IUPAC Ambiguous Nucleotides** +Identifies positions with more than one observed nucleotide across all sequences. Constant and ambiguous-only positions are excluded from the output. -SNPick considers IUPAC ambiguous nucleotides during the SNP extraction process. It can accurately identify variable positions even in the presence of ambiguous bases, ensuring that these nucleotides are appropriately handled. +### ASC bias correction support -### 4. **VCF Output Generation** +Reports constant site counts (`fconst`) directly, formatted for IQ-TREE's `+ASC` models: -SNPick can generate Variant Call Format (VCF) files from the extracted SNP positions. This feature is particularly useful for users who need to integrate SNP data into other bioinformatics pipelines, as VCF is a standard format widely used in genetic variation studies. +``` +[snpick] ASC fconst: 744123,1382922,1382180,743556 +``` -- To generate a VCF file, use the `--vcf` flag and specify an output VCF file with `--vcf-output`. The resulting VCF file will contain detailed information about each SNP position, including reference and alternate bases, and the base observed in each sample. +Use in IQ-TREE: +```bash +iqtree2 -s snps.fasta -m GTR+ASC -fconst 744123,1382922,1382180,743556 +``` -Example command to generate a VCF file: +### VCF generation -```sh -snpick --fasta input_alignment.fasta --output snp_alignment.fasta --vcf --vcf-output snp_output.vcf --threads 8 -``` +Optional VCF v4.2 output with per-sample genotypes. Reference allele taken from the first sequence. Ambiguous bases reported as missing (`.`). + +### IUPAC and gap handling -### 5. **Efficient Memory Usage** +- **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` -SNPick is designed to handle large-scale datasets without consuming excessive memory. By using efficient memory management and parallel processing, SNPick ensures that even extremely large alignments can be processed on standard computational infrastructure. -### 6. **Progress Monitoring** +### Parallel processing -The tool includes progress monitoring to keep users informed about the current status of the analysis. This is particularly useful for long-running processes on large datasets. SNPick provides a spinner and a progress bar to indicate the status of sequence reading, SNP identification, and writing processes. -Usage +Automatic multi-threaded scanning via Rayon when the dataset is large enough. Falls back to single-threaded for small inputs to avoid overhead. -SNPick can be easily integrated into existing pipelines for genomic data analysis. It is suitable for use in environments where memory efficiency and processing speed are crucial, such as analyzing large microbial datasets or conducting comparative genomic studies. +--- -## πŸ› οΈ Installation -You can install snpick via conda, mamba (for unix/mac) or downloading [the binary file](https://github.com/PathoGenOmics-Lab/snpick/releases/download/1.0.0/snpick) (unix): +## πŸ’Ύ Installation -### 🐍 Using conda -``` +### Bioconda (recommended) + +```bash conda install -c bioconda snpick -``` -### 🐍 Using mamba -``` +# or mamba install -c bioconda snpick ``` -### πŸ“¨ Using binary -``` -wget https://github.com/PathoGenOmics-Lab/snpick/releases/download/1.0.0/snpick + +### From source + +```bash +git clone https://github.com/PathoGenOmics-Lab/snpick.git +cd snpick +cargo build --release +# Binary at target/release/snpick ``` -## πŸ—ƒοΈ Arguments +### Pre-built binary (Linux) -SNPick provides a variety of arguments to customize the SNP extraction process: +```bash +wget https://github.com/PathoGenOmics-Lab/snpick/releases/latest/download/snpick +chmod +x snpick +``` -- `--fasta` (**required**): Specifies the input alignment file in FASTA format. This file contains the aligned sequences from which SNPs will be extracted. -- `--output` (**required**): Specifies the output file that will contain only the SNP sites extracted from the input alignment. -- `--threads` (**optional**): Specifies the number of threads to use for multi-threading. By default, SNPick uses 4 threads, but this can be adjusted to take advantage of additional computational resources for faster processing. -- `--vcf` (**optional**): If provided, generates a Variant Call Format (VCF) file with the SNP information. VCF is a standard format widely used for genetic variation studies. -- `--vcf-output` (**optional**): Specifies the output VCF file name. This argument is used in combination with `--vcf` to store detailed information about each SNP position, including reference and alternate bases, and the observed base for each sample. -- `--include-gaps` (**optional**): When this flag is used, SNPick will consider gaps (`-`) in the sequences as valid characters when determining variable positions. This can be useful for certain types of evolutionary analyses where gaps are informative. +--- -## 🎴 Example +## πŸ—ƒοΈ Usage -To extract SNPs from a given FASTA file: -```sh -snpick --fasta input_alignment.fasta --output snp_alignment.fasta --threads 8 ``` -**Input**: The `--fasta` parameter specifies the input alignment file in FASTA format. This file contains the aligned sequences from which SNPs will be extracted. +snpick [OPTIONS] --fasta --output +``` -**Output**: The `--output` parameter specifies the output file, which will contain only the SNP sites extracted from the input alignment. +| 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 | -This command will generate an output file containing only the SNP sites, optimized for downstream analysis. -### Input and Output Example +### Example -Input File (input_alignment.fasta): -```sh +**Input** (`alignment.fasta`): +``` >sequence1 ATGCTAGCTAGCTAGCTA >sequence2 @@ -126,35 +149,88 @@ ATGCTAGCTGGCTAGCTA >sequence3 ATGCTAGCTAGCTAGCTA ``` -Output File (snp_alignment.fasta): -```sh + +**Command:** +```bash +snpick -f alignment.fasta -o snps.fasta +``` + +**Output** (`snps.fasta`): +``` >sequence1 -T +A >sequence2 G >sequence3 -T +A +``` + +**stderr:** +``` +[snpick] Mapped 63 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. +``` + +--- + +## πŸ“Š Benchmarks + +Simulated *M. tuberculosis*-like genomes (4.4 Mbp, ~65% GC, 3.6% variable sites). + +### Scaling by number of sequences + +

+ Benchmark: sequence scaling +

+ +### Scaling by sequence length + +

+ Benchmark: length scaling +

+ +SNPick maintains **O(L)** memory regardless of sequence count, while snp-sites requires **O(NΓ—L)**. + +--- + +## πŸ—οΈ Architecture + +``` +Input FASTA ──mmap──▢ Index records ──▢ Pass 1: bitmask scan ──▢ Analyze + β”‚ (parallel) β”‚ + β”‚ β–Ό + └──────────▢ Pass 2: extract sites ──▢ FASTA + VCF + (sparse random access) ``` -In this example, the SNP site is at the position where sequence2 has a G while the others have a T. -### Another Example -If you have a large dataset and want to utilize all available CPU cores for faster processing, you can use the following command: -```sh -snpick --fasta large_dataset.fasta --output snp_large_output.fasta --threads 16 +- **Single memory-mapped file** shared across both passes β€” zero copies +- **Pass 1**: OR-based bitmask over all sequences (parallel with Rayon) +- **Pass 2**: only reads variable positions (sparse access via mmap) +- **Lookup tables**: 256-byte arrays for O(1) nucleotide classification and case conversion + +--- + +## πŸ“ 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} +} ``` -In this example, large_dataset.fasta is the input file, snp_large_output.fasta is the output file, and --threads 16 specifies the use of 16 CPU threads to speed up the SNP extraction process. --- -

-✨ [Contributors]((https://github.com/PathoGenOmics-Lab/AMAP/graphs/contributors)) -

+

✨ Contributors

- - -
-SNPick is developed with ❀️ by: + +
@@ -168,9 +244,9 @@ SNPick is developed with ❀️ by: πŸ”¬ πŸ€” πŸ”£ - 🎨 + 🎨 πŸ”§ - @@ -178,19 +254,15 @@ SNPick is developed with ❀️ by: Mireia Coscolla
- πŸ” + πŸ” πŸ€” πŸ§‘β€πŸ« πŸ”¬ πŸ““ -
This project follows the [all-contributors](https://github.com/all-contributors/all-contributors) specification ([emoji key](https://allcontributors.org/docs/en/emoji-key)). - - - - ---- +
diff --git a/benchmarks/benchmark.png b/benchmarks/benchmark.png new file mode 100644 index 0000000..0834b60 Binary files /dev/null and b/benchmarks/benchmark.png differ diff --git a/benchmarks/benchmark.svg b/benchmarks/benchmark.svg new file mode 100644 index 0000000..fdfbf56 --- /dev/null +++ b/benchmarks/benchmark.svg @@ -0,0 +1,3357 @@ + + + + + + + + 2026-03-31T10:13:33.880130 + image/svg+xml + + + Matplotlib v3.9.4, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/benchmarks/benchmark_length.png b/benchmarks/benchmark_length.png new file mode 100644 index 0000000..a48e7ce Binary files /dev/null and b/benchmarks/benchmark_length.png differ diff --git a/benchmarks/benchmark_length.svg b/benchmarks/benchmark_length.svg new file mode 100644 index 0000000..c4c0b98 --- /dev/null +++ b/benchmarks/benchmark_length.svg @@ -0,0 +1,3908 @@ + + + + + + + + 2026-03-31T10:24:12.519707 + image/svg+xml + + + Matplotlib v3.9.4, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/benchmarks/results.tsv b/benchmarks/results.tsv new file mode 100644 index 0000000..d115655 --- /dev/null +++ b/benchmarks/results.tsv @@ -0,0 +1,10 @@ +seqs snpick_time_s snpsites_time_s speedup snpick_mem_mb snpsites_mem_mb +10 0.09 0.21 2.3x 39 20 +25 0.16 0.49 3.1x 45 21 +50 0.28 1.03 3.7x 66 30 +100 0.54 2.44 4.5x 76 56 +150 0.97 5.06 5.2x 87 96 +200 1.36 7.48 5.5x 97 150 +250 1.72 9.38 5.5x 105 213 +500 3.80 32.49 8.6x 156 715 +1000 10.27 NA NA 217 NA diff --git a/src/extract.rs b/src/extract.rs new file mode 100644 index 0000000..7c907fd --- /dev/null +++ b/src/extract.rs @@ -0,0 +1,99 @@ +//! Pass 2: variable site extraction and output FASTA generation. +//! +//! Reads only the variable positions from each sequence (sparse access) +//! and writes a reduced FASTA. Optionally collects a genotype matrix for VCF. + +use std::fs::File; +use std::io::{self, BufWriter, Write}; + +use crate::fasta::FastaRecord; +use crate::types::*; + +/// Parameters for variable site extraction (pass 2). +pub struct ExtractParams<'a> { + pub records: &'a [FastaRecord<'a>], + pub output: &'a str, + pub collect_vcf: bool, + pub lookup: &'a [u8; 256], + pub upper: &'a [u8; 256], + pub layout: SeqLayout, +} + +/// Pass 2: extract variable sites from alignment and write output FASTA. +/// +/// For single-line FASTA: direct byte access via `data[seq_offset + pos]`. +/// For multi-line: linear scan per record, skipping newlines. +/// +/// Returns VCF genotype matrix if `collect_vcf` is true. +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 collect_vcf = *collect_vcf; + let layout = *layout; + 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 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() }; + let mut var_buf = vec![0u8; num_var]; + + for (si, rec) in records.iter().enumerate() { + if layout.single_line { + let base = rec.seq_offset; + for (vi, &p) in pos_indices.iter().enumerate() { + var_buf[vi] = upper[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] { + var_buf[var_idx] = upper[b as usize]; + var_idx += 1; + } + base_idx += 1; + } + } + + 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")?; + + if collect_vcf { + for (vi, &nuc) in var_buf.iter().enumerate() { + vcf_geno[vi * num_samples + si] = nuc; + if lookup[nuc as usize] != 0 { + ns_counts[vi] += 1; + } + } + } + } + + writer.flush()?; + + if collect_vcf { + for (vi, vp) in var_positions.iter_mut().enumerate() { + vp.ns = ns_counts[vi]; + } + } + + 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 new file mode 100644 index 0000000..08f21fd --- /dev/null +++ b/src/fasta.rs @@ -0,0 +1,115 @@ +//! Zero-copy FASTA parser over memory-mapped data. +//! +//! Records are indexed by scanning for `>` headers and tracking sequence offsets. +//! No data is copied β€” IDs and descriptions are `&[u8]` slices into the mmap. + +use std::io; + +use crate::types::{SeqLayout, MAX_SEQ_LENGTH}; + +/// A FASTA record as zero-copy slices into memory-mapped data. +pub struct FastaRecord<'a> { + pub id: &'a [u8], + pub desc: &'a [u8], + pub seq_offset: usize, +} + +/// Index FASTA records from memory-mapped data. Zero-copy: stores `&[u8]` slices. +/// +/// Returns `(records, seq_length, layout)`. +/// All sequences must have the same length (alignment requirement). +pub fn index_fasta(data: &[u8]) -> io::Result<(Vec>, usize, SeqLayout)> { + if data.is_empty() { + return Err(io::Error::new(io::ErrorKind::InvalidData, "Input FASTA is empty (0 bytes).")); + } + + let mut records = Vec::new(); + let mut pos = 0; + let len = data.len(); + let mut seq_length = 0usize; + let mut is_single_line = true; + + while pos < len { + while pos < len && (data[pos] == b'\n' || data[pos] == b'\r') { pos += 1; } + if pos >= len { break; } + + if data[pos] != b'>' { + return Err(io::Error::new(io::ErrorKind::InvalidData, + format!("Expected '>' at byte {}, got '{}'.", pos, data[pos] as char))); + } + pos += 1; + + // Header line + let header_start = pos; + while pos < len && data[pos] != b'\n' && data[pos] != b'\r' { pos += 1; } + let header = &data[header_start..pos]; + if pos < len && data[pos] == b'\r' { pos += 1; } + if pos < len && data[pos] == b'\n' { pos += 1; } + + let (id, desc) = if let Some(sp) = header.iter().position(|&b| b == b' ' || b == b'\t') { + (&header[..sp], &header[sp + 1..]) + } else { + (header, &data[0..0]) + }; + + let seq_offset = pos; + + // Count bases + let mut seq_len = 0usize; + let mut line_count = 0usize; + while pos < len && data[pos] != b'>' { + let line_start = pos; + while pos < len && data[pos] != b'\n' && data[pos] != b'\r' { pos += 1; } + let line_len = pos - line_start; + if line_len > 0 { + seq_len += line_len; + line_count += 1; + } + 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; } + + if records.is_empty() { + if seq_len == 0 { + return Err(io::Error::new(io::ErrorKind::InvalidData, "First sequence has length 0.")); + } + if seq_len > MAX_SEQ_LENGTH { + return Err(io::Error::new(io::ErrorKind::InvalidData, + format!("Sequence length {} exceeds maximum.", seq_len))); + } + seq_length = seq_len; + } else if seq_len != seq_length { + let id_str = std::str::from_utf8(id).unwrap_or("?"); + return Err(io::Error::new(io::ErrorKind::InvalidData, + format!("Sequence '{}' (#{}) has length {} but expected {}.", + id_str, records.len() + 1, seq_len, seq_length))); + } + + records.push(FastaRecord { id, desc, seq_offset }); + } + + if records.is_empty() { + return Err(io::Error::new(io::ErrorKind::InvalidData, "Input FASTA is empty.")); + } + + Ok((records, seq_length, SeqLayout { single_line: is_single_line })) +} + +/// Extract reference sequence from first record. +pub fn get_ref_seq(data: &[u8], rec: &FastaRecord, seq_length: usize, layout: SeqLayout) -> Vec { + if layout.single_line { + data[rec.seq_offset..rec.seq_offset + seq_length].to_vec() + } else { + let mut seq = Vec::with_capacity(seq_length); + let mut pos = rec.seq_offset; + let end = data.len(); + while seq.len() < seq_length && pos < end { + let b = data[pos]; + pos += 1; + if b != b'\n' && b != b'\r' { seq.push(b); } + } + seq + } +} diff --git a/src/main.rs b/src/main.rs index 6a1d51b..f1bb6cd 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,423 +1,493 @@ -// Necessary imports -use bio::io::fasta; +mod extract; +mod fasta; +mod scan; +mod types; +mod vcf; + use clap::Parser; -use indicatif::{ProgressBar, ProgressStyle}; -use rayon::prelude::*; +use memmap2::Mmap; use std::fs::File; -use std::io::{self, BufReader, BufWriter, Write}; -use std::sync::{Arc, Mutex}; -use std::sync::atomic::{AtomicUsize, Ordering}; -use std::time::Duration; -use chrono::Local; -use sysinfo::{System}; -use std::collections::{HashSet, HashMap}; - -/// snpick: A tool to extract variable sites from a FASTA alignment and generate a VCF with actual bases, including ambiguous bases and codons. +use std::io::{self, BufWriter, 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; + +// ============================================================================= +// CLI +// ============================================================================= + #[derive(Parser, Debug)] #[command( name = "snpick", - version = "1.0.0", + version = env!("CARGO_PKG_VERSION"), author = "Paula Ruiz-Rodriguez ", - about = "A fast and efficient tool for extracting variable sites and generating a VCF with actual bases, including ambiguous bases and codons." + about = "A fast, memory-efficient tool for extracting variable sites from FASTA alignments." )] struct Args { - /// Input FASTA alignment file - #[arg(short, long, help = "Input FASTA alignment file")] - fasta: String, - - /// Output FASTA file with variable sites - #[arg(short, long, help = "Output FASTA file with variable sites")] - output: String, - - /// Number of threads to use (optional) - #[arg(short, long, default_value_t = 4, help = "Number of threads to use (optional)")] - threads: usize, - - /// Consider the '-' (gap) symbol in variable site detection - #[arg(short = 'g', long, help = "Consider the '-' (gap) symbol in variable site detection")] - include_gaps: bool, + #[arg(short, long)] fasta: String, + #[arg(short, long)] output: String, + #[arg(short = 'g', long)] include_gaps: bool, + #[arg(long)] vcf: bool, + #[arg(long)] vcf_output: Option, +} - /// Generate VCF file with variable sites - #[arg(long, help = "Generate VCF file with variable sites")] - vcf: bool, +// ============================================================================= +// Path validation +// ============================================================================= - /// Output VCF file (optional) - #[arg(long, help = "Output VCF file (optional)")] - vcf_output: Option, +fn resolve_path(p: &str) -> io::Result { + let path = Path::new(p); + if path.exists() { + return std::fs::canonicalize(path); + } + let parent = path.parent().unwrap_or(Path::new(".")); + let parent_abs = std::fs::canonicalize(parent).map_err(|e| { + io::Error::new(e.kind(), format!("Cannot resolve parent of '{}': {}", p, e)) + })?; + Ok(parent_abs.join(path.file_name().unwrap_or_default())) } -/// Converts a nucleotide to a bitmask -fn nucleotide_to_bit(nuc: u8, include_gaps: bool) -> Option { - match nuc.to_ascii_uppercase() { - b'A' => Some(0b000001), - b'C' => Some(0b000010), - b'G' => Some(0b000100), - b'T' => Some(0b001000), - b'-' if include_gaps => Some(0b010000), - _ => None, // Exclude ambiguous bases +fn check_paths_differ(a: &str, b: &str) -> io::Result<()> { + let pa = resolve_path(a)?; + let pb = resolve_path(b)?; + if pa == pb { + return Err(io::Error::new(io::ErrorKind::InvalidInput, + format!("Paths resolve to same file: {}", pa.display()))); } + Ok(()) } -fn main() -> io::Result<()> { - // Parse command-line arguments - let args = Args::parse(); +// ============================================================================= +// Pipeline +// ============================================================================= - let input_filename = args.fasta; - let output_filename = args.output; - let num_threads = args.threads; - let include_gaps = args.include_gaps; - let generate_vcf = args.vcf; - let vcf_output_filename = args.vcf_output.unwrap_or_else(|| "output.vcf".to_string()); - - // Configure a local thread pool for Rayon - let pool = rayon::ThreadPoolBuilder::new() - .num_threads(num_threads) - .build() - .expect("Failed to build Rayon thread pool"); - - // Initialize system for RAM usage reporting - let system = Arc::new(Mutex::new(System::new_all())); - - // Execute processing within the thread pool - pool.install(|| { - // Step 1: Identify variable positions and extract individual genotypes - println!("Starting Step 1: Identifying variable positions..."); - let (variable_positions_info, total_sequences, sample_names) = - identify_variable_positions(&input_filename, &system, include_gaps) - .expect("Failed to identify variable positions"); - println!( - "Step 1 Completed: Found {} variable positions.", - variable_positions_info.len() - ); - - if variable_positions_info.is_empty() { - eprintln!("No variable positions found in the alignment."); - std::process::exit(0); +fn run() -> io::Result<()> { + let args = Args::parse(); + 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(); + + // Validate paths + check_paths_differ(&args.fasta, &args.output)?; + 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() + }); + check_paths_differ(&args.fasta, &vp)?; + check_paths_differ(&args.output, &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 { + return Err(io::Error::new(io::ErrorKind::InvalidData, + format!("Input file '{}' 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[..]; + + // Index records + let (records, seq_length, layout) = index_fasta(data)?; + let num_samples = records.len(); + + eprintln!("[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); + let t1 = start.elapsed().as_secs_f64(); + + 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.", + 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); + + // 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)?; } + w.flush()?; + return Ok(()); + } - // Step 2: Extract variable positions and write to output file - println!("Starting Step 2: Extracting variable positions and writing to output..."); - extract_and_write_variables( - &input_filename, - &output_filename, - &variable_positions_info, - total_sequences, - &system, - ) - .expect("Failed to extract and write variable positions"); - println!("Variable positions alignment written to {}", output_filename); - - // Step 3: Generate VCF if requested - if generate_vcf { - println!("Generating VCF file..."); - generate_vcf_file( - &variable_positions_info, - &vcf_output_filename, - &sample_names, - ) - .expect("Failed to generate VCF file"); - println!("VCF file generated: {}", vcf_output_filename); + // VCF size guard + if do_vcf { + let geno_bytes = num_var.saturating_mul(num_samples); + if geno_bytes > MAX_VCF_GENO_BYTES { + return Err(io::Error::new(io::ErrorKind::InvalidInput, + format!("VCF genotype matrix would require {} GB ({} vars Γ— {} samples). \ + Use without --vcf or reduce input.", + geno_bytes / 1_000_000_000, num_var, num_samples))); } - }); + } + // Pass 2: extract variable sites + let ep = ExtractParams { + records: &records, output: &args.output, + collect_vcf: do_vcf, lookup: &lookup, upper: &upper, layout, + }; + let vcf_geno = pass2_extract(data, &mut var_positions, &ep)?; + + // 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); + } + + eprintln!("[snpick] Done in {:.2}s. {} vars from {} seqs Γ— {} pos.", + start.elapsed().as_secs_f64(), num_var, num_samples, seq_length); Ok(()) } -/// Structure to store information about variable positions -struct VariablePositionInfo { - position: usize, - reference_base: u8, - alternate_bases: HashSet, - genotypes: Vec, // Bases at this position for each sample +fn main() { + if let Err(e) = run() { + eprintln!("[snpick] Error: {}", e); + std::process::exit(1); + } } -/// Step 1: Identify variable positions and extract individual genotypes -fn identify_variable_positions( - input_filename: &str, - _system: &Arc>, - include_gaps: bool, -) -> io::Result<(Vec, usize, Vec)> { - // Open the input FASTA file - let input_file = File::open(input_filename)?; - let reader = BufReader::with_capacity(16 * 1024 * 1024, input_file); - let fasta_reader = fasta::Reader::new(reader); - - // Read all sequences and store names and sequences - let mut sequences = Vec::new(); - let mut sample_names = Vec::new(); - - for result in fasta_reader.records() { - let record = result?; - let id = record.id().to_string(); - let seq = record.seq().to_owned(); - sample_names.push(id); - sequences.push(seq); +// ============================================================================= +// Tests +// ============================================================================= + +#[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; + + fn tmp(name: &str, c: &str) -> String { + let p = format!("/tmp/snpick_t_{}.fa", name); + std::fs::write(&p, c).unwrap(); p + } + fn setup(path: &str) -> Mmap { + let f = File::open(path).unwrap(); + unsafe { Mmap::map(&f).unwrap() } } - let total_sequences = sequences.len(); + #[test] fn test_lookup() { + let lk = build_lookup(false); + assert_eq!(lk[b'A' as usize], BIT_A); + assert_eq!(lk[b'N' as usize], 0); + assert_eq!(build_lookup(true)[b'-' as usize], BIT_GAP); + } - if total_sequences == 0 { - eprintln!("The input FASTA file is empty."); - return Ok((Vec::new(), 0, Vec::new())); + #[test] fn test_index() { + let p = tmp("idx3", ">s1\nATGC\n>s2\nATCC\n"); + let m = setup(&p); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + assert_eq!(recs.len(), 2); + assert_eq!(sl, 4); + assert!(layout.single_line); + assert_eq!(recs[0].id, b"s1"); + assert_eq!(recs[1].id, b"s2"); + std::fs::remove_file(&p).ok(); } - let seq_length = sequences[0].len(); + #[test] fn test_pass1() { + let p = tmp("p1g", ">s1\nATGC\n>s2\nATCC\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); + assert_eq!(bm[2], BIT_G | BIT_C); + std::fs::remove_file(&p).ok(); + } - // Verify that all sequences have the same length - for seq in &sequences { - if seq.len() != seq_length { - eprintln!("All sequences must have the same length."); - return Ok((Vec::new(), 0, Vec::new())); - } + #[test] fn test_e2e() { + let p = tmp("e2eg", ">ref\nATGCATGC\n>s1\nATGTATGC\n>s2\nACGCATGC\n"); + let o = "/tmp/snpick_t_e2eg_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); + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + 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], "TC"); assert_eq!(l[3], "TT"); assert_eq!(l[5], "CC"); + std::fs::remove_file(&p).ok(); std::fs::remove_file(o).ok(); } - println!( - "[{}] Processing {} sequences of length {}.", - Local::now().format("%Y-%m-%d %H:%M:%S"), - total_sequences, - seq_length - ); - - // Initialize a progress spinner - let pb = ProgressBar::new_spinner(); - pb.set_style( - ProgressStyle::default_spinner() - .template("{spinner:.green} Processing positions: {pos}/{len}") - .unwrap(), - ); - pb.enable_steady_tick(Duration::from_millis(100)); - - // Atomic counter for processed positions - let pos_counter = AtomicUsize::new(0); - - // Identify variable positions and extract genotypes - let variable_positions_info: Vec = (0..seq_length) - .into_par_iter() - .filter_map(|pos| { - let mut counts = HashMap::new(); - let mut genotypes = Vec::with_capacity(total_sequences); - - // Collect bases at this position for all samples - for seq in &sequences { - let nuc = seq[pos]; - genotypes.push(nuc); - - if nucleotide_to_bit(nuc, include_gaps).is_some() { - counts.entry(nuc).and_modify(|c| *c += 1).or_insert(1); - } - } - - // Count types of nucleotides A, C, G, T (and gap if included) - let nucleotide_types = counts.len(); - - // Update the progress spinner - let current = pos_counter.fetch_add(1, Ordering::SeqCst) + 1; - pb.set_position(current as u64); - pb.set_length(seq_length as u64); - - if nucleotide_types > 1 { - // Determine the reference base (most frequent) - let mut max_count = 0; - let mut reference_base = b'N'; - for (&nuc, &count) in &counts { - if count > max_count { - max_count = count; - reference_base = nuc; - } - } - - // Alternate bases - let alternate_bases: HashSet = counts - .keys() - .cloned() - .filter(|&nuc| nuc != reference_base) - .collect(); - - Some(VariablePositionInfo { - position: pos, - reference_base, - alternate_bases, - genotypes, - }) - } else { - None - } - }) - .collect(); - - pb.finish_with_message("Position processing completed."); - - println!( - "[{}] Variable position identification completed.", - Local::now().format("%Y-%m-%d %H:%M:%S") - ); - println!("Total sequences processed: {}", total_sequences); - - Ok((variable_positions_info, total_sequences, sample_names)) -} + #[test] fn test_vcf() { + let p = tmp("vcfg", ">ref\nATGC\n>s1\nCTGC\n>s2\nGTGC\n>s3\nTTGC\n"); + let fo = "/tmp/snpick_t_vcfg_out.fa"; let vo = "/tmp/snpick_t_vcfg.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 }; + let g = pass2_extract(&m, &mut v, &ep).unwrap().unwrap(); + write_vcf(&g, recs.len(), &v, vo, &recs, sl).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"); assert_eq!(f[4], "C,G,T"); + assert_eq!(f[7], "NS=4"); assert_eq!(f[9], "0"); assert_eq!(f[12], "3"); + std::fs::remove_file(&p).ok(); std::fs::remove_file(fo).ok(); std::fs::remove_file(vo).ok(); + } -/// Step 2: Extract variable positions and write to the output FASTA file -fn extract_and_write_variables( - input_filename: &str, - output_filename: &str, - variable_positions_info: &[VariablePositionInfo], - total_sequences: usize, - system: &Arc>, -) -> io::Result<()> { - // Open the input FASTA file again for the second step - let input_file = File::open(input_filename)?; - let reader = BufReader::with_capacity(16 * 1024 * 1024, input_file); - let fasta_reader = fasta::Reader::new(reader); - - // Open the output FASTA file for writing - let output_file = File::create(output_filename)?; - let writer = BufWriter::with_capacity(16 * 1024 * 1024, output_file); - let writer = Arc::new(Mutex::new(writer)); - - // Get the variable positions - let variable_positions: Vec = variable_positions_info.iter().map(|info| info.position).collect(); - - // Initialize a progress spinner for Step 2 - let pb_write = ProgressBar::new_spinner(); - pb_write.set_style( - ProgressStyle::default_spinner() - .template("{spinner:.green} Writing sequences: {pos}/{len}") - .unwrap(), - ); - pb_write.enable_steady_tick(Duration::from_millis(100)); - - // Atomic counter for written sequences - let seq_counter = AtomicUsize::new(0); - - // Process and write sequences in parallel - fasta_reader - .records() - .par_bridge() - .try_for_each(|result| -> io::Result<()> { - let record = result?; - let seq = record.seq(); - - // Extract variable nucleotides based on variable_positions - let var_seq: Vec = variable_positions.iter().map(|&pos| seq[pos]).collect(); - - // Convert to String - let var_seq_str = String::from_utf8_lossy(&var_seq).to_string(); - - // Safely write to the output file - { - let mut writer = writer.lock().unwrap(); - writeln!(writer, ">{}", record.id())?; - writeln!(writer, "{}", var_seq_str)?; - } - - // Update the progress spinner - let current = seq_counter.fetch_add(1, Ordering::SeqCst) + 1; - pb_write.set_position(current as u64); - pb_write.set_length(total_sequences as u64); - - // Report RAM usage every 100,000 sequences - if current % 100_000 == 0 { - if let Ok(mut sys) = system.lock() { - sys.refresh_memory(); - let total_memory = sys.total_memory(); - let used_memory = sys.used_memory(); - println!( - "[{}] Written {} sequences. RAM usage: {} KB used / {} KB total.", - Local::now().format("%Y-%m-%d %H:%M:%S"), - current, - used_memory, - total_memory - ); - } - } - - Ok(()) - })?; - - // Finish the progress spinner for Step 2 - pb_write.finish_with_message(format!( - "Completed: {} sequences written.", - total_sequences - )); - - println!( - "[{}] Total sequences written: {}", - Local::now().format("%Y-%m-%d %H:%M:%S"), - total_sequences - ); + #[test] fn test_ambig() { + let p = tmp("ambg", ">ref\nATGC\n>s1\nCTGC\n>s2\nNTGC\n"); + let fo = "/tmp/snpick_t_ambg_out.fa"; let vo = "/tmp/snpick_t_ambg.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 }; + let g = pass2_extract(&m, &mut v, &ep).unwrap().unwrap(); + write_vcf(&g, recs.len(), &v, vo, &recs, sl).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[7], "NS=2"); assert_eq!(f[11], "."); + std::fs::remove_file(&p).ok(); std::fs::remove_file(fo).ok(); std::fs::remove_file(vo).ok(); + } - 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"; + 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 }; + pass2_extract(&m, &mut v, &ep).unwrap(); + let c = std::fs::read_to_string(o).unwrap(); + assert!(c.contains(">s1 some description")); + assert!(c.contains(">s2 another desc")); + std::fs::remove_file(&p).ok(); std::fs::remove_file(o).ok(); + } -/// Step 3: Generate the VCF file with actual bases, including ambiguous bases and codons -fn generate_vcf_file( - variable_positions_info: &[VariablePositionInfo], - vcf_output_filename: &str, - sample_names: &[String], -) -> io::Result<()> { - // Open the output VCF file for writing - let output_file = File::create(vcf_output_filename)?; - let mut writer = BufWriter::new(output_file); - - // Write the VCF header - writeln!(writer, "##fileformat=VCFv4.2")?; - writeln!(writer, "##source=snpick")?; - writeln!(writer, "##reference=.")?; - writeln!(writer, "##contig=", variable_positions_info.len())?; - writeln!( - writer, - "##INFO=" - )?; - writeln!( - writer, - "##FORMAT=" - )?; - - // Write the header line with sample names - write!(writer, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT")?; - for sample_name in sample_names { - write!(writer, "\t{}", sample_name)?; + #[test] fn test_gaps() { + let p = tmp("gapg", ">ref\nATGC\n>s1\nA-GC\n"); + let m = setup(&p); + let lk_no = build_lookup(false); + let lk_yes = build_lookup(true); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + let bm1 = pass1_scan(&m, &recs, sl, layout, &lk_no); + let rs = get_ref_seq(&m, &recs[0], sl, layout); + let (v1, _) = analyze(&bm1, &rs, &lk_no, false); + assert!(v1.is_empty()); + let bm2 = pass1_scan(&m, &recs, sl, layout, &lk_yes); + let (v2, _) = analyze(&bm2, &rs, &lk_yes, true); + assert_eq!(v2.len(), 1); + std::fs::remove_file(&p).ok(); } - writeln!(writer)?; - - // Generate VCF entries - for info in variable_positions_info { - let chrom = "1"; // Change this if you have chromosome information - let pos = info.position + 1; // Positions in VCF are 1-based - let id = "."; - let ref_base = info.reference_base as char; - let mut alt_bases: Vec = info.alternate_bases.iter().map(|&b| b as char).collect(); - alt_bases.sort(); - let alt = alt_bases.iter().collect::().replace('-', "."); - let qual = "."; - let filter = "PASS"; - let info_field = format!("NS={}", sample_names.len()); - let format_field = "BASE"; // Using a custom field - - // Generate genotypes for each sample, using the actual bases - let genotypes: Vec = info - .genotypes - .iter() - .map(|&genotype_base| { - let base_char = genotype_base as char; - base_char.to_string() - }) - .collect(); - - // Write the VCF line - write!( - writer, - "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}", - chrom, pos, id, ref_base, alt, qual, filter, info_field, format_field - )?; - for gt in genotypes { - write!(writer, "\t{}", gt)?; - } - writeln!(writer)?; + + #[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()); + std::fs::remove_file(&p).ok(); } - // Ensure all data is written correctly - writer.flush()?; + #[test] fn test_multiline_pass1() { + let p = tmp("mlg1", ">s1\nAT\nGC\n>s2\nAT\nCC\n"); + 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(); + } - Ok(()) + #[test] fn test_multiline_pass2() { + let p = tmp("mlg2", ">s1\nAT\nGC\n>s2\nAT\nCC\n"); + let o = "/tmp/snpick_t_mlg2_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); + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + 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], "G"); assert_eq!(l[3], "C"); + std::fs::remove_file(&p).ok(); std::fs::remove_file(o).ok(); + } + + #[test] fn test_multiline_vcf() { + let p = tmp("mlg3", ">ref\nAT\nGC\n>s1\nAT\nCC\n>s2\nCT\nGC\n"); + let fo = "/tmp/snpick_t_mlg3_out.fa"; let vo = "/tmp/snpick_t_mlg3.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 }; + let g = pass2_extract(&m, &mut v, &ep).unwrap().unwrap(); + write_vcf(&g, recs.len(), &v, vo, &recs, sl).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"); + let vc = std::fs::read_to_string(vo).unwrap(); + let dl: Vec<&str> = vc.lines().filter(|l| !l.starts_with('#')).collect(); + assert_eq!(dl.len(), 2); + std::fs::remove_file(&p).ok(); std::fs::remove_file(fo).ok(); std::fs::remove_file(vo).ok(); + } + + #[test] fn test_multiline_wide() { + let p = tmp("mlgw", ">s1\nATG\nCAT\nG\n>s2\nATG\nCCT\nG\n"); + let o = "/tmp/snpick_t_mlgw_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, 7); 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, 4); + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + 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], "C"); + std::fs::remove_file(&p).ok(); std::fs::remove_file(o).ok(); + } + + #[test] fn test_mixed_wrapping() { + let p = tmp("mxwg", ">s1\nATGCATGC\n>s2\nATGC\nATCC\n"); + let o = "/tmp/snpick_t_mxwg_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, 8); 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, 6); + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + 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], "G"); assert_eq!(l[3], "C"); + std::fs::remove_file(&p).ok(); std::fs::remove_file(o).ok(); + } + + #[test] fn test_crlf_multiline() { + // Windows line endings (\r\n) with multi-line wrapping + let p = tmp("crlfml", ""); + std::fs::write(&p, b">s1\r\nAT\r\nGC\r\n>s2\r\nAT\r\nCC\r\n").unwrap(); + 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, 2); + let o = "/tmp/snpick_t_crlfml_out.fa"; + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + 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], "G"); assert_eq!(l[3], "C"); + 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"); + let m = setup(&p); + let lk = build_lookup(false); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + assert_eq!(recs.len(), 1); + 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!(v.is_empty()); + assert_eq!(sc.constant.total(), 4); + assert_eq!(sc.variable, 0); + std::fs::remove_file(&p).ok(); + } + + #[test] fn test_no_trailing_newline() { + let p = tmp("noeof", ""); + std::fs::write(&p, b">s1\nATGC\n>s2\nATCC").unwrap(); + let m = setup(&p); + let lk = build_lookup(false); + let up = build_upper(); + let (recs, sl, layout) = index_fasta(&m).unwrap(); + assert_eq!(recs.len(), 2); + 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); + let o = "/tmp/snpick_t_noeof_out.fa"; + let ep = ExtractParams { records: &recs, output: o, collect_vcf: false, lookup: &lk, upper: &up, layout }; + 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], "G"); assert_eq!(l[3], "C"); + std::fs::remove_file(&p).ok(); std::fs::remove_file(o).ok(); + } + + #[test] fn test_empty() { + let p = tmp("empg", ""); + let m = setup(&p); + assert!(index_fasta(&m).is_err()); + std::fs::remove_file(&p).ok(); + } } diff --git a/src/scan.rs b/src/scan.rs new file mode 100644 index 0000000..2e18df9 --- /dev/null +++ b/src/scan.rs @@ -0,0 +1,128 @@ +//! Pass 1: bitmask scan and site classification. +//! +//! Builds a per-position bitmask by OR-ing each sequence's nucleotide flags, +//! then classifies positions as variable (>1 allele), constant, or ambiguous. + +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). +#[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; + } + std::hint::black_box(sum); +} + +/// Pass 1: build bitmask of observed nucleotides at each position. +/// +/// Iterates all sequences, OR-ing each base's lookup value into the bitmask. +/// Prefaults mmap pages first, then scans sequentially per sequence. +/// For multi-line FASTA, scans byte-by-byte skipping newlines. +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; + } + } + } + + bitmask +} + +/// 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], +) { + 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]; + } + } + } else { + for rec in records { + 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; } + bitmask[i] |= lookup[b as usize]; + i += 1; + } + } + } +} + +/// Classify positions into variable, constant, or ambiguous-only. +pub fn analyze( + bitmask: &[u8], ref_seq: &[u8], lookup: &[u8; 256], include_gaps: bool, +) -> (Vec, SiteCounts) { + let mut vars = Vec::new(); + let mut cs = ConstantSiteCounts { a: 0, c: 0, g: 0, t: 0 }; + let mut ambiguous = 0usize; + + for (pos, &bits) in bitmask.iter().enumerate() { + let ones = bits.count_ones(); + if ones > 1 { + let rb = ref_seq[pos].to_ascii_uppercase(); + let ref_base = if lookup[rb as usize] != 0 { rb } else { bits_to_bases(bits, include_gaps)[0] }; + let alt_bases: Vec = bits_to_bases(bits, include_gaps) + .into_iter().filter(|&b| b != ref_base).collect(); + vars.push(VariablePosition { index: pos, ref_base, alt_bases, ns: 0 }); + } else if ones == 1 { + if bits & BIT_A != 0 { cs.a += 1; } + 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; + } + } + + let num_variable = vars.len(); + (vars, SiteCounts { constant: cs, variable: num_variable, ambiguous }) +} diff --git a/src/types.rs b/src/types.rs new file mode 100644 index 0000000..a8dab5a --- /dev/null +++ b/src/types.rs @@ -0,0 +1,92 @@ +//! Shared types, constants, and lookup tables for snpick. +//! +//! This module defines the core data structures (variable positions, site counts) +//! and the bitmask-based nucleotide classification used throughout the pipeline. + +/// Bitmask flags for nucleotide classification. +pub const BIT_A: u8 = 0b00001; +pub const BIT_C: u8 = 0b00010; +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 VCF genotype matrix size in bytes. +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; + +/// 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. +#[derive(Clone, Copy)] +pub struct SeqLayout { + pub single_line: bool, +} + +/// A variable position detected in the alignment. +pub struct VariablePosition { + pub index: usize, + pub ref_base: u8, + pub alt_bases: Vec, + pub ns: usize, +} + +/// Counts of constant sites by nucleotide. +pub struct ConstantSiteCounts { + pub a: usize, + pub c: usize, + pub g: usize, + pub t: usize, +} + +impl std::fmt::Display for ConstantSiteCounts { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "A:{} C:{} G:{} T:{}", self.a, self.c, self.g, self.t) + } +} + +impl ConstantSiteCounts { + pub fn total(&self) -> usize { self.a + self.c + self.g + self.t } + pub fn fconst(&self) -> String { format!("{},{},{},{}", self.a, self.c, self.g, self.t) } +} + +/// Summary of site classification across the alignment. +pub struct SiteCounts { + pub constant: ConstantSiteCounts, + pub variable: usize, + pub ambiguous: usize, +} + +/// Build nucleotide β†’ bitmask lookup table. +pub fn build_lookup(include_gaps: bool) -> [u8; 256] { + let mut t = [0u8; 256]; + t[b'A' as usize] = BIT_A; t[b'a' as usize] = BIT_A; + t[b'C' as usize] = BIT_C; t[b'c' as usize] = BIT_C; + t[b'G' as usize] = BIT_G; t[b'g' as usize] = BIT_G; + t[b'T' as usize] = BIT_T; t[b't' as usize] = BIT_T; + if include_gaps { t[b'-' as usize] = BIT_GAP; } + t +} + +/// Build lowercase β†’ uppercase lookup table. +pub fn build_upper() -> [u8; 256] { + let mut t = [0u8; 256]; + for (i, v) in t.iter_mut().enumerate() { *v = i as u8; } + for c in b'a'..=b'z' { t[c as usize] = c - 32; } + t +} + +/// Convert bitmask to sorted list of bases it represents. +pub fn bits_to_bases(bits: u8, include_gaps: bool) -> Vec { + let mut v = Vec::with_capacity(5); + if bits & BIT_A != 0 { v.push(b'A'); } + if bits & BIT_C != 0 { v.push(b'C'); } + if bits & BIT_G != 0 { v.push(b'G'); } + if bits & BIT_T != 0 { v.push(b'T'); } + if include_gaps && bits & BIT_GAP != 0 { v.push(b'-'); } + v +} diff --git a/src/vcf.rs b/src/vcf.rs new file mode 100644 index 0000000..e69b877 --- /dev/null +++ b/src/vcf.rs @@ -0,0 +1,64 @@ +//! VCF output writer. +//! +//! 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::fasta::FastaRecord; +use crate::types::VariablePosition; + +/// Write VCF output from genotype matrix and variable positions. +pub fn write_vcf( + vcf_geno: &[u8], num_samples: usize, var_positions: &[VariablePosition], + vcf_path: &str, records: &[FastaRecord], seq_length: usize, +) -> 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); + + // 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, "##INFO=")?; + writeln!(w, "##FORMAT=")?; + write!(w, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT")?; + for rec in records { + write!(w, "\t")?; + w.write_all(rec.id)?; + } + writeln!(w)?; + + // Data rows + let mut lut = [255u8; 256]; + 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(","); + + // 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)?; + + let row = 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)?; } + } + writeln!(w)?; + + // Reset LUT entries + lut[vp.ref_base as usize] = 255; + for &ab in &vp.alt_bases { lut[ab as usize] = 255; } + } + + w.flush() +}