diff --git a/Cargo.lock b/Cargo.lock index baa3b38..3d3d3b9 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -8,6 +8,12 @@ version = "1.0.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f26201604c87b1e01bd3d98f8d5d9a8fcbb815e8cedb41ffccbeb4bf593a35fe" +[[package]] +name = "adler2" +version = "2.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "512761e0bb2578dd7380c6baaa0f4ce03e84f95e960231d1dec8bf4d7d6e2627" + [[package]] name = "aho-corasick" version = "1.1.3" @@ -34,57 +40,58 @@ dependencies = [ [[package]] name = "anstream" -version = "0.6.13" +version = "0.6.15" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d96bd03f33fe50a863e394ee9718a706f988b9079b20c3784fb726e7678b62fb" +checksum = "64e15c1ab1f89faffbf04a634d5e1962e9074f2741eef6d97f3c4e322426d526" dependencies = [ "anstyle", "anstyle-parse", "anstyle-query", "anstyle-wincon", "colorchoice", + "is_terminal_polyfill", "utf8parse", ] [[package]] name = "anstyle" -version = "1.0.6" +version = "1.0.8" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8901269c6307e8d93993578286ac0edf7f195079ffff5ebdeea6a59ffb7e36bc" +checksum = "1bec1de6f59aedf83baf9ff929c98f2ad654b97c9510f4e70cf6f661d49fd5b1" [[package]] name = "anstyle-parse" -version = "0.2.3" +version = "0.2.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c75ac65da39e5fe5ab759307499ddad880d724eed2f6ce5b5e8a26f4f387928c" +checksum = "eb47de1e80c2b463c735db5b217a0ddc39d612e7ac9e2e96a5aed1f57616c1cb" dependencies = [ "utf8parse", ] [[package]] name = "anstyle-query" -version = "1.0.2" +version = "1.1.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e28923312444cdd728e4738b3f9c9cac739500909bb3d3c94b43551b16517648" +checksum = "6d36fc52c7f6c869915e99412912f22093507da8d9e942ceaf66fe4b7c14422a" dependencies = [ - "windows-sys", + "windows-sys 0.52.0", ] [[package]] name = "anstyle-wincon" -version = "3.0.2" +version = "3.0.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1cd54b81ec8d6180e24654d0b371ad22fc3dd083b6ff8ba325b72e00c87660a7" +checksum = "5bf74e1b6e971609db8ca7a9ce79fd5768ab6ae46441c572e46cf596f59e57f8" dependencies = [ "anstyle", - "windows-sys", + "windows-sys 0.52.0", ] [[package]] name = "anyhow" -version = "1.0.82" +version = "1.0.86" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f538837af36e6f6a9be0faa67f9a314f8119e4e4b5867c6ab40ed60360142519" +checksum = "b3d1d046238990b9cf5bcde22a3fb3584ee5cf65fb2765f454ed428c7a0063da" [[package]] name = "approx" @@ -97,21 +104,21 @@ dependencies = [ [[package]] name = "arrayref" -version = "0.3.7" +version = "0.3.8" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6b4930d2cb77ce62f89ee5d5289b4ac049559b1c45539271f5ed4fdc7db34545" +checksum = "9d151e35f61089500b617991b791fc8bfd237ae50cd5950803758a179b41e67a" [[package]] name = "arrayvec" -version = "0.7.4" +version = "0.7.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "96d30a06541fbafbc7f82ed10c06164cfbd2c401138f6addd8404629c4b16711" +checksum = "7c02d123df017efcdfbd739ef81735b36c5ba83ec3c59c80a9d7ecc718f92e50" [[package]] name = "autocfg" -version = "1.2.0" +version = "1.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f1fdabc7756949593fe60f30ec81974b613357de856987752631dea1e3394c80" +checksum = "0c4b4d0bd25bd0b74681c0ad21497610ce1b7c91b1022cd21c80c6fbdd9476b0" [[package]] name = "base64" @@ -152,7 +159,7 @@ dependencies = [ "serde_derive", "statrs", "strum", - "strum_macros", + "strum_macros 0.25.3", "thiserror", "triple_accel", "vec_map", @@ -160,14 +167,14 @@ dependencies = [ [[package]] name = "bio-types" -version = "1.0.1" +version = "1.0.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9d45749b87f21808051025e9bf714d14ff4627f9d8ca967eade6946ea769aa4a" +checksum = "f4dcf54f8b7f51450207d54780bab09c05f30b8b0caa991545082842e466ad7e" dependencies = [ - "derive-new", + "derive-new 0.6.0", "lazy_static", "regex", - "strum_macros", + "strum_macros 0.26.4", "thiserror", ] @@ -194,9 +201,9 @@ checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" [[package]] name = "bitflags" -version = "2.5.0" +version = "2.6.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cf4b9d6a944f767f8e5e0db018570623c85f3d925ac718db4e06d0187adb21c1" +checksum = "b048fb63fd8b5923fc5aa7b340d8e156aec7ec02f0c78fa8a6ddc2613f6f71de" [[package]] name = "bumpalo" @@ -216,15 +223,15 @@ dependencies = [ [[package]] name = "bytecount" -version = "0.6.7" +version = "0.6.8" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e1e5f035d16fc623ae5f74981db80a439803888314e3a555fd6f04acd51a3205" +checksum = "5ce89b21cab1437276d2650d57e971f9d548a2d9037cc231abdc0562b97498ce" [[package]] name = "bytemuck" -version = "1.15.0" +version = "1.17.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5d6d68c57235a3a081186990eca2867354726650f42f7516ca50c28d6281fd15" +checksum = "6fd4c6dcc3b0aea2f5c0b4b82c2b15fe39ddbc76041a310848f4706edf76bb31" [[package]] name = "byteorder" @@ -234,12 +241,13 @@ checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" [[package]] name = "cc" -version = "1.0.94" +version = "1.1.13" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "17f6e324229dc011159fcc089755d1e2e216a90d43a7dea6853ca740b84f35e7" +checksum = "72db2f7947ecee9b03b510377e8bb9077afa27176fdbff55c51027e976fdcc48" dependencies = [ "jobserver", "libc", + "shlex", ] [[package]] @@ -264,9 +272,9 @@ dependencies = [ [[package]] name = "clap" -version = "4.5.4" +version = "4.5.16" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "90bc066a67923782aa8515dbaea16946c5bcc5addbd668bb80af688e53e548a0" +checksum = "ed6719fffa43d0d87e5fd8caeab59be1554fb028cd30edc88fc4369b17971019" dependencies = [ "clap_builder", "clap_derive", @@ -274,9 +282,9 @@ dependencies = [ [[package]] name = "clap_builder" -version = "4.5.2" +version = "4.5.15" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ae129e2e766ae0ec03484e609954119f123cc1fe650337e155d03b022f24f7b4" +checksum = "216aec2b177652e3846684cbfe25c9964d18ec45234f0f5da5157b207ed1aab6" dependencies = [ "anstream", "anstyle", @@ -286,27 +294,27 @@ dependencies = [ [[package]] name = "clap_derive" -version = "4.5.4" +version = "4.5.13" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "528131438037fd55894f62d6e9f068b8f45ac57ffa77517819645d10aed04f64" +checksum = "501d359d5f3dcaf6ecdeee48833ae73ec6e42723a1e52419c79abf9507eec0a0" dependencies = [ "heck 0.5.0", "proc-macro2", "quote", - "syn 2.0.59", + "syn 2.0.75", ] [[package]] name = "clap_lex" -version = "0.7.0" +version = "0.7.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "98cc8fbded0c607b7ba9dd60cd98df59af97e84d24e49c8557331cfc26d301ce" +checksum = "1462739cb27611015575c0c11df5df7601141071f07518d56fcc1be504cbec97" [[package]] name = "cmake" -version = "0.1.50" +version = "0.1.51" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a31c789563b815f77f4250caee12365734369f942439b7defd71e18a48197130" +checksum = "fb1e43aa7fd152b1f968787f7dbcdeb306d1867ff373c69955211876c053f91a" dependencies = [ "cc", ] @@ -319,30 +327,30 @@ checksum = "3d7b894f5411737b7867f4827955924d7c254fc9f4d91a6aad6b097804b1018b" [[package]] name = "colorchoice" -version = "1.0.0" +version = "1.0.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "acbf1af155f9b9ef647e42cdc158db4b64a1b61f743629225fde6f3e0be2a7c7" +checksum = "d3fd119d74b830634cea2a0f58bbd0d54540518a14397557951e79340abc28c0" [[package]] name = "core-foundation-sys" -version = "0.8.6" +version = "0.8.7" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "06ea2b9bc92be3c2baa9334a323ebca2d6f074ff852cd1d7b11064035cd3868f" +checksum = "773648b94d0e5d620f64f280777445740e61fe701025087ec8b57f45c791888b" [[package]] name = "crc32fast" -version = "1.4.0" +version = "1.4.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b3855a8a784b474f333699ef2bbca9db2c4a1f6d9088a90a2d25b1eb53111eaa" +checksum = "a97769d94ddab943e4510d138150169a2758b5ef3eb191a9ee688de3e23ef7b3" dependencies = [ "cfg-if", ] [[package]] name = "crossbeam-channel" -version = "0.5.12" +version = "0.5.13" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ab3db02a9c5b5121e1e42fbdb1aeb65f5e02624cc58c43f2884c6ccac0b82f95" +checksum = "33480d6946193aa8033910124896ca395333cae7e2d1113d1fef6c3272217df2" dependencies = [ "crossbeam-utils", ] @@ -368,9 +376,9 @@ dependencies = [ [[package]] name = "crossbeam-utils" -version = "0.8.19" +version = "0.8.20" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "248e3bacc7dc6baa3b21e405ee045c3047101a49145e7e9eca583ab4c2ca5345" +checksum = "22ec99545bb0ed0ea7bb9b8e1e9122ea386ff8a48c0922e43f36d45ab09e0e80" [[package]] name = "csv" @@ -425,6 +433,17 @@ dependencies = [ "syn 1.0.109", ] +[[package]] +name = "derive-new" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d150dea618e920167e5973d70ae6ece4385b7164e0d799fe7c122dd0a5d912ad" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.75", +] + [[package]] name = "editdistancek" version = "1.0.2" @@ -433,9 +452,9 @@ checksum = "3e02df23d5b1c6f9e69fa603b890378123b93073df998a21e6e33b9db0a32613" [[package]] name = "either" -version = "1.11.0" +version = "1.13.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a47c1c47d2f5964e29c61246e81db715514cd532db6b5116a25ea3c03d6780a2" +checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0" [[package]] name = "enum-map" @@ -454,7 +473,7 @@ checksum = "f282cfdfe92516eb26c2af8589c274c7c17681f5ecc03c18255fe741c6aa64eb" dependencies = [ "proc-macro2", "quote", - "syn 2.0.59", + "syn 2.0.75", ] [[package]] @@ -478,19 +497,19 @@ checksum = "5443807d6dff69373d433ab9ef5378ad8df50ca6298caf15de6e52e24aaf54d5" [[package]] name = "errno" -version = "0.3.8" +version = "0.3.9" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a258e46cdc063eb8519c00b9fc845fc47bcfca4130e2f08e88665ceda8474245" +checksum = "534c5cf6194dfab3db3242765c03bbe257cf92f22b38f6bc0c58d59108a820ba" dependencies = [ "libc", - "windows-sys", + "windows-sys 0.52.0", ] [[package]] name = "fastrand" -version = "2.0.2" +version = "2.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "658bd65b1cf4c852a3cc96f18a8ce7b5640f6b703f905c7d74532294c2a63984" +checksum = "9fc0510504f03c51ada170672ac806f1f105a88aa97a5281117e1ddc3368e51a" [[package]] name = "fdeflate" @@ -515,12 +534,12 @@ checksum = "0ce7134b9999ecaf8bcd65542e436736ef32ddca1b3e06094cb6ec5755203b80" [[package]] name = "flate2" -version = "1.0.28" +version = "1.0.32" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "46303f565772937ffe1d394a4fac6f411c6013172fadde9dcdb1e147a086940e" +checksum = "9c0596c1eac1f9e04ed902702e9878208b336edc9d6fddc8a48387349bab3666" dependencies = [ "crc32fast", - "miniz_oxide", + "miniz_oxide 0.8.0", ] [[package]] @@ -531,11 +550,11 @@ checksum = "98de4bbd547a563b716d8dfa9aad1cb19bfab00f4fa09a6a4ed21dbcf44ce9c4" [[package]] name = "fontconfig-parser" -version = "0.5.6" +version = "0.5.7" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6a595cb550439a117696039dfc69830492058211b771a2a165379f2a1a53d84d" +checksum = "c1fcfcd44ca6e90c921fee9fa665d530b21ef1327a4c1a6c5250ea44b776ada7" dependencies = [ - "roxmltree 0.19.0", + "roxmltree 0.20.0", ] [[package]] @@ -581,9 +600,9 @@ dependencies = [ [[package]] name = "getrandom" -version = "0.2.14" +version = "0.2.15" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "94b22e06ecb0110981051723910cbf0b5f5e09a2062dd7663334ee79a9d1286c" +checksum = "c4567c8db10ae91089c99af84c68c38da3ec2f087c3f82960bcdbf3656b6f4d7" dependencies = [ "cfg-if", "libc", @@ -618,9 +637,9 @@ checksum = "d2fabcfbdc87f4758337ca535fb41a6d701b65693ce38287d856d1674551ec9b" [[package]] name = "hashbrown" -version = "0.14.3" +version = "0.14.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "290f1a1d9242c78d09ce40a5e87e7554ee637af1351968159f4952f028f75604" +checksum = "e5274423e17b7c9fc20b6e7e208532f9b19825d82dfd615708b70edd83df41f1" [[package]] name = "heck" @@ -636,15 +655,15 @@ checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" [[package]] name = "hermit-abi" -version = "0.3.9" +version = "0.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d231dfb89cfffdbc30e7fc41579ed6066ad03abda9e567ccafae602b97ec5024" +checksum = "fbf6a919d6cf397374f7dfeeea91d974c7c0a7221d0d0f4f20d859d329e53fcc" [[package]] name = "hts-sys" -version = "2.1.1" +version = "2.1.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "deebfb779c734d542e7f14c298597914b9b5425e4089aef482eacb5cab941915" +checksum = "e9f348d14cb4e50444e39fcd6b00302fe2ed2bc88094142f6278391d349a386d" dependencies = [ "cc", "fs-utils", @@ -720,9 +739,9 @@ checksum = "029d73f573d8e8d63e6d5020011d3255b28c3ba85d6cf870a07184ed23de9284" [[package]] name = "indexmap" -version = "2.2.6" +version = "2.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "168fb715dda47215e360912c096649d23d58bf392ac62f73919e831745e40f26" +checksum = "93ead53efc7ea8ed3cfb0c79fc8023fbb782a5432b52830b6518941cebe6505c" dependencies = [ "equivalent", "hashbrown", @@ -730,15 +749,21 @@ dependencies = [ [[package]] name = "is-terminal" -version = "0.4.12" +version = "0.4.13" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f23ff5ef2b80d608d61efee834934d862cd92461afc0560dedf493e4c033738b" +checksum = "261f68e344040fbd0edea105bef17c66edf46f984ddb1115b775ce31be948f4b" dependencies = [ "hermit-abi", "libc", - "windows-sys", + "windows-sys 0.52.0", ] +[[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.11.0" @@ -774,9 +799,9 @@ checksum = "49f1f14873335454500d59611f1cf4a4b0f786f9ac11f4312a78e4cf2566695b" [[package]] name = "jobserver" -version = "0.1.30" +version = "0.1.32" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "685a7d121ee3f65ae4fddd72b25a04bb36b6af81bc0828f7d5434c0fe60fa3a2" +checksum = "48d1dbcbbeb6a7fec7e059840aa538bd62aaccf972c7346c4d9d2059312853d0" dependencies = [ "libc", ] @@ -789,9 +814,9 @@ checksum = "f5d4a7da358eff58addd2877a45865158f0d78c911d43a5784ceb7bbf52833b0" [[package]] name = "js-sys" -version = "0.3.69" +version = "0.3.70" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "29c15563dc2726973df627357ce0c9ddddbea194836909d655df6a75d2cf296d" +checksum = "1868808506b929d7b0cfa8f75951347aa71bb21144b7791bae35d9bccfcfe37a" dependencies = [ "wasm-bindgen", ] @@ -813,15 +838,15 @@ dependencies = [ [[package]] name = "lazy_static" -version = "1.4.0" +version = "1.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" +checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe" [[package]] name = "libc" -version = "0.2.153" +version = "0.2.158" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9c198f91728a82281a64e1f4f9eeb25d82cb32a5de251c6bd1b5154d63a8e7bd" +checksum = "d8adc4bb1803a324070e64a98ae98f38934d91957a99cfb3a43dcbc01bc56439" [[package]] name = "libm" @@ -831,9 +856,9 @@ checksum = "4ec2a862134d2a7d32d7983ddcdd1c4923530833c9f2ea1a44fc5fa473989058" [[package]] name = "libz-sys" -version = "1.1.16" +version = "1.1.19" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5e143b5e666b2695d28f6bca6497720813f699c9602dd7f5cac91008b8ada7f9" +checksum = "fdc53a7799a7496ebc9fd29f31f7df80e83c9bda5299768af5f9e59eeea74647" dependencies = [ "cc", "cmake", @@ -850,21 +875,21 @@ checksum = "bfae20f6b19ad527b550c223fddc3077a547fc70cda94b9b566575423fd303ee" [[package]] name = "linux-raw-sys" -version = "0.4.13" +version = "0.4.14" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "01cda141df6706de531b6c46c3a33ecca755538219bd484262fa09410c13539c" +checksum = "78b3ae25bc7c8c38cec158d1f2757ee79e9b3740fbc7ccf0e59e4b08d793fa89" [[package]] name = "log" -version = "0.4.21" +version = "0.4.22" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "90ed8c1e510134f979dbc4f070f87d4313098b704861a105fe34231c70a3901c" +checksum = "a7a70ba024b9dc04c27ea2f0c0548feb474ec5c54bba33a7f72f873a39d07b24" [[package]] name = "matrixmultiply" -version = "0.3.8" +version = "0.3.9" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7574c1cf36da4798ab73da5b215bbf444f50718207754cb522201d78d1cd0ff2" +checksum = "9380b911e3e96d10c1f415da0876389aaf1b56759054eeb0de7df940c456ba1a" dependencies = [ "autocfg", "rawpointer", @@ -872,9 +897,9 @@ dependencies = [ [[package]] name = "memchr" -version = "2.7.2" +version = "2.7.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6c8640c5d730cb13ebd907d8d04b52f55ac9a2eec55b440c8892f40d56c76c1d" +checksum = "78ca9ab1a0babb1e7d5695e3530886289c18cf2f87ec19a575a0abdce112e3a3" [[package]] name = "memmap2" @@ -887,14 +912,23 @@ dependencies = [ [[package]] name = "miniz_oxide" -version = "0.7.2" +version = "0.7.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9d811f3e15f28568be3407c8e7fdb6514c1cda3cb30683f15b6a1a1dc4ea14a7" +checksum = "b8a240ddb74feaf34a79a7add65a741f3167852fba007066dcac1ca548d89c08" dependencies = [ "adler", "simd-adler32", ] +[[package]] +name = "miniz_oxide" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e2d80299ef12ff69b16a84bb182e3b9df68b5a91574d3d4fa6e41b65deec4df1" +dependencies = [ + "adler2", +] + [[package]] name = "multimap" version = "0.9.1" @@ -957,9 +991,9 @@ dependencies = [ [[package]] name = "num-complex" -version = "0.4.5" +version = "0.4.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "23c6602fda94a57c990fe0df199a035d83576b496aa29f4e634a8ac6004e68a6" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" dependencies = [ "num-traits", ] @@ -981,20 +1015,19 @@ dependencies = [ [[package]] name = "num-rational" -version = "0.4.1" +version = "0.4.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0" +checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" dependencies = [ - "autocfg", "num-integer", "num-traits", ] [[package]] name = "num-traits" -version = "0.2.18" +version = "0.2.19" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "da0df0e5185db44f69b44f26786fe401b6c293d1907744beaa7fa62b2e5a517a" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" dependencies = [ "autocfg", "libm", @@ -1026,15 +1059,15 @@ dependencies = [ [[package]] name = "paste" -version = "1.0.14" +version = "1.0.15" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "de3145af08024dea9fa9914f381a17b8fc6034dfb00f3a84013f7ff43f29ed4c" +checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" [[package]] name = "pdf-writer" -version = "0.9.2" +version = "0.9.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "644b654f2de28457bf1e25a4905a76a563d1128a33ce60cf042f721f6818feaf" +checksum = "24e9127455063c816e661caac9ecd9043ad2871f55be93014e6838a8ced2332b" dependencies = [ "bitflags 1.3.2", "itoa", @@ -1050,9 +1083,9 @@ checksum = "e3148f5046208a5d56bcfc03053e3ca6334e51da8dfb19b6cdc8b306fae3283e" [[package]] name = "petgraph" -version = "0.6.4" +version = "0.6.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e1d3afd2628e69da2be385eb6f2fd57c8ac7977ceeff6dc166ff1657b0e386a9" +checksum = "b4c5cc86750666a3ed20bdaf5ca2a0344f9c67674cae0515bec2da16fbaa47db" dependencies = [ "fixedbitset", "indexmap", @@ -1064,6 +1097,17 @@ version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5be167a7af36ee22fe3115051bc51f6e6c7054c9348e28deb4f49bd6f705a315" +[[package]] +name = "pipeplot" +version = "0.1.0" +dependencies = [ + "resvg", + "svg2pdf", + "tempfile", + "tiny-skia", + "usvg", +] + [[package]] name = "pkg-config" version = "0.3.30" @@ -1080,7 +1124,7 @@ dependencies = [ "crc32fast", "fdeflate", "flate2", - "miniz_oxide", + "miniz_oxide 0.7.4", ] [[package]] @@ -1091,15 +1135,18 @@ checksum = "439ee305def115ba05938db6eb1644ff94165c5ab5e9420d1c1bcedbba909391" [[package]] name = "ppv-lite86" -version = "0.2.17" +version = "0.2.20" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5b40af805b3121feab8a3c29f04d8ad262fa8e0561883e7653e024ae4479e6de" +checksum = "77957b295656769bb8ad2b6a6b09d897d94f05c41b069aede1fcdaa675eaea04" +dependencies = [ + "zerocopy", +] [[package]] name = "proc-macro2" -version = "1.0.80" +version = "1.0.86" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a56dea16b0a29e94408b9aa5e2940a4eedbd128a1ba20e8f7ae60fd3d465af0e" +checksum = "5e719e8df665df0d1c8fbfd238015744736151d4445ec0836b8e628aae103b77" dependencies = [ "unicode-ident", ] @@ -1193,9 +1240,9 @@ checksum = "3b42e27ef78c35d3998403c1d26f3efd9e135d3e5121b0a4845cc5cc27547f4f" [[package]] name = "regex" -version = "1.10.4" +version = "1.10.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c117dbdfde9c8308975b6a18d71f3f385c89461f7b3fb054288ecf2a2058ba4c" +checksum = "4219d74c6b67a3654a9fbebc4b419e22126d13d2f3c4a07ee0cb61ff79a79619" dependencies = [ "aho-corasick", "memchr", @@ -1205,9 +1252,9 @@ dependencies = [ [[package]] name = "regex-automata" -version = "0.4.6" +version = "0.4.7" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "86b83b8b9847f9bf95ef68afb0b8e6cdb80f498442f5179a29fad448fcc1eaea" +checksum = "38caf58cc5ef2fed281f89292ef23f6365465ed9a41b7a7754eb4e26496c92df" dependencies = [ "aho-corasick", "memchr", @@ -1216,9 +1263,9 @@ dependencies = [ [[package]] name = "regex-syntax" -version = "0.8.3" +version = "0.8.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "adad44e29e4c806119491a7f06f03de4d1af22c3a680dd47f1e6e179439d1f56" +checksum = "7a66a03ae7c801facd77a29370b4faec201768915ac14a721ba36f20bc9c209b" [[package]] name = "resvg" @@ -1239,9 +1286,9 @@ dependencies = [ [[package]] name = "rgb" -version = "0.8.37" +version = "0.8.48" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "05aaa8004b64fd573fc9d002f4e632d51ad4f026c2b5ba95fcb6c2f32c2c47d8" +checksum = "0f86ae463694029097b846d8f99fd5536740602ae00022c0c50c5600720b2f71" dependencies = [ "bytemuck", ] @@ -1257,9 +1304,9 @@ dependencies = [ [[package]] name = "roxmltree" -version = "0.19.0" +version = "0.20.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3cd14fd5e3b777a7422cca79358c57a8f6e3a703d9ac187448d0daf220c2407f" +checksum = "6c20b6793b5c2fa6553b250154b78d6d0db37e72700ae35fad9387a46f487c97" [[package]] name = "rust-htslib" @@ -1270,7 +1317,7 @@ dependencies = [ "bio-types", "byteorder", "custom_derive", - "derive-new", + "derive-new 0.5.9", "hts-sys", "ieee754", "lazy_static", @@ -1294,22 +1341,22 @@ dependencies = [ [[package]] name = "rustix" -version = "0.38.32" +version = "0.38.34" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "65e04861e65f21776e67888bfbea442b3642beaa0138fdb1dd7a84a52dffdb89" +checksum = "70dc5ec042f7a43c4a73241207cecc9873a06d45debb38b329f8541d85c2730f" dependencies = [ - "bitflags 2.5.0", + "bitflags 2.6.0", "errno", "libc", "linux-raw-sys", - "windows-sys", + "windows-sys 0.52.0", ] [[package]] name = "rustversion" -version = "1.0.15" +version = "1.0.17" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "80af6f9131f277a45a3fba6ce8e2258037bb0477a67e610d3c1fe046ab31de47" +checksum = "955d28af4278de8121b7ebeb796b6a45735dc01436d898801014aced2773a3d6" [[package]] name = "rustybuzz" @@ -1329,15 +1376,15 @@ dependencies = [ [[package]] name = "ryu" -version = "1.0.17" +version = "1.0.18" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e86697c916019a8588c99b5fac3cead74ec0b4b819707a682fd4d23fa0ce1ba1" +checksum = "f3cb5ba0dc43242ce17de99c180e96db90b235b8a9fdc9543c96d2209116bd9f" [[package]] name = "safe_arch" -version = "0.7.1" +version = "0.7.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f398075ce1e6a179b46f51bd88d0598b92b00d3551f1a2d4ac49e771b56ac354" +checksum = "c3460605018fdc9612bce72735cba0d27efbcd9904780d44c7e3a9948f96148a" dependencies = [ "bytemuck", ] @@ -1356,24 +1403,30 @@ checksum = "61697e0a1c7e512e84a621326239844a24d8207b4669b41bc18b32ea5cbf988b" [[package]] name = "serde" -version = "1.0.197" +version = "1.0.208" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3fb1c873e1b9b056a4dc4c0c198b24c3ffa059243875552b2bd0933b1aee4ce2" +checksum = "cff085d2cb684faa248efb494c39b68e522822ac0de72ccf08109abde717cfb2" dependencies = [ "serde_derive", ] [[package]] name = "serde_derive" -version = "1.0.197" +version = "1.0.208" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7eb0b34b42edc17f6b7cac84a52a1c5f0e1bb2227e997ca9011ea3dd34e8610b" +checksum = "24008e81ff7613ed8e5ba0cfaf24e2c2f1e5b8a0495711e44fcd4882fca62bcf" dependencies = [ "proc-macro2", "quote", - "syn 2.0.59", + "syn 2.0.75", ] +[[package]] +name = "shlex" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" + [[package]] name = "simba" version = "0.6.0" @@ -1425,9 +1478,9 @@ checksum = "3c5e1a9a646d36c3599cd173a41282daf47c44583ad367b8e6837255952e5c67" [[package]] name = "statrs" -version = "0.16.0" +version = "0.16.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2d08e5e1748192713cc281da8b16924fb46be7b0c2431854eadc785823e5696e" +checksum = "b35a062dbadac17a42e0fc64c27f419b25d6fae98572eb43c8814c9e873d7721" dependencies = [ "approx", "lazy_static", @@ -1467,7 +1520,20 @@ dependencies = [ "proc-macro2", "quote", "rustversion", - "syn 2.0.59", + "syn 2.0.75", +] + +[[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.75", ] [[package]] @@ -1477,7 +1543,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a81da66842e426278f20062cd249779565e13f9ab4bfe0ac9e94eb476bc3a0f3" dependencies = [ "image", - "miniz_oxide", + "miniz_oxide 0.7.4", "once_cell", "pdf-writer", "usvg", @@ -1506,9 +1572,9 @@ dependencies = [ [[package]] name = "syn" -version = "2.0.59" +version = "2.0.75" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4a6531ffc7b071655e4ce2e04bd464c4830bb585a61cabb96cf808f05172615a" +checksum = "f6af063034fc1935ede7be0122941bafa9bacb949334d090b77ca98b5817c7d9" dependencies = [ "proc-macro2", "quote", @@ -1517,14 +1583,15 @@ dependencies = [ [[package]] name = "tempfile" -version = "3.10.1" +version = "3.12.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "85b77fafb263dd9d05cbeac119526425676db3784113aa9295c88498cbf8bff1" +checksum = "04cbcdd0c794ebb0d4cf35e88edd2f7d2c4c3e9a5a6dab322839b321c6a87a64" dependencies = [ "cfg-if", "fastrand", + "once_cell", "rustix", - "windows-sys", + "windows-sys 0.59.0", ] [[package]] @@ -1538,22 +1605,22 @@ dependencies = [ [[package]] name = "thiserror" -version = "1.0.58" +version = "1.0.63" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "03468839009160513471e86a034bb2c5c0e4baae3b43f79ffc55c4a5427b3297" +checksum = "c0342370b38b6a11b6cc11d6a805569958d54cfa061a29969c3b5ce2ea405724" dependencies = [ "thiserror-impl", ] [[package]] name = "thiserror-impl" -version = "1.0.58" +version = "1.0.63" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c61f3ba182994efc43764a46c018c347bc492c79f024e705f46567b418f6d4f7" +checksum = "a4558b58466b9ad7ca0f102865eccc95938dca1a74a856f2b57b6629050da261" dependencies = [ "proc-macro2", "quote", - "syn 2.0.59", + "syn 2.0.75", ] [[package]] @@ -1617,9 +1684,9 @@ dependencies = [ [[package]] name = "tinyvec" -version = "1.6.0" +version = "1.8.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "87cc5ceb3875bb20c2890005a4e226a4651264a5c75edb2421b52861a0a0cb50" +checksum = "445e881f4f6d382d5f27c034e25eb92edd7c784ceab92a0937db7f2e9471b938" dependencies = [ "tinyvec_macros", ] @@ -1632,7 +1699,7 @@ checksum = "1f3ccbac311fea05f86f61904b462b55fb3df8837a366dfc601a0161d0532f20" [[package]] name = "trgt" -version = "1.2.0" +version = "1.3.0" dependencies = [ "arrayvec", "bio", @@ -1645,15 +1712,12 @@ dependencies = [ "kodama", "log", "once_cell", + "pipeplot", "rand", "rayon", - "resvg", "rust-htslib", "semver 1.0.23", - "svg2pdf", "tempfile", - "tiny-skia", - "usvg", "vergen", ] @@ -1710,9 +1774,9 @@ dependencies = [ [[package]] name = "unicode-properties" -version = "0.1.1" +version = "0.1.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e4259d9d4425d9f0661581b804cb85fe66a4c631cadd8f490d1c13a35d5d9291" +checksum = "52ea75f83c0137a9b98608359a5f1af8144876eb67bcb1ce837368e906a9f524" [[package]] name = "unicode-script" @@ -1728,9 +1792,9 @@ checksum = "b1d386ff53b415b7fe27b50bb44679e2cc4660272694b7b6f3326d8480823a94" [[package]] name = "url" -version = "2.5.0" +version = "2.5.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "31e6302e3bb753d46e83516cae55ae196fc0c309407cf11ab35cc51a4c2a4633" +checksum = "22784dbdf76fdde8af1aeda5622b546b422b6fc585325248a2bf9f5e41e94d6c" dependencies = [ "form_urlencoded", "idna", @@ -1800,9 +1864,9 @@ dependencies = [ [[package]] name = "utf8parse" -version = "0.2.1" +version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "711b9620af191e0cdc7468a8d14e709c3dcdb115b36f838e601583af800a370a" +checksum = "06abde3611657adf66d383f00b093d7faecc7fa57071cce2578660c9f1010821" [[package]] name = "vcpkg" @@ -1821,9 +1885,9 @@ dependencies = [ [[package]] name = "vergen" -version = "8.3.1" +version = "8.3.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e27d6bdd219887a9eadd19e1c34f32e47fa332301184935c6d9bca26f3cca525" +checksum = "2990d9ea5967266ea0ccf413a4aa5c42a93dbcfda9cb49a97de6931726b12566" dependencies = [ "anyhow", "cfg-if", @@ -1833,9 +1897,9 @@ dependencies = [ [[package]] name = "version_check" -version = "0.9.4" +version = "0.9.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "49874b5167b65d7193b8aba1567f5c7d93d001cafc34600cee003eda787e483f" +checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wasi" @@ -1845,34 +1909,35 @@ checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423" [[package]] name = "wasm-bindgen" -version = "0.2.92" +version = "0.2.93" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4be2531df63900aeb2bca0daaaddec08491ee64ceecbee5076636a3b026795a8" +checksum = "a82edfc16a6c469f5f44dc7b571814045d60404b55a0ee849f9bcfa2e63dd9b5" dependencies = [ "cfg-if", + "once_cell", "wasm-bindgen-macro", ] [[package]] name = "wasm-bindgen-backend" -version = "0.2.92" +version = "0.2.93" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "614d787b966d3989fa7bb98a654e369c762374fd3213d212cfc0251257e747da" +checksum = "9de396da306523044d3302746f1208fa71d7532227f15e347e2d93e4145dd77b" dependencies = [ "bumpalo", "log", "once_cell", "proc-macro2", "quote", - "syn 2.0.59", + "syn 2.0.75", "wasm-bindgen-shared", ] [[package]] name = "wasm-bindgen-macro" -version = "0.2.92" +version = "0.2.93" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a1f8823de937b71b9460c0c34e25f3da88250760bec0ebac694b49997550d726" +checksum = "585c4c91a46b072c92e908d99cb1dcdf95c5218eeb6f3bf1efa991ee7a68cccf" dependencies = [ "quote", "wasm-bindgen-macro-support", @@ -1880,22 +1945,22 @@ dependencies = [ [[package]] name = "wasm-bindgen-macro-support" -version = "0.2.92" +version = "0.2.93" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e94f17b526d0a461a191c78ea52bbce64071ed5c04c9ffe424dcb38f74171bb7" +checksum = "afc340c74d9005395cf9dd098506f7f44e38f2b4a21c6aaacf9a105ea5e1e836" dependencies = [ "proc-macro2", "quote", - "syn 2.0.59", + "syn 2.0.75", "wasm-bindgen-backend", "wasm-bindgen-shared", ] [[package]] name = "wasm-bindgen-shared" -version = "0.2.92" +version = "0.2.93" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "af190c94f2773fdb3729c55b007a722abb5384da03bc0986df4c289bf5567e96" +checksum = "c62a0a307cb4a311d3a07867860911ca130c3494e8c2719593806c08bc5d0484" [[package]] name = "weezl" @@ -1905,45 +1970,23 @@ checksum = "53a85b86a771b1c87058196170769dd264f66c0782acf1ae6cc51bfd64b39082" [[package]] name = "wide" -version = "0.7.15" +version = "0.7.28" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "89beec544f246e679fc25490e3f8e08003bc4bf612068f325120dad4cea02c1c" +checksum = "b828f995bf1e9622031f8009f8481a85406ce1f4d4588ff746d872043e855690" dependencies = [ "bytemuck", "safe_arch", ] -[[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-util" -version = "0.1.6" +version = "0.1.9" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f29e6f9198ba0d26b4c9f07dbe6f9ed633e1f3d5b8b414090084349e46a52596" +checksum = "cf221c93e13a30d793f7645a0e7762c55d169dbb0a49671918a2319d289b10bb" dependencies = [ - "winapi", + "windows-sys 0.59.0", ] -[[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" @@ -1962,11 +2005,20 @@ dependencies = [ "windows-targets", ] +[[package]] +name = "windows-sys" +version = "0.59.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e38bc4d79ed67fd075bcc251a1c39b32a1776bbe92e5bef1f0bf1f8c531853b" +dependencies = [ + "windows-targets", +] + [[package]] name = "windows-targets" -version = "0.52.5" +version = "0.52.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6f0713a46559409d202e70e28227288446bf7841d3211583a4b53e3f6d96e7eb" +checksum = "9b724f72796e036ab90c1021d4780d4d3d648aca59e491e6b98e725b84e99973" dependencies = [ "windows_aarch64_gnullvm", "windows_aarch64_msvc", @@ -1980,51 +2032,51 @@ dependencies = [ [[package]] name = "windows_aarch64_gnullvm" -version = "0.52.5" +version = "0.52.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7088eed71e8b8dda258ecc8bac5fb1153c5cffaf2578fc8ff5d61e23578d3263" +checksum = "32a4622180e7a0ec044bb555404c800bc9fd9ec262ec147edd5989ccd0c02cd3" [[package]] name = "windows_aarch64_msvc" -version = "0.52.5" +version = "0.52.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9985fd1504e250c615ca5f281c3f7a6da76213ebd5ccc9561496568a2752afb6" +checksum = "09ec2a7bb152e2252b53fa7803150007879548bc709c039df7627cabbd05d469" [[package]] name = "windows_i686_gnu" -version = "0.52.5" +version = "0.52.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "88ba073cf16d5372720ec942a8ccbf61626074c6d4dd2e745299726ce8b89670" +checksum = "8e9b5ad5ab802e97eb8e295ac6720e509ee4c243f69d781394014ebfe8bbfa0b" [[package]] name = "windows_i686_gnullvm" -version = "0.52.5" +version = "0.52.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "87f4261229030a858f36b459e748ae97545d6f1ec60e5e0d6a3d32e0dc232ee9" +checksum = "0eee52d38c090b3caa76c563b86c3a4bd71ef1a819287c19d586d7334ae8ed66" [[package]] name = "windows_i686_msvc" -version = "0.52.5" +version = "0.52.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "db3c2bf3d13d5b658be73463284eaf12830ac9a26a90c717b7f771dfe97487bf" +checksum = "240948bc05c5e7c6dabba28bf89d89ffce3e303022809e73deaefe4f6ec56c66" [[package]] name = "windows_x86_64_gnu" -version = "0.52.5" +version = "0.52.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4e4246f76bdeff09eb48875a0fd3e2af6aada79d409d33011886d3e1581517d9" +checksum = "147a5c80aabfbf0c7d901cb5895d1de30ef2907eb21fbbab29ca94c5b08b1a78" [[package]] name = "windows_x86_64_gnullvm" -version = "0.52.5" +version = "0.52.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "852298e482cd67c356ddd9570386e2862b5673c85bd5f88df9ab6802b334c596" +checksum = "24d5b23dc417412679681396f2b49f3de8c1473deb516bd34410872eff51ed0d" [[package]] name = "windows_x86_64_msvc" -version = "0.52.5" +version = "0.52.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "bec47e5bfd1bff0eeaf6d8b485cc1074891a197ab4225d504cb7a1ab88b02bf0" +checksum = "589f6da84c646204747d1270a2a5661ea66ed1cced2631d546fdfb155959f9ec" [[package]] name = "xmlparser" @@ -2037,3 +2089,24 @@ name = "xmlwriter" version = "0.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "ec7a2a501ed189703dba8b08142f057e887dfc4b2cc4db2d343ac6376ba3e0b9" + +[[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.75", +] diff --git a/Cargo.toml b/Cargo.toml index 21dbd3d..a627399 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,9 +1,12 @@ [package] name = "trgt" -version = "1.2.0" +version = "1.3.0" edition = "2021" build = "build.rs" +[workspace] +members = ["crates/pipeplot"] + [build-dependencies] vergen = { version = "8.2", features = ["git", "gitcl"] } @@ -20,11 +23,8 @@ kodama = "0.3" rand = "0.8" once_cell = "1.18" flate2 = "1.0" -resvg = "0.36" -usvg = "0.36" -tiny-skia = "0.11" -svg2pdf = "0.9" -tempfile = "3" rayon = "1.10" crossbeam-channel = "0.5" -semver = "1.0" \ No newline at end of file +pipeplot = { path = "crates/pipeplot" } +semver = "1.0" +tempfile = "3" diff --git a/README.md b/README.md index a54ff71..e1bc25b 100644 --- a/README.md +++ b/README.md @@ -136,8 +136,12 @@ tandem repeats at genome scale. 2024](https://www.nature.com/articles/s41587-023 - Merging now skips and logs problematic loci by default. Use the `--quit-on-errors` flag to terminate on errors. Statistics are logged post-merge, including counts of failed and skipped TRs. - `trgt validate` - Always outputs statistics directly to stdout and stderr instead of logging them. - - Bug fix: + - Bug fix: - Resolved issue with handling bgzip-compressed BED files. +- 1.3.0 + - Plotting code has been refactored as we prepare to revamp repeat visualizations + - The maximum number of reads per allele to plot can now be specified by `--max-allele-reads` + - bugfix: repeat identifiers are now permitted to contain commas ### DISCLAIMER diff --git a/crates/pipeplot/Cargo.toml b/crates/pipeplot/Cargo.toml new file mode 100644 index 0000000..aafb820 --- /dev/null +++ b/crates/pipeplot/Cargo.toml @@ -0,0 +1,13 @@ +[package] +name = "pipeplot" +version = "0.1.0" +edition = "2021" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] +resvg = "0.36" +usvg = "0.36" +tiny-skia = "0.11" +svg2pdf = "0.9" +tempfile = "3" diff --git a/crates/pipeplot/src/image.rs b/crates/pipeplot/src/image.rs new file mode 100644 index 0000000..eb467c0 --- /dev/null +++ b/crates/pipeplot/src/image.rs @@ -0,0 +1,46 @@ +use crate::{pdf, png, svg, PipePlot}; +use std::path::{Path, PathBuf}; +use tempfile::NamedTempFile; + +pub fn generate(plot: &PipePlot, path: &PathBuf) -> Result<(), String> { + let extension = Path::new(path) + .extension() + .and_then(|ext| ext.to_str()) + .ok_or_else(|| format!("Failed to get extension from path: {path:?}"))?; + let file_type = FileType::from_extension(extension) + .ok_or_else(|| format!("Unsupported file extension: {extension:?}"))?; + + let temp_dir = PathBuf::from(&path); + let temp_dir = temp_dir + .parent() + .ok_or_else(|| format!("Invalid path: {path:?}"))?; + let svg_temp_path = NamedTempFile::new_in(temp_dir) + .map_err(|e| format!("Failed to create temporary file: {e}"))? + .into_temp_path(); + + svg::generate(plot, &svg_temp_path.to_path_buf()); + match file_type { + FileType::Svg => svg_temp_path.persist(path).map_err(|e| e.to_string())?, + FileType::Png => png::render(&svg_temp_path.to_path_buf(), path)?, + FileType::Pdf => pdf::render(&svg_temp_path.to_path_buf(), path)?, + } + Ok(()) +} + +#[derive(Debug, PartialEq)] +enum FileType { + Svg, + Png, + Pdf, +} + +impl FileType { + fn from_extension(extension: &str) -> Option { + match extension.to_lowercase().as_str() { + "svg" => Some(FileType::Svg), + "png" => Some(FileType::Png), + "pdf" => Some(FileType::Pdf), + _ => None, + } + } +} diff --git a/crates/pipeplot/src/lib.rs b/crates/pipeplot/src/lib.rs new file mode 100644 index 0000000..3bb422e --- /dev/null +++ b/crates/pipeplot/src/lib.rs @@ -0,0 +1,18 @@ +/*! +This crate provides functionality to generate the so-called "pipe plots" +consisting of stacked horizontal pipes (bars). Each pipe consists of segments of +specified width, shape, and color. Pipe plots can be annotated with scales, +labels, and legends. The crate supports rendering of pipe plots as SVG, PNG, and +PDF images. + +Pipe plots are useful for representing pileups of sequenced reads. +*/ + +mod image; +mod pdf; +mod pipeplot; +mod png; +mod svg; + +pub use image::generate as generate_image; +pub use pipeplot::{Band, Color, Legend, Pipe, PipePlot, Seg, Shape}; diff --git a/src/trvz/pdf.rs b/crates/pipeplot/src/pdf.rs similarity index 58% rename from src/trvz/pdf.rs rename to crates/pipeplot/src/pdf.rs index a54185e..b65aee9 100644 --- a/src/trvz/pdf.rs +++ b/crates/pipeplot/src/pdf.rs @@ -1,12 +1,14 @@ +use std::path::PathBuf; use usvg::{TreeParsing, TreeTextToPath}; -pub fn render(svg_path: &str, pdf_path: &str) { - let svg = std::fs::read_to_string(svg_path).unwrap(); +pub fn render(svg_path: &PathBuf, pdf_path: &PathBuf) -> Result<(), String> { + let svg = std::fs::read_to_string(svg_path).map_err(|e| e.to_string())?; let options = usvg::Options::default(); let mut tree = usvg::Tree::from_str(&svg, &options).expect("Error parsing SVG"); let mut db = usvg::fontdb::Database::new(); db.load_system_fonts(); tree.convert_text(&db); let pdf = svg2pdf::convert_tree(&tree, svg2pdf::Options::default()); - std::fs::write(pdf_path, pdf).unwrap(); + std::fs::write(pdf_path, pdf).map_err(|e| e.to_string())?; + Ok(()) } diff --git a/crates/pipeplot/src/pipeplot.rs b/crates/pipeplot/src/pipeplot.rs new file mode 100644 index 0000000..8bb483e --- /dev/null +++ b/crates/pipeplot/src/pipeplot.rs @@ -0,0 +1,48 @@ +#[derive(Debug, PartialEq)] +pub enum Shape { + Rect, + HLine, + VLine, + None, + Tick(Option), +} + +pub type Color = String; + +#[derive(Debug, PartialEq)] +pub struct Seg { + pub width: u32, + pub color: Color, + pub shape: Shape, +} + +#[derive(Debug)] +pub struct Band { + pub pos: u32, // Position relative to pipe's start + pub width: u32, + pub color: Color, +} + +#[derive(Debug)] +pub struct Pipe { + pub xpos: u32, + pub ypos: u32, + pub height: u32, + pub segs: Vec, + pub bands: Vec, + pub outline: bool, +} + +#[derive(Debug)] +pub struct Legend { + pub xpos: u32, + pub ypos: u32, + pub height: u32, + pub labels: Vec<(String, String)>, +} + +#[derive(Debug)] +pub struct PipePlot { + pub pipes: Vec, + pub legend: Legend, +} diff --git a/src/trvz/png.rs b/crates/pipeplot/src/png.rs similarity index 51% rename from src/trvz/png.rs rename to crates/pipeplot/src/png.rs index b3ad496..f709219 100644 --- a/src/trvz/png.rs +++ b/crates/pipeplot/src/png.rs @@ -1,15 +1,18 @@ +use std::path::PathBuf; use usvg::{TreeParsing, TreeTextToPath}; -pub fn render(svg_path: &str, png_path: &str) { - let svg = std::fs::read(svg_path).unwrap(); +pub fn render(svg_path: &PathBuf, png_path: &PathBuf) -> Result<(), String> { + let svg = std::fs::read(svg_path).map_err(|e| e.to_string())?; let options = usvg::Options::default(); - let mut tree = usvg::Tree::from_data(&svg, &options).unwrap(); + let mut tree = usvg::Tree::from_data(&svg, &options).map_err(|e| e.to_string())?; let mut db = usvg::fontdb::Database::new(); db.load_system_fonts(); tree.convert_text(&db); let pixmap_size = tree.size.to_int_size(); - let mut pixmap = tiny_skia::Pixmap::new(pixmap_size.width(), pixmap_size.height()).unwrap(); + let mut pixmap = tiny_skia::Pixmap::new(pixmap_size.width(), pixmap_size.height()) + .ok_or("Unable to init image".to_string())?; let tree = resvg::Tree::from_usvg(&tree); tree.render(resvg::usvg::Transform::identity(), &mut pixmap.as_mut()); - pixmap.save_png(png_path).unwrap(); + pixmap.save_png(png_path).map_err(|e| e.to_string())?; + Ok(()) } diff --git a/crates/pipeplot/src/svg.rs b/crates/pipeplot/src/svg.rs new file mode 100644 index 0000000..e98ae98 --- /dev/null +++ b/crates/pipeplot/src/svg.rs @@ -0,0 +1,243 @@ +use crate::pipeplot::Color; +use crate::pipeplot::Legend; +use crate::pipeplot::Pipe; +use crate::pipeplot::PipePlot; +use crate::pipeplot::Shape; +use std::path::PathBuf; +use std::{fs::File, io::Write}; + +const DEFAULT_X_SCALE: f64 = 750.0; +const DEFAULT_Y_SCALE: f64 = 3.0; +const DEFAULT_PADDING: f64 = 12.0; + +pub fn generate(genotype_plot: &PipePlot, path: &PathBuf) { + let longest_pipe = get_longest_pipe(genotype_plot); + let x_scale = DEFAULT_X_SCALE / longest_pipe as f64; + let scale = (x_scale, DEFAULT_Y_SCALE); + let file = File::create(path).unwrap(); + let mut generator = Generator::new(scale, file, DEFAULT_PADDING); + generator.generate(genotype_plot); +} + +fn get_longest_pipe(plot: &PipePlot) -> u32 { + plot.pipes + .iter() + .map(|pipe| pipe.segs.iter().map(|seg| seg.width).sum()) + .max() + .unwrap_or(0) +} + +struct Generator { + scale: (f64, f64), + pad: f64, + file: File, +} + +impl Generator { + fn new(scale: (f64, f64), file: File, pad: f64) -> Self { + Self { scale, file, pad } + } + + pub fn generate(&mut self, pipe_plot: &PipePlot) { + let (width, height) = self.get_dimensions(pipe_plot); + self.start_svg(width, height); + self.add_background(); + + for pipe in &pipe_plot.pipes { + self.plot_pipe(pipe); + if pipe.outline { + self.plot_outline(pipe); + } + } + + self.plot_legend(&pipe_plot.legend); + self.end_svg(); + } + + fn plot_legend(&mut self, legend: &Legend) { + let base_x = self.to_x(legend.xpos) + self.pad; + let base_y = self.to_y(legend.ypos) + self.pad; + let height = self.to_y(legend.height); + let mut x = base_x; + for (label, color) in &legend.labels { + self.add_rect((x, base_y), (height, height), color, false); + x += height + 2.0; + let point = format!("x=\"{}\" y=\"{}\"", x, base_y + height - 1.0); + let height = r#"font-size="14px""#; + let style = r#"font-family="monospace" font-weight="bold""#; + let line = format!("{}", point, style, height, label); + writeln!(self.file, "{}", line).unwrap(); + x += 5.0 * (2 * label.len() as u32 + 1) as f64; + } + } + + fn plot_pipe(&mut self, pipe: &Pipe) { + let x = self.to_x(pipe.xpos) + self.pad; + let y = self.to_y(pipe.ypos) + self.pad; + let pipe_height = self.to_y(pipe.height); + + // Plot the main pipe + let mut x_cur = x; + + for seg in &pipe.segs { + let dims = (self.to_x(seg.width), pipe_height); + match seg.shape { + Shape::Rect => self.add_rect((x_cur, y), dims, &seg.color, true), + Shape::HLine => self.add_hline((x_cur, y), dims, &seg.color, 1.5), + Shape::Tick(label) => self.add_tick((x_cur, y), dims, &seg.color, label), + Shape::None | Shape::VLine => {} + } + + x_cur += self.to_x(seg.width); + } + + // Plot insertions + let mut x_cur = x; + for seg in &pipe.segs { + let dims = (self.to_x(seg.width), pipe_height); + + if seg.shape == Shape::VLine { + self.add_vline((x_cur, y), dims, &seg.color); + } + + x_cur += self.to_x(seg.width); + } + + // Plot bands + for band in &pipe.bands { + let beta_x = x + self.to_x(band.pos); + let dims = (self.to_x(1), pipe_height); + self.add_rect((beta_x, y), dims, &band.color, false); + } + } + + fn plot_outline(&mut self, pipe: &Pipe) { + let height = self.to_y(pipe.height); + let width = self.to_x(pipe.segs.iter().map(|seg| seg.width).sum()); + + let x = self.to_x(pipe.xpos) + self.pad; + let y = self.to_y(pipe.ypos) + self.pad; + + let dimensions = format!("width=\"{width}\" height=\"{height}\""); + let pos = format!("x=\"{x}\" y=\"{y}\""); + let style = r##"stroke="#000000" stroke-width="1.5" fill="transparent""##; + let line = format!("", dimensions, pos, style); + writeln!(self.file, "{}", line).unwrap(); + } + + fn add_rect(&mut self, pos: (f64, f64), dims: (f64, f64), color: &Color, add_highlight: bool) { + let (x, y) = pos; + let (w, h) = dims; + + let pos = format!("x=\"{}\" y=\"{}\"", x, y); + let dim = format!("height=\"{}\" width=\"{}\"", h, w); + + let style = format!("fill=\"{}\" stroke=\"{}\" stroke-width=\"0\"", color, color); + + let rect = format!("", pos, dim, style); + writeln!(self.file, "{}", rect).unwrap(); + + if add_highlight { + let pos = format!("x=\"{}\" y=\"{}\"", x, y + h * 0.18); + let dim = format!("height=\"{}\" width=\"{}\"", h / 3.0, w); + let highlight = format!("", pos, dim); + writeln!(self.file, "{}", highlight).unwrap(); + } + } + + fn add_hline(&mut self, pos: (f64, f64), dims: (f64, f64), color: &Color, stroke: f64) { + let x1 = pos.0; + let x2 = pos.0 + dims.0; + let y1 = pos.1 + dims.1 / 2.0; + let y2 = y1; + + let x1y1 = format!("x1=\"{}\" y1=\"{}\"", x1, y1); + let x2y2 = format!("x2=\"{}\" y2=\"{}\"", x2, y2); + + let style = format!("stroke=\"{}\" stroke-width=\"{}\"", color, stroke); + + let line = format!("", x1y1, x2y2, style); + writeln!(self.file, "{}", line).unwrap(); + } + + fn add_vline(&mut self, pos: (f64, f64), dims: (f64, f64), color: &Color) { + let x1 = pos.0; + let x2 = pos.0; + let y1 = pos.1; + let y2 = pos.1 + dims.1; + + let x1y1 = format!("x1=\"{}\" y1=\"{}\"", x1, y1); + let x2y2 = format!("x2=\"{}\" y2=\"{}\"", x2, y2); + + let stroke_width = 2.0_f64.min(self.to_x(1)); + let style = format!("stroke=\"{}\" stroke-width=\"{}\"", color, stroke_width); + + let line = format!("", x1y1, x2y2, style); + writeln!(self.file, "{}", line).unwrap(); + } + + fn add_tick(&mut self, pos: (f64, f64), dims: (f64, f64), color: &Color, label: Option) { + let x1 = pos.0; + let x2 = pos.0; + let y1 = pos.1; + let y2 = pos.1 + dims.1; + + let x1y1 = format!("x1=\"{}\" y1=\"{}\"", x1, y1); + let x2y2 = format!("x2=\"{}\" y2=\"{}\"", x2, y2); + + let stroke_width = 1.5; + let style = format!("stroke=\"{}\" stroke-width=\"{}\"", color, stroke_width); + + let line = format!("", x1y1, x2y2, style); + writeln!(self.file, "{}", line).unwrap(); + + if let Some(label) = label { + let point = format!("x=\"{}\" y=\"{}\"", x1, y1 - 2.0); // 2.0 is padding + let height = r#"font-size="14px""#; + let style = r#"font-family="monospace" font-weight="bold" text-anchor="middle""#; + let line = format!("{}", point, style, height, label); + writeln!(self.file, "{}", line).unwrap(); + } + } + + fn get_dimensions(&self, pipe_plot: &PipePlot) -> (f64, f64) { + let width = pipe_plot + .pipes + .iter() + .map(|p| p.xpos + p.segs.iter().map(|s| s.width).sum::()) + .max() + .unwrap_or(0); + let height = pipe_plot.legend.ypos + pipe_plot.legend.height; + + let xdim = self.to_x(width) + 2.0 * self.pad; + let ydim = self.to_y(height) + 2.0 * self.pad; + (xdim, ydim) + } + + fn start_svg(&mut self, width: f64, height: f64) { + writeln!(self.file, r#""#).unwrap(); + let line = r#"", width, height).unwrap(); + } + + fn end_svg(&mut self) { + writeln!(self.file, "").unwrap(); + } + + fn add_background(&mut self) { + writeln!( + self.file, + r#""# + ) + .unwrap(); + } + + fn to_x(&self, x: u32) -> f64 { + x as f64 * self.scale.0 + } + + fn to_y(&self, y: u32) -> f64 { + y as f64 * self.scale.1 + } +} diff --git a/docs/cli.md b/docs/cli.md index 00d978e..7c10428 100644 --- a/docs/cli.md +++ b/docs/cli.md @@ -59,6 +59,7 @@ Plotting: Advanced: - `--flank-len ` Length of flanking regions [default: 50]. +- `--max-allele-reads ` Max number of reads per allele to plot ## Merge command-line @@ -84,4 +85,4 @@ Options: - `-b, --repeats ` BED file with repeat coordinates and structure. Advanced: -- `--flank-len ` Length of flanking regions [default: 50]. \ No newline at end of file +- `--flank-len ` Length of flanking regions [default: 50]. diff --git a/src/cli.rs b/src/cli.rs index 1c44344..dafefa4 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -79,7 +79,7 @@ pub struct MergeArgs { #[clap(value_name = "FILE")] #[clap(help = "Write output to a file [standard output]")] #[arg(value_parser = check_prefix_path)] - pub output: Option, + pub output: Option, #[clap(help_heading("Advanced"))] #[clap(short = 'O')] @@ -184,7 +184,7 @@ pub struct GenotypeArgs { #[clap(help = "Prefix for output files")] #[clap(value_name = "OUTPUT_PREFIX")] #[arg(value_parser = check_prefix_path)] - pub output_prefix: String, + pub output_prefix: PathBuf, #[clap(long = "karyotype")] #[clap(short = 'k')] @@ -324,7 +324,7 @@ pub struct PlotArgs { #[clap(help = "Output image path")] #[clap(value_name = "IMAGE")] #[arg(value_parser = check_image_path)] - pub output_path: String, + pub output_path: PathBuf, #[clap(help_heading("Plotting"))] #[clap(long = "plot-type")] @@ -348,6 +348,12 @@ pub struct PlotArgs { #[clap(help = "Length of flanking regions")] #[clap(default_value = "50")] pub flank_len: usize, + + #[clap(help_heading("Advanced"))] + #[clap(long = "max-allele-reads")] + #[clap(value_name = "MAX_READS")] + #[clap(help = "Max number of reads per allele to plot")] + pub max_allele_reads: Option, } #[derive(Parser, Debug)] @@ -410,17 +416,17 @@ pub fn init_verbose(args: &Cli) { .init(); } -fn check_prefix_path(s: &str) -> Result { +fn check_prefix_path(s: &str) -> Result { let path = Path::new(s); if let Some(parent_dir) = path.parent() { if !parent_dir.as_os_str().is_empty() && !parent_dir.exists() { return Err(format!("Path does not exist: {}", parent_dir.display())); } } - Ok(s.to_string()) + Ok(PathBuf::from(s)) } -fn check_image_path(s: &str) -> Result { +fn check_image_path(s: &str) -> Result { let prefix_check = check_prefix_path(s)?; let path = Path::new(s); match path.extension().and_then(|ext| ext.to_str()) { diff --git a/src/commands/plot.rs b/src/commands/plot.rs index 028250c..8760adc 100644 --- a/src/commands/plot.rs +++ b/src/commands/plot.rs @@ -1,75 +1,24 @@ use crate::cli::PlotArgs; -use crate::trvz::{ - align::{align_reads, align_to_flanks}, - genotype_plot::plot_genotype, - input, pdf, - pipe_plot::{DisplayParams, PipePlot}, - png, svg, - waterfall_plot::plot_waterfall, -}; +use crate::trvz::color::pick_colors; +use crate::trvz::waterfall_plot::plot_waterfall; +use crate::trvz::{allele_plot::plot_alleles, input}; use crate::utils::{open_catalog_reader, open_genome_reader, Result}; -use std::path::{Path, PathBuf}; -use tempfile::NamedTempFile; - -#[derive(Debug, PartialEq)] -enum FileType { - Svg, - Png, - Pdf, -} - -impl FileType { - fn from_extension(extension: &str) -> Option { - match extension.to_lowercase().as_str() { - "svg" => Some(FileType::Svg), - "png" => Some(FileType::Png), - "pdf" => Some(FileType::Pdf), - _ => None, - } - } -} - -fn create_image(plot: PipePlot, path: String) -> Result<()> { - let extension = Path::new(&path) - .extension() - .and_then(|ext| ext.to_str()) - .ok_or("Failed to get file extension")?; - let file_type = FileType::from_extension(extension).ok_or("Unsupported file extension")?; - - let temp_dir = PathBuf::from(&path); - let temp_dir = temp_dir.parent().unwrap(); - let svg_temp_path = NamedTempFile::new_in(temp_dir).unwrap().into_temp_path(); - - svg::generate(&plot, svg_temp_path.to_str().unwrap()); - match file_type { - FileType::Svg => svg_temp_path.persist(path).unwrap(), - FileType::Png => png::render(svg_temp_path.to_str().unwrap(), &path), - FileType::Pdf => pdf::render(svg_temp_path.to_str().unwrap(), &path), - } - Ok(()) -} +use pipeplot::generate_image; pub fn trvz(args: PlotArgs) -> Result<()> { - let display_params = DisplayParams { - what_to_show: args.what_to_show, - max_width: 750, - }; - let catalog_reader = open_catalog_reader(&args.repeats_path)?; let genome_reader = open_genome_reader(&args.genome_path)?; - let locus = input::get_locus(catalog_reader, genome_reader, &args.tr_id, args.flank_len)?; - let reads = input::get_reads(&args.reads_path, &locus)?; + let reads = input::get_reads(&args.reads_path, &locus, args.max_allele_reads)?; + let colors = pick_colors(&locus.motifs); let pipe_plot = if args.plot_type == "allele" { - let genotype = input::get_genotype(&args.bcf_path, &locus)?; - let aligns = align_reads(&genotype, reads); - plot_genotype(&locus, &display_params, &genotype, &aligns) + let allele_seqs = input::get_alleles(&args.bcf_path, &locus)?; + plot_alleles(&locus, &args.what_to_show, &allele_seqs, &reads, colors) } else { - let aligns = align_to_flanks(&locus, reads); - plot_waterfall(&locus, &display_params, &aligns) + plot_waterfall(&locus, &args.what_to_show, &reads, &colors) }; - create_image(pipe_plot, args.output_path)?; + generate_image(&pipe_plot, &args.output_path)?; Ok(()) } diff --git a/src/merge/vcf_writer.rs b/src/merge/vcf_writer.rs index 1594765..9958c12 100644 --- a/src/merge/vcf_writer.rs +++ b/src/merge/vcf_writer.rs @@ -1,3 +1,5 @@ +use std::path::PathBuf; + use crate::utils::Result; use rust_htslib::bcf; @@ -22,11 +24,11 @@ impl VcfWriter { pub fn new( header: &bcf::Header, output_type: &Option, - output: Option<&String>, + output: Option<&PathBuf>, ) -> Result { let output_type = match (output_type, output) { (Some(output_type), _) => output_type.clone(), - (None, Some(path)) => Self::infer_output_type_from_extension(path)?, + (None, Some(path)) => Self::infer_output_type_from_extension(path.to_str().unwrap())?, (None, None) => OutputType::Vcf { is_uncompressed: true, level: None, diff --git a/src/trgt/genotype/genotype_flank.rs b/src/trgt/genotype/genotype_flank.rs index d1d2eca..3cc99be 100644 --- a/src/trgt/genotype/genotype_flank.rs +++ b/src/trgt/genotype/genotype_flank.rs @@ -106,6 +106,11 @@ fn get_trs_with_clustering<'a>( .unwrap() .0; + // No SNPs were called since both haplotypes are identical + if top_gt.0 == top_gt.1 { + return None; + } + let mut allele_assignment = Vec::new(); let mut assignment_tie_breaker = 1; let mut trs_by_allele = [Vec::new(), Vec::new()]; diff --git a/src/trgt/writers/write_bam.rs b/src/trgt/writers/write_bam.rs index e7eea7d..b253d97 100644 --- a/src/trgt/writers/write_bam.rs +++ b/src/trgt/writers/write_bam.rs @@ -113,6 +113,8 @@ impl BamWriter { let tr_tag = Aux::String(&locus.id); rec.push_aux(b"TR", tr_tag).unwrap(); + rec.push_aux(b"rq", Aux::Double(read.read_qual.unwrap_or(-1.0))) + .unwrap(); if let Some(meth) = &read.meth { let mc_tag = Aux::ArrayU8(meth.into()); rec.push_aux(b"MC", mc_tag).unwrap(); diff --git a/src/trvz/align.rs b/src/trvz/align.rs index a6fdc88..123aec5 100644 --- a/src/trvz/align.rs +++ b/src/trvz/align.rs @@ -1,88 +1,39 @@ -use crate::trvz::{ - locus::{Allele, Locus}, - read::Read, -}; -use bio::alignment::{pairwise::Aligner, Alignment}; +use super::read::Betas; #[derive(Debug, Clone)] -pub struct AlignInfo { - pub seq: String, - pub meth: Option>, - pub align: Alignment, - pub allele_index: usize, +pub struct AlleleAlign { + pub seq: (Align, Vec), + pub reads: Vec<(Align, Betas)>, } -pub fn align_reads(genotype: &Vec, reads: Vec) -> Vec { - let likely_max_len = get_longest_allele(genotype) + 50; - let mut aligner = get_aligner(likely_max_len); - let mut align_infos = Vec::new(); - - for read in reads { - let i = read.allele as usize; - if i < genotype.len() { - align_infos.push(AlignInfo { - align: aligner.global(read.seq.as_bytes(), genotype[i].seq.as_bytes()), - seq: read.seq, - meth: read.meth, - allele_index: i, - }); - } - } - - align_infos -} - -pub fn align_to_flanks(locus: &Locus, reads: Vec) -> Vec { - let max_read_len = reads.iter().map(|r| r.seq.len()).max().unwrap(); - let flank_len = locus.left_flank.len(); - - let mut aligns = Vec::new(); - let mut aligner = get_aligner(max_read_len); - for read in reads { - let mut ref_seq = locus.left_flank[locus.left_flank.len() - flank_len..].to_string(); - ref_seq += &read.seq[read.left_flank..read.seq.len() - read.right_flank]; - ref_seq += &locus.right_flank[..flank_len]; - - let align = aligner.global(read.seq.as_bytes(), ref_seq.as_bytes()); - - aligns.push(AlignInfo { - seq: read.seq, - meth: read.meth, - align, - allele_index: 0, - }); - } - aligns +#[derive(Debug, Clone)] +pub struct MotifBound { + pub start: usize, + pub end: usize, + pub motif_index: usize, } -type ScoreFunc = fn(u8, u8) -> i32; +pub type Align = Vec; -fn score(a: u8, b: u8) -> i32 { - if a == b { - 1i32 - } else { - -1i32 - } +#[derive(Debug, Clone)] +pub struct AlignSeg { + pub width: usize, + pub op: AlignOp, + pub seg_type: SegType, } -fn get_aligner<'a>(likely_max_len: usize) -> Aligner<&'a ScoreFunc> { - Aligner::with_capacity( - likely_max_len, - likely_max_len, - -5, - -1, - &(score as ScoreFunc), - ) +#[derive(Debug, Clone, PartialEq)] +pub enum AlignOp { + Match, + Subst, + Ins, + Del, } -fn get_longest_allele(genotype: &Vec) -> usize { - let mut max_len = 0; - - for allele in genotype { - if allele.seq.len() > max_len { - max_len = allele.seq.len(); - } - } - - max_len +#[derive(Debug, Clone, Copy, PartialEq, Hash, Eq)] +pub enum SegType { + // The integer is the motif index or motifs.len() for unsegmented regions + Tr(usize), + LeftFlank, + RightFlank, } diff --git a/src/trvz/align_allele.rs b/src/trvz/align_allele.rs new file mode 100644 index 0000000..5d0420d --- /dev/null +++ b/src/trvz/align_allele.rs @@ -0,0 +1,14 @@ +use super::align::AlleleAlign; +use super::align_consensus::align_consensus; +use super::align_reads::align_reads; +use super::locus::Locus; +use super::read::Read; + +pub fn get_allele_align(locus: &Locus, consensus: &str, reads: &[&Read]) -> AlleleAlign { + let (consensus_align, motif_bounds) = align_consensus(locus, consensus); + let read_aligns = align_reads(consensus, &consensus_align, reads); + AlleleAlign { + seq: (consensus_align, motif_bounds), + reads: read_aligns, + } +} diff --git a/src/trvz/align_consensus.rs b/src/trvz/align_consensus.rs new file mode 100644 index 0000000..9573146 --- /dev/null +++ b/src/trvz/align_consensus.rs @@ -0,0 +1,143 @@ +use super::align::{AlignOp, AlignSeg, MotifBound}; +use super::locus::Locus; +use crate::hmm::{build_hmm, get_events, HmmEvent}; +use crate::trvz::align::SegType; +use itertools::Itertools; + +/// Aligns a given allele to a perfect repeat as specified by the locus +/// definition +pub fn align_consensus(locus: &Locus, consensus: &str) -> (Vec, Vec) { + let mut align = vec![AlignSeg { + width: locus.left_flank.len(), + op: AlignOp::Match, + seg_type: SegType::LeftFlank, + }]; + + let motifs = locus + .motifs + .iter() + .map(|m| m.as_bytes().to_vec()) + .collect_vec(); + let query = &consensus[locus.left_flank.len()..consensus.len() - locus.right_flank.len()]; + let (query_align, motif_bounds) = align_motifs(&motifs, query); + let motif_bounds = motif_bounds + .into_iter() + .map(|bound| MotifBound { + start: locus.left_flank.len() + bound.start, + end: locus.left_flank.len() + bound.end, + motif_index: bound.motif_index, + }) + .collect_vec(); + align.extend(query_align); + align.push(AlignSeg { + width: locus.right_flank.len(), + op: AlignOp::Match, + seg_type: SegType::RightFlank, + }); + (align, motif_bounds) +} + +/// Aligns the given sequence to a perfect repeat composed of the given set of +/// motifs +fn align_motifs(motifs: &[Vec], seq: &str) -> (Vec, Vec) { + if seq.is_empty() { + return (Vec::new(), Vec::new()); + } + + let hmm = build_hmm(motifs); + let states = hmm.label(seq); + let motif_spans = hmm.label_motifs(&states); + let mut motif_by_base = vec![motifs.len(); seq.len()]; + for span in motif_spans { + motif_by_base[span.start..span.end].fill(span.motif_index); + } + + let events = get_events(&hmm, motifs, &states, seq.as_bytes()); + let mut align = Vec::new(); + let mut bound_events = Vec::new(); + let mut base_pos = 0; + + for (event, group) in &events.iter().group_by(|e| *e) { + let seg_type = if base_pos < motif_by_base.len() { + SegType::Tr(motif_by_base[base_pos]) + } else { + assert_eq!(base_pos, motif_by_base.len()); + SegType::Tr(motif_by_base[base_pos.saturating_sub(1)]) + }; + + let width = group.count(); + match event { + HmmEvent::Trans => (), + HmmEvent::MotifStart | HmmEvent::MotifEnd => { + if let SegType::Tr(index) = seg_type { + bound_events.push((event, base_pos, index)); + } + } + HmmEvent::Del => align.push(AlignSeg { + width: 0, + op: AlignOp::Ins, + seg_type, + }), + HmmEvent::Ins => align.push(AlignSeg { + width, + op: AlignOp::Del, + seg_type, + }), + HmmEvent::Match => align.push(AlignSeg { + width, + op: AlignOp::Match, + seg_type, + }), + HmmEvent::Mismatch => align.push(AlignSeg { + width, + op: AlignOp::Subst, + seg_type, + }), + HmmEvent::Skip => { + assert_eq!(seg_type, SegType::Tr(motifs.len())); + align.push(AlignSeg { + width, + op: AlignOp::Match, + seg_type, + }) + } + } + + base_pos += match event { + HmmEvent::Match | HmmEvent::Mismatch | HmmEvent::Ins | HmmEvent::Skip => width, + HmmEvent::MotifStart | HmmEvent::MotifEnd | HmmEvent::Del | HmmEvent::Trans => 0, + }; + } + + assert_eq!(base_pos, seq.len()); + + let mut merged_align = Vec::new(); + let iter = align + .into_iter() + .group_by(|seg| (seg.op.clone(), seg.seg_type)); + for ((op, segment), group) in &iter { + let width: usize = group.into_iter().map(|s| s.width).sum(); + merged_align.push(AlignSeg { + width, + op, + seg_type: segment, + }); + } + + let motif_bounds = bound_events + .chunks(2) + .map(|rec| { + let (start, end) = rec + .iter() + .collect_tuple() + .expect("Unexpected motif segmentation"); + MotifBound { + start: start.1, + end: end.1, + motif_index: start.2, + } + }) + .collect(); + + (merged_align, motif_bounds) +} diff --git a/src/trvz/align_reads.rs b/src/trvz/align_reads.rs new file mode 100644 index 0000000..4ddba3d --- /dev/null +++ b/src/trvz/align_reads.rs @@ -0,0 +1,139 @@ +use super::align::Align; +use super::align::{AlignOp, AlignSeg}; +use super::read::{project_betas, Betas, Read}; +use bio::alignment::pairwise::Aligner; +use bio::alignment::Alignment as BioAlign; +use itertools::Itertools; + +type BioOp = bio::alignment::AlignmentOperation; + +/// Align reads to the consensus sequence +pub fn align_reads( + consensus: &str, + consensus_align: &[AlignSeg], + reads: &[&Read], +) -> Vec<(Align, Betas)> { + let likely_max_len = consensus.len() + 50; + let mut aligner = get_aligner(likely_max_len); + let mut read_aligns = Vec::new(); + + for read in reads { + let bio_align = aligner.global(read.seq.as_bytes(), consensus.as_bytes()); + let align = convert(consensus_align, &bio_align); + let betas = project_betas(&bio_align, &read.betas); + read_aligns.push((align, betas)); + } + + read_aligns +} + +/// Convert a rust-bio alignment into an internal alignment +fn convert(consensus_align: &[AlignSeg], bio_align: &BioAlign) -> Align { + assert!(bio_align.xstart == 0); + assert!(bio_align.ystart == 0); + + let mut seg_type_by_ref = Vec::new(); + for align_seg in consensus_align { + let seg_type = align_seg.seg_type; + seg_type_by_ref.extend(match align_seg.op { + AlignOp::Del | AlignOp::Match | AlignOp::Subst => vec![seg_type; align_seg.width], + AlignOp::Ins => vec![], + }); + } + + let mut ref_pos = 0; + let mut ops_and_segs = Vec::new(); + for op in &bio_align.operations { + let seg_type = if ref_pos == seg_type_by_ref.len() { + // Handle trailing insertion + assert_eq!(*op, BioOp::Ins); + seg_type_by_ref[ref_pos - 1] + } else { + seg_type_by_ref[ref_pos] + }; + ops_and_segs.push((*op, seg_type)); + ref_pos += match *op { + BioOp::Match => 1, + BioOp::Subst => 1, + BioOp::Del => 1, + BioOp::Ins => 0, + _ => panic!("Unhandled operation {op:?}"), + }; + } + + let mut align = Vec::new(); + let mut ref_pos = 0; + let mut seq_pos = 0; + + for ((bio_op, seg_type), group) in &ops_and_segs.iter().group_by(|rec| *rec) { + let seg_type = *seg_type; + let run_len = group.count(); + + let align_seg = match *bio_op { + BioOp::Match => AlignSeg { + width: run_len, + op: AlignOp::Match, + seg_type, + }, + BioOp::Subst => AlignSeg { + width: run_len, + op: AlignOp::Subst, + seg_type, + }, + BioOp::Del => AlignSeg { + width: run_len, + op: AlignOp::Del, + seg_type, + }, + BioOp::Ins => AlignSeg { + width: 0, + op: AlignOp::Ins, + seg_type, + }, + _ => panic!("No logic to handle {bio_op:?}"), + }; + + align.push(align_seg); + + ref_pos += match *bio_op { + BioOp::Match => run_len, + BioOp::Subst => run_len, + BioOp::Del => run_len, + BioOp::Ins => 0, + _ => panic!("Unhandled operation {:?}", *bio_op), + }; + + seq_pos += match *bio_op { + BioOp::Match => run_len, + BioOp::Subst => run_len, + BioOp::Del => 0, + BioOp::Ins => run_len, + _ => panic!("Unhandled operation {:?}", *bio_op), + }; + } + + assert_eq!(bio_align.x_aln_len(), seq_pos); + assert_eq!(bio_align.y_aln_len(), ref_pos); + + align +} + +type ScoreFunc = fn(u8, u8) -> i32; + +fn score(a: u8, b: u8) -> i32 { + if a == b { + 1i32 + } else { + -1i32 + } +} + +fn get_aligner<'a>(likely_max_len: usize) -> Aligner<&'a ScoreFunc> { + Aligner::with_capacity( + likely_max_len, + likely_max_len, + -5, + -1, + &(score as ScoreFunc), + ) +} diff --git a/src/trvz/allele_plot.rs b/src/trvz/allele_plot.rs new file mode 100644 index 0000000..10b0c36 --- /dev/null +++ b/src/trvz/allele_plot.rs @@ -0,0 +1,240 @@ +use super::align::{MotifBound, SegType}; +use super::align_allele::get_allele_align; +use super::color::{get_meth_colors, Color, ColorMap}; +use super::read::{Betas, Read}; +use crate::trvz::align::{Align, AlignOp}; +use crate::trvz::locus::Locus; +use itertools::Itertools; +use pipeplot::{Band, Legend, Pipe, PipePlot, Seg, Shape}; + +pub fn plot_alleles( + locus: &Locus, + what_to_show: &str, + allele_seqs: &[String], + reads: &[Read], + colors: ColorMap, +) -> PipePlot { + let aligns_by_allele = allele_seqs + .iter() + .enumerate() + .map(|(index, allele_seq)| { + let allele_reads = reads + .iter() + .filter(|r| r.allele == index as i32) + .collect_vec(); + get_allele_align(locus, allele_seq, &allele_reads) + }) + .collect_vec(); + let bounds_by_allele = aligns_by_allele + .iter() + .map(|align| &align.seq.1) + .collect_vec(); + let tick_spacing = get_tick_spacing(&bounds_by_allele); + + let height = 4; + let xpos = 0; + let mut ypos = 0; + let mut pipes = Vec::new(); + for (allele_index, allele) in aligns_by_allele.iter().enumerate() { + let pipe = get_scales( + xpos, + ypos, + height / 3, + tick_spacing, + locus.motifs.len(), + &allele.seq.1, + ); + pipes.push(pipe); + ypos += height / 3; + let pipe = get_pipe( + xpos, + ypos, + height, + &allele.seq.0, + &Vec::new(), + &colors, + true, + ); + pipes.push(pipe); + ypos += 5; + + // TODO: Confirm that this allele / index correspondence is always correct + for (align, betas) in &allele.reads { + let (colors, betas) = if what_to_show == "meth" { + (get_meth_colors(&locus.motifs), betas.clone()) + } else { + (colors.clone(), Vec::new()) + }; + + let pipe = get_pipe(xpos, ypos, height, align, &betas, &colors, false); + pipes.push(pipe); + ypos += 5; + } + + if allele_index + 1 != aligns_by_allele.len() { + ypos += 7; + } + } + + let mut labels = Vec::new(); + for (index, motif) in locus.motifs.iter().enumerate() { + let color = colors.get(&SegType::Tr(index)).unwrap().to_string(); + labels.push((motif.clone(), color)); + } + if what_to_show == "meth" { + labels.push(("Methylated".to_string(), Color::Grad(1.0).to_string())); + labels.push(("Unmethylated".to_string(), Color::Grad(0.0).to_string())); + } + + ypos += 1; + + let legend = Legend { + xpos, + ypos, + height, + labels, + }; + + PipePlot { pipes, legend } +} + +fn get_scales( + mut xpos: u32, + ypos: u32, + height: u32, + tick_spacing: usize, + motif_count: usize, + bounds: &[MotifBound], +) -> Pipe { + xpos += bounds.first().unwrap().start as u32; + let mut segs = Vec::new(); + + for (motif_index, group) in &bounds.iter().group_by(|bound| bound.motif_index) { + let mut tick_index = 0; + for bound in group { + if tick_index % tick_spacing == 0 { + let tick_label = if motif_index < motif_count { + Some(tick_index as u32) + } else { + None + }; + segs.push(Seg { + width: 0, + color: Color::Black.to_string(), + shape: Shape::Tick(tick_label), + }); + } + segs.push(Seg { + width: (bound.end - bound.start) as u32, + color: Color::Black.to_string(), + shape: Shape::None, + }); + tick_index += 1; + } + if tick_index % tick_spacing == 0 { + let tick_label = if motif_index < motif_count { + Some(tick_index as u32) + } else { + None + }; + segs.push(Seg { + width: 0, + color: Color::Black.to_string(), + shape: Shape::Tick(tick_label), + }); + } + } + + let mut culled_segs = Vec::new(); + for (is_tick, group) in &segs + .into_iter() + .group_by(|seg| matches!(seg.shape, Shape::Tick(_))) + { + if is_tick { + culled_segs.push(group.last().unwrap()); + } else { + culled_segs.extend(group); + } + } + + Pipe { + xpos, + ypos, + height, + segs: culled_segs, + bands: Vec::new(), + outline: false, + } +} + +fn get_pipe( + xpos: u32, + ypos: u32, + height: u32, + align: &Align, + betas: &Betas, + colors: &ColorMap, + outline: bool, +) -> Pipe { + let mut segs = Vec::new(); + for align_seg in align { + let shape = match align_seg.op { + AlignOp::Del => Shape::HLine, + AlignOp::Ins => Shape::VLine, + AlignOp::Match | AlignOp::Subst => Shape::Rect, + }; + let color = if align_seg.op == AlignOp::Match { + colors.get(&align_seg.seg_type).unwrap() + } else if align_seg.op == AlignOp::Subst { + &Color::Gray + } else { + &Color::Black + }; + segs.push(Seg { + width: align_seg.width as u32, + color: color.to_string(), + shape, + }); + } + + let mut bands = Vec::new(); + + for beta in betas { + // Band width is 2 to cover the entire CpG + let color = Color::Grad(beta.value); + bands.push(Band { + pos: beta.pos as u32, + width: 2, + color: color.to_string(), + }); + } + + Pipe { + xpos, + ypos, + height, + segs, + bands, + outline, + } +} + +fn get_tick_spacing(bounds_by_allele: &[&Vec]) -> usize { + let max_motif_count = bounds_by_allele + .iter() + .map(|bounds| bounds.len()) + .max() + .unwrap_or(1); + + match max_motif_count { + 0..=20 => 1, + 21..=100 => 5, + 101..=200 => 10, + 201..=500 => 20, + 501..=1000 => 50, + 1001..=1500 => 100, + 1501..=2000 => 200, + 2001..=3000 => 250, + _ => 500, + } +} diff --git a/src/trvz/color.rs b/src/trvz/color.rs new file mode 100644 index 0000000..bb96326 --- /dev/null +++ b/src/trvz/color.rs @@ -0,0 +1,94 @@ +use super::align::SegType; +use std::{collections::HashMap, fmt}; + +pub type ColorMap = HashMap; + +pub fn pick_colors(motifs: &[String]) -> ColorMap { + let tr_colors = [ + Color::Blue, + Color::Purple, + Color::Orange, + Color::Pink, + Color::Yellow, + Color::Green, + Color::Red, + Color::Khaki, + Color::PaleRed, + Color::PaleBlue, + ]; + + let mut colors = HashMap::new(); + colors.insert(SegType::LeftFlank, Color::Teal); + colors.insert(SegType::RightFlank, Color::Teal); + + for index in 0..motifs.len() { + let color = tr_colors[index % tr_colors.len()].clone(); + colors.insert(SegType::Tr(index), color); + } + colors.insert(SegType::Tr(motifs.len()), Color::Gray); + + colors +} + +pub fn get_meth_colors(motifs: &[String]) -> ColorMap { + let mut colors = HashMap::new(); + colors.insert(SegType::LeftFlank, Color::Teal); + colors.insert(SegType::RightFlank, Color::Teal); + + for index in 0..motifs.len() + 1 { + colors.insert(SegType::Tr(index), Color::Gray); + } + + colors +} + +#[derive(Debug, PartialEq, Clone)] +pub enum Color { + Purple, + Blue, + Orange, + Teal, + Gray, + LightGray, + Black, + Green, + Pink, + Yellow, + Red, + Khaki, + PaleRed, + PaleBlue, + Grad(f64), +} + +impl fmt::Display for Color { + fn fmt(&self, formatter: &mut fmt::Formatter) -> fmt::Result { + match self { + Color::Purple => write!(formatter, "#814ED1"), + Color::Blue => write!(formatter, "#1383C6"), + Color::Orange => write!(formatter, "#E16A2C"), + Color::Teal => write!(formatter, "#009CA2"), + Color::Gray => write!(formatter, "#BABABA"), + Color::LightGray => write!(formatter, "#D1D1D1"), + Color::Black => write!(formatter, "#000000"), + Color::Pink => write!(formatter, "#ED3981"), + Color::Yellow => write!(formatter, "#EFCD17"), + Color::Green => write!(formatter, "#009D4E"), + Color::Red => write!(formatter, "#E3371E"), + Color::Khaki => write!(formatter, "#F0E68C"), + Color::PaleRed => write!(formatter, "#FF4858"), + Color::PaleBlue => write!(formatter, "#46B2E8"), + Color::Grad(value) => write!(formatter, "{}", get_gradient(*value)), + } + } +} + +fn get_gradient(value: f64) -> String { + let blue: (u8, u8, u8) = (0, 73, 255); + let red: (u8, u8, u8) = (255, 0, 0); + let mix_red = (blue.0 as f64 * (1.0 - value) + red.0 as f64 * value).round() as u8; + let mix_green = (blue.1 as f64 * (1.0 - value) + red.1 as f64 * value).round() as u8; + let mix_blue = (blue.2 as f64 * (1.0 - value) + red.2 as f64 * value).round() as u8; + + format!("#{:02X}{:02X}{:02X}", mix_red, mix_green, mix_blue) +} diff --git a/src/trvz/genotype_plot.rs b/src/trvz/genotype_plot.rs deleted file mode 100644 index 22225f5..0000000 --- a/src/trvz/genotype_plot.rs +++ /dev/null @@ -1,94 +0,0 @@ -use crate::trvz::{ - align::AlignInfo, - locus::{Allele, BaseLabel, Locus}, - pipe_plot::{ - get_allele_pipe, get_color, get_pipe, AlignOp, DisplayParams, Legend, PipePlot, PlotPanel, - }, - struc::RegionLabel, -}; - -pub fn plot_genotype( - locus: &Locus, - params: &DisplayParams, - genotype: &[Allele], - aligns: &[AlignInfo], -) -> PipePlot { - let tick_spacing = get_tick_spacing(locus, genotype); - - let mut allele_plots = Vec::new(); - for (allele_index, allele) in genotype.iter().enumerate() { - let allele_aligns: Vec<&AlignInfo> = aligns - .iter() - .filter(|info| info.allele_index == allele_index) - .collect(); - allele_plots.push(plot_allele( - locus, - allele, - &allele_aligns, - params.what_to_show == "meth", - tick_spacing, - )); - } - - let mut labels = Vec::new(); - for motif in &locus.motifs { - let color = get_color(locus, AlignOp::Match, &RegionLabel::Tr(0, 0, motif.clone())); - labels.push((motif.clone(), color)); - } - - let height = 4; - let legend = Legend { labels, height }; - PipePlot { - panels: allele_plots, - legend, - } -} - -fn plot_allele( - locus: &Locus, - allele: &Allele, - aligns: &Vec<&AlignInfo>, - show_betas: bool, - tick_spacing: usize, -) -> PlotPanel { - let labels = if show_betas { - &allele.flank_labels - } else { - &allele.region_labels - }; - - let mut allele_plot = Vec::new(); - - allele_plot.push(get_allele_pipe( - locus, - tick_spacing, - &allele.region_labels, - &allele.base_labels, - )); - for align in aligns { - allele_plot.push(get_pipe(locus, labels, align, show_betas)); - } - - allele_plot -} - -fn get_tick_spacing(_locus: &Locus, genotype: &[Allele]) -> usize { - let max_motif_count = genotype - .iter() - .map(|a| { - a.base_labels - .iter() - .filter(|l| **l == BaseLabel::MotifBound) - .count() - }) - .max() - .unwrap(); - - if max_motif_count <= 150 { - 1 - } else if max_motif_count <= 250 { - 5 - } else { - 10 - } -} diff --git a/src/trvz/input.rs b/src/trvz/input.rs index fcb2832..0372292 100644 --- a/src/trvz/input.rs +++ b/src/trvz/input.rs @@ -1,23 +1,15 @@ -use crate::hmm::{build_hmm, get_events, Hmm, HmmEvent}; -use crate::trvz::{ - locus::{self, Allele, BaseLabel, Locus}, - read::Read, - struc::RegionLabel, -}; +use crate::trvz::locus::{self, Locus}; use crate::utils::Result; use itertools::Itertools; -use rust_htslib::{ - bam::{self, record::Aux, Read as BamRead}, - bcf::{self, record::GenotypeAllele::UnphasedMissing, Read as BcfRead, Record}, - faidx, -}; -use std::{ - io::{BufRead, BufReader, Read as ioRead}, - path::PathBuf, - str, -}; +use rust_htslib::bam::{self, record::Aux, Read as BamRead}; +use rust_htslib::bcf::{self, record::GenotypeAllele::UnphasedMissing, Read as BcfRead, Record}; +use rust_htslib::faidx; +use std::collections::HashSet; +use std::io::{BufRead, BufReader, Read as ioRead}; +use std::path::PathBuf; +use std::str; -type RegionLabels = Vec; +use super::read::{Beta, Betas, Read}; #[derive(Debug)] pub struct Span { @@ -26,13 +18,16 @@ pub struct Span { pub end: usize, } -pub fn get_genotype(bcf_path: &PathBuf, locus: &Locus) -> Result> { +pub fn get_alleles(bcf_path: &PathBuf, locus: &Locus) -> Result> { let mut bcf = bcf::Reader::from_path(bcf_path).unwrap(); for record in bcf.records() { - let record = record.unwrap(); + let record = record.map_err(|e| e.to_string())?; let tr_id = record.info(b"TRID").string().unwrap().unwrap(); - let tr_id = str::from_utf8(tr_id.to_vec()[0]).unwrap(); + let tr_id = tr_id + .iter() + .map(|elem| str::from_utf8(elem).unwrap()) + .join(","); if tr_id != locus.id { continue; @@ -43,22 +38,8 @@ pub fn get_genotype(bcf_path: &PathBuf, locus: &Locus) -> Result> { return Err(format!("Missing genotype for TRID={}", tr_id)); } - let allele_seqs = get_allele_seqs(locus, &record); - let region_labels_by_allele = get_region_labels(locus, &allele_seqs, &record); - let flank_labels_by_allele = get_flank_labels(locus, ®ion_labels_by_allele); - let base_labels_by_allele = label_with_hmm(locus, &allele_seqs); - - let mut genotype = Vec::new(); - for (index, seq) in allele_seqs.into_iter().enumerate() { - genotype.push(Allele { - seq, - region_labels: region_labels_by_allele[index].clone(), - flank_labels: flank_labels_by_allele[index].clone(), - base_labels: base_labels_by_allele[index].clone(), - }); - } - - return Ok(genotype); + let alleles = get_allele_seqs(locus, &record); + return Ok(alleles); } Err(format!("TRID={} missing", &locus.id)) } @@ -79,7 +60,11 @@ pub fn get_locus( Err(format!("Unable to find locus {}", tr_id)) } -pub fn get_reads(bam_path: &PathBuf, locus: &Locus) -> Result> { +pub fn get_reads( + bam_path: &PathBuf, + locus: &Locus, + max_allele_reads: Option, +) -> Result> { let mut reads = bam::IndexedReader::from_path(bam_path).unwrap(); // This assumes that TRGT outputs flanks shorter than 1Kbps in length. We may want // to implement a more flexible mechanism for handling flank lengths here and elsewhere. @@ -108,22 +93,7 @@ pub fn get_reads(bam_path: &PathBuf, locus: &Locus) -> Result> { continue; } - let meth = match read.aux(b"MC") { - Ok(Aux::ArrayU8(value)) => { - if !value.is_empty() { - Some(value.iter().collect::>()) - } else { - None - } - } - Ok(_) => { - return Err(format!( - "malformed MC tag in read {:?}.", - String::from_utf8(read.qname().to_vec()).unwrap() - )) - } - Err(_) => None, - }; + let meth = parse_meth(&read, &seq)?; let allele = match read.aux(b"AL") { Ok(Aux::I32(value)) => value, @@ -135,7 +105,7 @@ pub fn get_reads(bam_path: &PathBuf, locus: &Locus) -> Result> { } Err(_) => { return Err(format!( - "malformatted read. Expected AL tag not found: {:?}", + "Malformed read. Expected AL tag not found: {:?}", String::from_utf8(read.qname().to_vec()).unwrap() )) } @@ -155,13 +125,13 @@ pub fn get_reads(bam_path: &PathBuf, locus: &Locus) -> Result> { } Ok(_) => { return Err(format!( - "malformatted FL tag in read {:?}.", + "Malformed FL tag in read {:?}.", String::from_utf8(read.qname().to_vec()).unwrap() )) } Err(_) => { return Err(format!( - "malformatted read. Expected FL tag not found: {:?}", + "Malformed read. Expected FL tag not found: {:?}", String::from_utf8(read.qname().to_vec()).unwrap() )) } @@ -172,9 +142,34 @@ pub fn get_reads(bam_path: &PathBuf, locus: &Locus) -> Result> { left_flank, right_flank, allele, - meth, + betas: meth, }); } + let seqs = if let Some(max_reads) = max_allele_reads { + let unique_alleles: Vec = seqs + .iter() + .map(|x| x.allele) + .collect::>() + .into_iter() + .collect(); + unique_alleles + .into_iter() + .map(|allele| { + let num_allele_reads = seqs.iter().filter(|x| x.allele == allele).count(); + + // GS: this only works if spanning.bam is sorted by TR length + let step = std::cmp::max(1_usize, num_allele_reads / max_reads); + seqs.clone() + .into_iter() + .filter(|seq| seq.allele == allele) + .step_by(step) + .take(max_reads) + .collect() + }) + .concat() + } else { + seqs + }; Ok(seqs) } @@ -198,130 +193,46 @@ fn get_allele_seqs(locus: &Locus, record: &Record) -> Vec { alleles } -fn get_region_labels(locus: &Locus, alleles: &[String], record: &Record) -> Vec { - let lf_len = locus.left_flank.len(); - let rf_len = locus.right_flank.len(); - - let mut labels_by_hap = Vec::new(); - let ms_field = record.format(b"MS").string().unwrap(); - let ms_field = str::from_utf8(ms_field.to_vec()[0]).unwrap(); - for (allele_index, spans) in ms_field.split(',').enumerate() { - let allele_len = alleles[allele_index].len(); - if spans == "." { - let tr_start = lf_len; - let tr_end = allele_len - rf_len; - - labels_by_hap.push(vec![ - RegionLabel::Flank(0, tr_start), - RegionLabel::Other(tr_start, tr_end), - RegionLabel::Flank(tr_end, allele_len), - ]); - continue; - } - let mut labels = vec![RegionLabel::Flank(0, locus.left_flank.len())]; - let mut last_seg_end = locus.left_flank.len(); - for span in spans.split('_') { - let (motif_index, start, end) = span - .trim_end_matches(')') - .split(&['(', '-']) - .map(|s| s.parse::().unwrap()) - .collect_tuple() - .unwrap(); - let motif = locus.motifs[motif_index].clone(); - let start = start + locus.left_flank.len(); - let end = end + locus.left_flank.len(); - - if start != last_seg_end { - labels.push(RegionLabel::Seq(last_seg_end, start)); +fn parse_meth(read: &bam::Record, seq: &str) -> Result { + let values = match read.aux(b"MC") { + Ok(Aux::ArrayU8(value)) => { + if !value.is_empty() { + value.iter().collect::>() + } else { + return Ok(Vec::new()); } - labels.push(RegionLabel::Tr(start, end, motif)); - last_seg_end = end; } - - if last_seg_end != allele_len - rf_len { - let seg_end = allele_len - rf_len; - labels.push(RegionLabel::Seq(last_seg_end, seg_end)); - last_seg_end = seg_end; + Ok(_) => { + let read_name = String::from_utf8(read.qname().to_vec()).unwrap(); + return Err(format!("Bad MC tag in read {read_name}")); + } + Err(_) => return Ok(Vec::new()), + }; + + let mut betas = Vec::new(); + let mut cpg_count = 0; + for pos in 0..seq.len() { + let is_cpg = pos + 1 < seq.len() && &seq[pos..pos + 2] == "CG"; + if is_cpg { + if cpg_count == values.len() { + let read_name = String::from_utf8(read.qname().to_vec()).unwrap(); + return Err(format!("Bad MC tag in read {read_name}")); + } + let level = values[cpg_count] as f64 / 255.0; + let beta = Beta { pos, value: level }; + betas.push(beta); + cpg_count += 1; } - - labels.push(RegionLabel::Flank( - last_seg_end, - last_seg_end + locus.right_flank.len(), - )); - labels_by_hap.push(labels); - } - - labels_by_hap -} - -fn get_flank_labels(locus: &Locus, all_labels_by_allele: &Vec) -> Vec { - let mut flank_labels_by_allele = Vec::new(); - for all_labels in all_labels_by_allele { - let tr_len = all_labels - .iter() - .map(|l| match l { - RegionLabel::Flank(_, _) => 0, - RegionLabel::Other(start, end) => end - start, - RegionLabel::Seq(start, end) => end - start, - RegionLabel::Tr(start, end, _) => end - start, - }) - .sum::(); - - let tr_start = locus.left_flank.len(); - let tr_end = tr_start + tr_len; - let allele_end = tr_end + locus.right_flank.len(); - - flank_labels_by_allele.push(vec![ - RegionLabel::Flank(0, tr_start), - RegionLabel::Other(tr_start, tr_end), - RegionLabel::Flank(tr_end, allele_end), - ]); } - flank_labels_by_allele -} -fn label_with_hmm(locus: &Locus, alleles: &[String]) -> Vec> { - let motifs = locus - .motifs - .iter() - .map(|m| m.as_bytes().to_vec()) - .collect_vec(); - let hmm = build_hmm(&motifs); + let all_cpgs_present = cpg_count == values.len(); + let last_cpg_split = + seq.bytes().last().unwrap_or(b'X') == b'C' && cpg_count + 1 == values.len(); - let mut labels_by_allele = Vec::new(); - - for allele in alleles { - let query = &allele[locus.left_flank.len()..allele.len() - locus.right_flank.len()]; - let mut labels = vec![BaseLabel::Match; locus.left_flank.len()]; - let states = hmm.label(query); - labels.extend(get_base_labels( - &locus.motifs, - &hmm, - &states, - query.as_bytes(), - )); - labels.extend(vec![BaseLabel::Match; locus.right_flank.len()]); - labels_by_allele.push(labels); + if !(all_cpgs_present || last_cpg_split) { + let read_name = String::from_utf8(read.qname().to_vec()).unwrap(); + return Err(format!("Bad MC tag in read {read_name}")); } - labels_by_allele -} - -fn get_base_labels(motifs: &[String], hmm: &Hmm, states: &[usize], query: &[u8]) -> Vec { - let motifs = motifs.iter().map(|m| m.as_bytes().to_vec()).collect_vec(); - let events = get_events(hmm, &motifs, states, query); - - let mut base_labels = Vec::new(); - for event in events { - match event { - HmmEvent::Match | HmmEvent::Skip => base_labels.push(BaseLabel::Match), - HmmEvent::Ins => base_labels.push(BaseLabel::NoMatch), - HmmEvent::Del => base_labels.push(BaseLabel::Skip), - HmmEvent::Mismatch => base_labels.push(BaseLabel::Mismatch), - HmmEvent::MotifStart => base_labels.push(BaseLabel::MotifBound), - HmmEvent::Trans | HmmEvent::MotifEnd => {} - } - } - base_labels.push(BaseLabel::MotifBound); - base_labels + Ok(betas) } diff --git a/src/trvz/locus.rs b/src/trvz/locus.rs index e9d43c4..9f1f8a0 100644 --- a/src/trvz/locus.rs +++ b/src/trvz/locus.rs @@ -14,21 +14,12 @@ pub struct Locus { pub region: GenomicRegion, } -#[derive(Debug, PartialEq, Clone)] -pub enum BaseLabel { - Match, - Mismatch, - NoMatch, - Skip, - MotifBound, -} - #[derive(Debug)] pub struct Allele { pub seq: String, pub region_labels: Vec, pub flank_labels: Vec, - pub base_labels: Vec, + //pub base_labels: Vec, } fn get_flanks( diff --git a/src/trvz/mod.rs b/src/trvz/mod.rs index 28fc785..316ef6c 100644 --- a/src/trvz/mod.rs +++ b/src/trvz/mod.rs @@ -1,11 +1,11 @@ pub mod align; -pub mod genotype_plot; +pub mod align_allele; +pub mod align_consensus; +pub mod align_reads; +pub mod allele_plot; +pub mod color; pub mod input; pub mod locus; -pub mod pdf; -pub mod pipe_plot; -pub mod png; pub mod read; pub mod struc; -pub mod svg; pub mod waterfall_plot; diff --git a/src/trvz/pipe_plot.rs b/src/trvz/pipe_plot.rs deleted file mode 100644 index ac990ae..0000000 --- a/src/trvz/pipe_plot.rs +++ /dev/null @@ -1,484 +0,0 @@ -use crate::trvz::{ - align::AlignInfo, - locus::{BaseLabel, Locus}, - struc::RegionLabel, -}; -use bio::alignment::AlignmentOperation; -use itertools::Itertools; - -pub type AlignOp = AlignmentOperation; - -#[derive(Debug, PartialEq, Clone)] -pub enum Color { - Purple, - Blue, - Orange, - Teal, - Gray, - LightGray, - Black, - Green, - Pink, - Yellow, - Red, - Khaki, - PaleRed, - PaleBlue, - Grad(f64), -} - -#[derive(Debug, PartialEq)] -pub enum Shape { - Rect, - HLine, - VLine, -} - -#[derive(Debug)] -pub struct PipeSeg { - pub width: u32, - pub color: Color, - pub shape: Shape, -} - -#[derive(Debug)] -pub struct Beta { - pub value: f64, - pub pos: usize, -} - -#[derive(Debug)] -pub struct Pipe { - pub segs: Vec, - pub betas: Vec, - pub height: u32, - pub outline: Vec, - pub scale: Vec<(u32, Option)>, -} - -pub struct Legend { - pub labels: Vec<(String, Color)>, - pub height: u32, -} - -#[derive(Debug)] -pub struct DisplayParams { - pub what_to_show: String, - pub max_width: usize, -} - -pub type PlotPanel = Vec; - -pub struct PipePlot { - pub panels: Vec, - pub legend: Legend, - // pub has_meth: bool, -} - -pub fn encode_color(color: &Color) -> String { - match color { - Color::Purple => "#814ED1".to_string(), - Color::Blue => "#1383C6".to_string(), - Color::Orange => "#E16A2C".to_string(), - Color::Teal => "#009CA2".to_string(), - Color::Gray => "#BABABA".to_string(), - Color::LightGray => "#D1D1D1".to_string(), - Color::Black => "#000000".to_string(), - Color::Pink => "#ED3981".to_string(), - Color::Yellow => "#EFCD17".to_string(), - Color::Green => "#009D4E".to_string(), - Color::Red => "#E3371E".to_string(), - Color::Khaki => "#F0E68C".to_string(), - Color::PaleRed => "#FF4858".to_string(), - Color::PaleBlue => "#46B2E8".to_string(), - Color::Grad(value) => get_gradient(*value), - } -} - -fn get_gradient(value: f64) -> String { - let blue: (u8, u8, u8) = (0, 73, 255); - let red: (u8, u8, u8) = (255, 0, 0); - let mix_red = (blue.0 as f64 * (1.0 - value) + red.0 as f64 * value).round() as u8; - let mix_green = (blue.1 as f64 * (1.0 - value) + red.1 as f64 * value).round() as u8; - let mix_blue = (blue.2 as f64 * (1.0 - value) + red.2 as f64 * value).round() as u8; - - format!("#{:02X}{:02X}{:02X}", mix_red, mix_green, mix_blue) -} - -/* -pub fn get_base_pipe(locus: &Locus, allele: &Allele) -> Pipe { - let mut segs = Vec::new(); - for label in &allele.all_labels { - let color = get_color(&locus, AlignOp::Match, label); - if let RegionLabel::Flank(start, end) = label { - segs.push(PipeSeg { - width: (*end - *start) as u32, - color: color, - shape: Shape::Rect, - }); - } else if let RegionLabel::Seq(start, end) = label { - segs.push(PipeSeg { - width: (*end - *start) as u32, - color: color, - shape: Shape::Rect, - }); - } else if let RegionLabel::Other(start, end) = label { - segs.push(PipeSeg { - width: (*end - *start) as u32, - color: Color::Gray, - shape: Shape::Rect, - }); - } else if let RegionLabel::Tr(start, end, motif) = label { - let tr_seq = &allele.seq[*start..*end]; - let mut index = 0; - let mut last_matched_index = 0; - while index as i32 <= tr_seq.len() as i32 - motif.len() as i32 { - let observed_motif = &tr_seq[index..index + motif.len()]; - if observed_motif == motif { - if last_matched_index != index { - segs.push(PipeSeg { - width: (index - last_matched_index) as u32, - color: Color::Gray, - shape: Shape::Rect, - }); - } - segs.push(PipeSeg { - width: motif.len() as u32, - color: color.clone(), - shape: Shape::Rect, - }); - index += motif.len(); - last_matched_index = index; - } else { - index += 1; - } - } - if last_matched_index != index { - segs.push(PipeSeg { - width: (index - last_matched_index) as u32, - color: Color::Gray, - shape: Shape::Rect, - }); - } - } else { - panic!(); - } - } - - let height = 3; - Pipe { - segs, - betas: Vec::new(), - height, - outline: Vec::new(), - scale: Vec::new(), - } -} -*/ - -pub fn get_allele_pipe( - locus: &Locus, - tick_spacing: usize, - spans: &[RegionLabel], - allele_labels: &Vec, -) -> Pipe { - let mut span_and_label = Vec::new(); - let mut base_pos = 0; - for label in allele_labels - .iter() - .filter(|l| **l != BaseLabel::MotifBound) - { - let span_index = get_label_index(spans, &AlignOp::Match, base_pos); - span_and_label.push((span_index, label)); - - base_pos += match label { - BaseLabel::Match => 1, - BaseLabel::Mismatch => 1, - BaseLabel::MotifBound => 0, - BaseLabel::NoMatch => 1, - BaseLabel::Skip => 0, - }; - } - - let mut segs = Vec::new(); - let groups = span_and_label.iter().group_by(|(a, b)| (a, b)); - for (key, group) in &groups { - let length = group.count(); - let (span_index, allele_label) = key; - - let op = match allele_label { - BaseLabel::Match => AlignOp::Match, - BaseLabel::Mismatch => AlignOp::Subst, - BaseLabel::MotifBound => AlignOp::Ins, - BaseLabel::Skip => AlignOp::Ins, - BaseLabel::NoMatch => AlignOp::Del, - }; - - let width = match op { - AlignOp::Ins => 0, - _ => length as u32, - }; - - let color = get_color(locus, op, &spans[*span_index]); - - let shape = match op { - AlignOp::Del => Shape::HLine, - AlignOp::Ins => Shape::VLine, - _ => Shape::Rect, - }; - - segs.push(PipeSeg { - width, - color, - shape, - }); - } - - let outline = generate_outline(spans); - let scale = generate_scale(tick_spacing, spans, allele_labels); - - Pipe { - segs, - betas: Vec::new(), - height: 3, - outline, - scale, - } -} - -fn generate_scale( - tick_spacing: usize, - regions: &[RegionLabel], - bases: &Vec, -) -> Vec<(u32, Option)> { - let label_spacing = tick_spacing * 5; - let mut bounds_by_tr = Vec::new(); - let mut allele_pos: u32 = 0; - - let mut motif_index = 0; - for base in bases { - if *base == BaseLabel::MotifBound { - let region_index = get_region_id(regions, allele_pos); - let is_first_motif = match ®ions[region_index] { - RegionLabel::Tr(start, _, _) => allele_pos == *start as u32, - _ => continue, - }; - - if is_first_motif { - motif_index = 0; - } else { - motif_index += 1; - } - - if motif_index % tick_spacing != 0 { - continue; - } - - let label = if motif_index % label_spacing == 0 { - Some(motif_index as u32) - } else { - None - }; - - bounds_by_tr.push((allele_pos, label)); - } - - allele_pos += match base { - BaseLabel::Match => 1, - BaseLabel::Mismatch => 1, - BaseLabel::NoMatch => 1, - _ => 0, - }; - } - - bounds_by_tr -} - -fn get_region_id(regions: &[RegionLabel], pos: u32) -> usize { - for (index, region) in regions.iter().enumerate() { - let (start, end) = match region { - RegionLabel::Flank(start, end) => (start, end), - RegionLabel::Other(start, end) => (start, end), - RegionLabel::Seq(start, end) => (start, end), - RegionLabel::Tr(start, end, _) => (start, end), - }; - - if *start as u32 <= pos && pos < *end as u32 { - return index; - } - } - - panic!() -} - -fn generate_outline(regions: &[RegionLabel]) -> Vec { - regions - .iter() - .map(|r| { - let width = match r { - RegionLabel::Flank(start, end) => end - start, - RegionLabel::Other(start, end) => end - start, - RegionLabel::Seq(start, end) => end - start, - RegionLabel::Tr(start, end, _) => end - start, - } as u32; - - PipeSeg { - width, - color: Color::Black, - shape: Shape::Rect, - } - }) - .collect_vec() -} - -pub fn get_pipe( - locus: &Locus, - labels: &[RegionLabel], - align: &AlignInfo, - show_betas: bool, -) -> Pipe { - assert!(align.align.xstart == 0); - assert!(align.align.ystart == 0); - - let mut segs = Vec::new(); - let mut betas = Vec::new(); - - let mut ops_and_locs = Vec::new(); - let mut allele_pos = 0; - let mut seq_pos = 0; - let mut cpg_index = 0; - for op in &align.align.operations { - let label_index = get_label_index(labels, op, allele_pos); - ops_and_locs.push((*op, label_index)); - - let is_cpg = seq_pos + 1 < align.seq.len() && &align.seq[seq_pos..seq_pos + 2] == "CG"; - - if show_betas && *op != AlignOp::Del && align.meth.is_some() && is_cpg { - if *op == AlignOp::Match || *op == AlignOp::Subst { - let meth = align.meth.as_ref().unwrap(); - let value = meth[cpg_index] as f64 / 255.0; - betas.push(Beta { - value, - pos: allele_pos, - }); - } - cpg_index += 1; - } - - allele_pos += match *op { - AlignOp::Match => 1, - AlignOp::Subst => 1, - AlignOp::Del => 1, - AlignOp::Ins => 0, - _ => panic!("Unhandled operation {:?}", *op), - }; - - seq_pos += match *op { - AlignOp::Match => 1, - AlignOp::Subst => 1, - AlignOp::Del => 0, - AlignOp::Ins => 1, - _ => panic!("Unhandled operation {:?}", *op), - }; - } - - assert_eq!(align.seq.len(), seq_pos); - - let groups = ops_and_locs.iter().group_by(|(a, l)| (a, l)); - for (key, group) in &groups { - let group = group.collect_vec(); - segs.push(get_seg(locus, labels, key, group)); - } - - let height = 3; - Pipe { - segs, - betas, - height, - outline: Vec::new(), - scale: Vec::new(), - } -} - -// NOTE: ref_pos is not needed because we are working with global alignment -fn get_seg( - locus: &Locus, - labels: &[RegionLabel], - op_and_seg_index: (&AlignOp, &usize), - group: Vec<&(AlignOp, usize)>, -) -> PipeSeg { - let (op, seg_index) = op_and_seg_index; - let width = match op { - AlignOp::Ins => 0, - _ => group.len() as u32, - }; - - let color = get_color(locus, *op, &labels[*seg_index]); - - let shape = match op { - AlignOp::Del => Shape::HLine, - AlignOp::Ins => Shape::VLine, - _ => Shape::Rect, - }; - - PipeSeg { - width, - color, - shape, - } -} - -fn get_label_index(labels: &[RegionLabel], op: &AlignOp, pos: usize) -> usize { - for (index, label) in labels.iter().enumerate() { - let (start, end) = match label { - RegionLabel::Flank(start, end) => (*start, *end), - RegionLabel::Seq(start, end) => (*start, *end), - RegionLabel::Tr(start, end, _) => (*start, *end), - RegionLabel::Other(start, end) => (*start, *end), - }; - - if start <= pos && pos < end { - return index; - } - - if pos == end && *op == AlignOp::Ins { - return index; - } - } - - println!("labels = {:?}, pos = {}", labels, pos); - panic!(); -} - -pub fn get_color(locus: &Locus, op: AlignmentOperation, label: &RegionLabel) -> Color { - let tr_colors = [ - Color::Blue, - Color::Purple, - Color::Orange, - Color::Pink, - Color::Yellow, - Color::Green, - Color::Red, - Color::Khaki, - Color::PaleRed, - Color::PaleBlue, - ]; - - if op == AlignOp::Subst { - return Color::Gray; - } - - if op == AlignOp::Ins { - return Color::Black; - } - - match label { - RegionLabel::Flank(_, _) => Color::Teal, - RegionLabel::Seq(_, _) => Color::LightGray, - RegionLabel::Tr(_, _, motif) => { - let index = locus.motifs.iter().position(|m| m == motif).unwrap(); - tr_colors[index % tr_colors.len()].clone() - } - RegionLabel::Other(_, _) => Color::LightGray, - } -} diff --git a/src/trvz/read.rs b/src/trvz/read.rs index bac9c9c..b0d0f9b 100644 --- a/src/trvz/read.rs +++ b/src/trvz/read.rs @@ -1,7 +1,68 @@ +use bio::alignment::Alignment as BioAlign; +use bio::alignment::AlignmentOperation as BioOp; + +#[derive(Clone)] pub struct Read { pub seq: String, pub left_flank: usize, pub right_flank: usize, pub allele: i32, - pub meth: Option>, + pub betas: Betas, +} + +pub type Betas = Vec; + +#[derive(Debug, Clone)] +pub struct Beta { + pub pos: usize, + pub value: f64, +} + +/// Convert betas from read coordinates to allele coordinates +/// according to the provided alignment between the read and +/// the allele consensus sequence +pub fn project_betas(bio_align: &BioAlign, betas: &Betas) -> Betas { + if betas.is_empty() { + return Vec::new(); + } + let mut ref_pos = 0; + let mut seq_pos = 0; + let mut beta_index = 0; + + let mut proj_betas = Vec::new(); + for op in &bio_align.operations { + let at_pos = betas[beta_index].pos == seq_pos; + let is_visible = *op == BioOp::Match || *op == BioOp::Subst; + if at_pos && is_visible { + let beta = Beta { + pos: ref_pos, + value: betas[beta_index].value, + }; + proj_betas.push(beta); + } + if at_pos { + beta_index += 1; + } + if beta_index == betas.len() { + break; + } + + seq_pos += match *op { + BioOp::Match => 1, + BioOp::Subst => 1, + BioOp::Del => 0, + BioOp::Ins => 1, + _ => panic!("Unhandled operation {:?}", *op), + }; + + ref_pos += match *op { + BioOp::Match => 1, + BioOp::Subst => 1, + BioOp::Del => 1, + BioOp::Ins => 0, + _ => panic!("Unhandled operation {op:?}"), + }; + } + + proj_betas } diff --git a/src/trvz/svg.rs b/src/trvz/svg.rs deleted file mode 100644 index 23369e0..0000000 --- a/src/trvz/svg.rs +++ /dev/null @@ -1,374 +0,0 @@ -use crate::trvz::pipe_plot::{encode_color, Color, Legend, Pipe, PipePlot, PlotPanel, Shape}; -use std::{fs::File, io::Write}; - -pub fn generate(genotype_plot: &PipePlot, path: &str) { - let start = (10.0, 14.0); - - let longest_pipe = get_longest_pipe(genotype_plot); - let x_scale = 750.0 / longest_pipe as f64; - - let scale = (x_scale, 3.0); - let base_spacing = 3.0; - let read_spacing = 2.0; - let allele_spacing = 20.0; - let legend_spacing = 4.0; - let file = File::create(path).unwrap(); - let mut generator = Generator { - start, - scale, - base_spacing, - read_spacing, - allele_spacing, - legend_spacing, - file, - }; - generator.generate(genotype_plot); -} - -fn get_longest_pipe(plot: &PipePlot) -> u32 { - let mut widths = Vec::new(); - for panel in &plot.panels { - for pipe in panel { - let mut width = 0; - for seg in &pipe.segs { - width += seg.width; - } - widths.push(width); - } - } - if widths.is_empty() { - 0 - } else { - *widths.iter().max().unwrap() - } -} - -struct Generator { - start: (f64, f64), - scale: (f64, f64), - base_spacing: f64, - read_spacing: f64, - allele_spacing: f64, - legend_spacing: f64, - file: File, -} - -impl Generator { - pub fn generate(&mut self, genotype_plot: &PipePlot) { - let (width, height) = self.get_dimensions(genotype_plot); - self.start_svg(width, height); - self.add_background(); - - let base_x = self.start.0; - let mut base_y = self.start.1; - for (index, allele_plot) in genotype_plot.panels.iter().enumerate() { - self.plot_allele(base_x, &mut base_y, allele_plot); - if index + 1 != genotype_plot.panels.len() { - base_y += self.allele_spacing; - } - } - - base_y += self.legend_spacing; - self.plot_legend(base_x, base_y, &genotype_plot.legend); - - self.end_svg(); - } - - fn plot_allele(&mut self, base_x: f64, base_y: &mut f64, allele_plot: &PlotPanel) { - for pipe in allele_plot { - self.plot_pipe(pipe, base_x, *base_y); - if !pipe.outline.is_empty() { - self.plot_outline(pipe, base_x, *base_y); - self.plot_scale(pipe, base_x, *base_y); - *base_y += self.to_y(pipe.height) + self.base_spacing; - } else { - *base_y += self.to_y(pipe.height) + self.read_spacing; - } - } - } - - fn plot_legend(&mut self, base_x: f64, base_y: f64, legend: &Legend) { - let height = self.to_y(legend.height); - let mut x = base_x; - for (label, color) in &legend.labels { - self.add_rect((x, base_y), (height, height), color, false); - x += height + 2.0; //self.to_y(1); - let point = format!("x=\"{}\" y=\"{}\"", x, base_y + height - 1.0); - let height = r#"font-size="14px""#; - let style = r#"font-family="monospace" font-weight="bold""#; - let line = format!("{}", point, style, height, label); - writeln!(self.file, "{}", line).unwrap(); - x += 5.0 * (2 * label.len() as u32 + 1) as f64; - } - } - - fn plot_pipe(&mut self, pipe: &Pipe, x: f64, y: f64) { - let pipe_height = self.to_y(pipe.height); - - // Plot the main pipe - let mut x_cur = x; - - self.add_left_cap((x_cur, y), pipe_height); - for seg in &pipe.segs { - let dims = (self.to_x(seg.width), pipe_height); - if seg.shape == Shape::Rect { - self.add_rect((x_cur, y), dims, &seg.color, true); - } - - if seg.shape == Shape::HLine { - self.add_hline((x_cur, y), dims, &seg.color, 1.5); - } - - x_cur += self.to_x(seg.width); - } - self.add_right_cap((x_cur, y), pipe_height); - - // Plot insertions - let mut x_cur = x; - for seg in &pipe.segs { - let dims = (self.to_x(seg.width), pipe_height); - - if seg.shape == Shape::VLine { - self.add_vline((x_cur, y), dims, &seg.color); - } - - x_cur += self.to_x(seg.width); - } - - // Plot betas - for beta in &pipe.betas { - let beta_x = x + self.to_x(beta.pos as u32); - let dims = (self.to_x(1), pipe_height); - self.add_rect((beta_x, y), dims, &Color::Grad(beta.value), false); - } - } - - fn plot_outline(&mut self, pipe: &Pipe, x: f64, y: f64) { - let height = self.to_y(pipe.height); - let mut width = 0; - for seg in &pipe.outline { - width += seg.width; - } - let width = self.to_x(width); - - // Draw segment boundaries - let mut bar_x = x + self.to_x(pipe.outline.first().unwrap().width); - for segment in pipe.outline.iter().skip(1) { - let point1 = format!("x1=\"{}\" y1=\"{}\"", bar_x, y); - let point2 = format!("x2=\"{}\" y2=\"{}\"", bar_x, y + height); - let style = r##"stroke="#000000" stroke-width="1.5""##; - let line = format!("", point1, point2, style); - writeln!(self.file, "{}", line).unwrap(); - bar_x += self.to_x(segment.width); - } - - // Draw top horizontal line - let point1 = format!("x1=\"{}\" y1=\"{}\"", x, y); - let point2 = format!("x2=\"{}\" y2=\"{}\"", x + width, y); - let style = r##"stroke="#000000" stroke-width="1.5""##; - let line = format!("", point1, point2, style); - writeln!(self.file, "{}", line).unwrap(); - - // Draw bottom horizontal line - let point1 = format!("x1=\"{}\" y1=\"{}\"", x, y + height); - let point2 = format!("x2=\"{}\" y2=\"{}\"", x + width, y + height); - let style = r##"stroke="#000000" stroke-width="1.5""##; - let line = format!("", point1, point2, style); - writeln!(self.file, "{}", line).unwrap(); - - // Draw left cap - let point1 = format!("{} {}", x, y); - let point2 = format!("{} {}", x - height / 2.0, y); - let point3 = format!("{} {}", x - height / 2.0, y + height); - let point4 = format!("{} {}", x, y + height); - - let path = format!("d=\"M {} C {}, {}, {}\"", point1, point2, point3, point4); - let style = r##"stroke="#000000" stroke-width="1.5" fill="none" opacity="0.9""##; - - let line = format!("", path, style); - writeln!(self.file, "{}", line).unwrap(); - - // Draw right cap - let point1 = format!("{} {}", x + width, y); - let point2 = format!("{} {}", x + width + height / 2.0, y); - let point3 = format!("{} {}", x + width + height / 2.0, y + height); - let point4 = format!("{} {}", x + width, y + height); - - let path = format!("d=\"M {} C {}, {}, {}\"", point1, point2, point3, point4); - let style = r##"stroke="#000000" stroke-width="1.5" fill="none" opacity="0.9""##; - - let line = format!("", path, style); - writeln!(self.file, "{}", line).unwrap(); - } - - fn plot_scale(&mut self, pipe: &Pipe, x: f64, y: f64) { - let height = self.to_y(pipe.height) / 2.0; - - // Draw segment boundaries - for (x_pos, label) in pipe.scale.iter() { - let tick_x = x + self.to_x(*x_pos); - let point1 = format!("x1=\"{}\" y1=\"{}\"", tick_x, y - height / 2.0); - let point2 = format!("x2=\"{}\" y2=\"{}\"", tick_x, y + height / 2.0); - let style = r##"stroke="#000000" stroke-width="1.5""##; - let line = format!("", point1, point2, style); - writeln!(self.file, "{}", line).unwrap(); - - if let Some(label) = *label { - let point = format!("x=\"{}\" y=\"{}\"", tick_x, y - height / 2.0 - 2.0); // 2.0 is padding - let height = r#"font-size="10px""#; - let style = r#"font-family="monospace" font-weight="bold" text-anchor="middle""#; - let line = format!("{}", point, style, height, label); - writeln!(self.file, "{}", line).unwrap(); - } - } - } - - fn add_rect(&mut self, pos: (f64, f64), dims: (f64, f64), color: &Color, add_highlight: bool) { - let (x, y) = pos; - let (w, h) = dims; - - let pos = format!("x=\"{}\" y=\"{}\"", x, y); - let dim = format!("height=\"{}\" width=\"{}\"", h, w); - - let color = encode_color(color); - let style = format!("fill=\"{}\" stroke=\"{}\" stroke-width=\"0\"", color, color); - - let rect = format!("", pos, dim, style); - writeln!(self.file, "{}", rect).unwrap(); - - if add_highlight { - let pos = format!("x=\"{}\" y=\"{}\"", x, y + h * 0.18); - let dim = format!("height=\"{}\" width=\"{}\"", h / 3.0, w); - let highlight = format!("", pos, dim); - writeln!(self.file, "{}", highlight).unwrap(); - } - } - - fn add_hline(&mut self, pos: (f64, f64), dims: (f64, f64), color: &Color, stroke: f64) { - let x1 = pos.0; - let x2 = pos.0 + dims.0; - let y1 = pos.1 + dims.1 / 2.0; - let y2 = y1; - - let x1y1 = format!("x1=\"{}\" y1=\"{}\"", x1, y1); - let x2y2 = format!("x2=\"{}\" y2=\"{}\"", x2, y2); - - let style = format!( - "stroke=\"{}\" stroke-width=\"{}\"", - encode_color(color), - stroke - ); - - let line = format!("", x1y1, x2y2, style); - writeln!(self.file, "{}", line).unwrap(); - } - - fn add_vline(&mut self, pos: (f64, f64), dims: (f64, f64), color: &Color) { - let x1 = pos.0; - let x2 = pos.0; - let y1 = pos.1; - let y2 = pos.1 + dims.1; - - let x1y1 = format!("x1=\"{}\" y1=\"{}\"", x1, y1); - let x2y2 = format!("x2=\"{}\" y2=\"{}\"", x2, y2); - - let stroke_width = 2.0_f64.min(self.to_x(1)); - let style = format!( - "stroke=\"{}\" stroke-width=\"{}\"", - encode_color(color), - stroke_width - ); - - let line = format!("", x1y1, x2y2, style); - writeln!(self.file, "{}", line).unwrap(); - } - - fn add_left_cap(&mut self, pos: (f64, f64), height: f64) { - let (x, y) = pos; - let point1 = format!("{} {}", x, y); - let point2 = format!("{} {}", x - height / 2.0, y); - let point3 = format!("{} {}", x - height / 2.0, y + height); - let point4 = format!("{} {}", x, y + height); - - let path = format!("d=\"M {} C {}, {}, {}\"", point1, point2, point3, point4); - let fill = r##"fill="#009CA2""##; - let opacity = r#"opacity="0.9""#; - - let line = format!("", path, fill, opacity); - writeln!(self.file, "{}", line).unwrap(); - } - - fn add_right_cap(&mut self, pos: (f64, f64), height: f64) { - let (x, y) = pos; - let point1 = format!("{} {}", x, y); - let point2 = format!("{} {}", x + height / 2.0, y); - let point3 = format!("{} {}", x + height / 2.0, y + height); - let point4 = format!("{} {}", x, y + height); - - let path = format!("d=\"M {} C {}, {}, {}\"", point1, point2, point3, point4); - let fill = "fill=\"#009CA2\""; - - let opacity = r#"opacity="0.9""#; - - let line = format!("", path, fill, opacity); - writeln!(self.file, "{}", line).unwrap(); - } - - fn get_dimensions(&self, genotype_plot: &PipePlot) -> (f64, f64) { - let mut widths = Vec::new(); - let mut height = self.start.1; - - for allele_plot in &genotype_plot.panels { - for (index, pipe) in allele_plot.iter().enumerate() { - if index == 0 { - height += self.to_y(pipe.height) + self.base_spacing; - } else { - height += self.to_y(pipe.height) + self.read_spacing; - } - - let mut pipe_width = 2.0 * self.start.0; - for seg in &pipe.segs { - pipe_width += self.to_x(seg.width); - } - widths.push(pipe_width); - } - } - - height += self.allele_spacing * (genotype_plot.panels.len() as f64 - 1.0); - height += - self.legend_spacing + self.to_y(genotype_plot.legend.height) + self.legend_spacing; - - let width = *widths - .iter() - .max_by(|a, b| a.partial_cmp(b).unwrap()) - .unwrap(); - (width, height) - } - - fn start_svg(&mut self, width: f64, height: f64) { - writeln!(self.file, r#""#).unwrap(); - let line = r#"", width, height).unwrap(); - } - - fn end_svg(&mut self) { - writeln!(self.file, "").unwrap(); - } - - fn add_background(&mut self) { - writeln!( - self.file, - r#""# - ) - .unwrap(); - } - - fn to_x(&self, x: u32) -> f64 { - x as f64 * self.scale.0 - } - - fn to_y(&self, y: u32) -> f64 { - y as f64 * self.scale.1 - } -} diff --git a/src/trvz/waterfall_plot.rs b/src/trvz/waterfall_plot.rs index be87834..fefd2d1 100644 --- a/src/trvz/waterfall_plot.rs +++ b/src/trvz/waterfall_plot.rs @@ -1,185 +1,345 @@ -use crate::trvz::{ - align::AlignInfo, - locus::Locus, - pipe_plot::{ - get_color, get_pipe, AlignOp, Beta, Color, DisplayParams, Legend, Pipe, PipePlot, PipeSeg, - Shape, - }, - struc::RegionLabel, -}; +use super::color::get_meth_colors; +use super::color::Color; +use super::color::ColorMap; +use super::read::project_betas; +use super::read::Betas; +use super::{align::Align, locus::Locus, read::Read}; +use crate::trvz::align::AlignOp; +use crate::trvz::align::AlignSeg; +use crate::trvz::align::SegType; +use crate::trvz::read::Beta; +use bio::alignment::pairwise::Aligner; +use bio::alignment::Alignment as BioAlign; +use bio::alignment::AlignmentOperation as BioOp; use itertools::Itertools; +use pipeplot::Band; +use pipeplot::Legend; +use pipeplot::Pipe; +use pipeplot::PipePlot; +use pipeplot::Seg; +use pipeplot::Shape; -pub fn plot_waterfall(locus: &Locus, params: &DisplayParams, aligns: &[AlignInfo]) -> PipePlot { - let flank_len = locus.left_flank.len(); - let mut aligns = aligns.iter().collect_vec(); - aligns.sort_by_key(|a| get_tr_len(flank_len, a)); +pub fn plot_waterfall( + locus: &Locus, + what_to_show: &str, + reads: &[Read], + colors: &ColorMap, +) -> PipePlot { + let reads = reads + .iter() + .sorted_by(|r1, r2| r1.seq.len().cmp(&r2.seq.len())) + .collect_vec(); - let max_tr_len = aligns + let longest_read = reads.iter().map(|r| r.seq.len()).max().unwrap(); + let reads = reads .iter() - .map(|a| get_tr_len(flank_len, a)) - .max() - .unwrap(); - let mut waterfall = Vec::new(); - for align in aligns { - let labels = if params.what_to_show == "motifs" { - label_motifs(locus, align) - } else { - label_flanks(locus, align) - }; - let pipe = get_pipe(locus, &labels, align, params.what_to_show == "meth"); - let pipe = add_deletion(max_tr_len, &labels, pipe); - waterfall.push(pipe); - } + .map(|r| align(locus, longest_read, r)) + .collect_vec(); - let mut labels = Vec::new(); - for motif in &locus.motifs { - let color = get_color(locus, AlignOp::Match, &RegionLabel::Tr(0, 0, motif.clone())); - labels.push((motif.clone(), color)); - } + plot(locus, what_to_show, &reads, colors) +} - let legend = Legend { labels, height: 4 }; - PipePlot { - panels: vec![waterfall], - legend, - } +fn align(locus: &Locus, longest_read: usize, read: &Read) -> (Align, Vec) { + let lf_ref = locus.left_flank.as_bytes(); + let rf_ref = locus.right_flank.as_bytes(); + let lf_read = read.seq[..lf_ref.len()].as_bytes(); + let rf_read = read.seq[read.seq.len() - locus.right_flank.len()..].as_bytes(); + let lf_bio_align = get_flank_align(lf_ref, lf_read); + let mut align = convert(&lf_bio_align, SegType::LeftFlank); + // Placeholder for TR alignment + let tr = &read.seq[locus.left_flank.len()..read.seq.len() - locus.right_flank.len()]; + align.extend(label_motifs(&locus.motifs, tr)); + // Add deletion that lines up right flanks + align.push(AlignSeg { + width: longest_read - read.seq.len(), + op: AlignOp::Del, + seg_type: SegType::RightFlank, + }); + let rf_bio_align = get_flank_align(rf_ref, rf_read); + align.extend(convert(&rf_bio_align, SegType::RightFlank)); + + let mut proj_betas = Vec::new(); + + let lf_betas = &read + .betas + .iter() + .filter(|beta| beta.pos < lf_read.len()) + .cloned() + .collect_vec(); + proj_betas.extend(project_betas(&lf_bio_align, lf_betas)); + + let tr_betas = &read + .betas + .iter() + .filter(|beta| lf_read.len() <= beta.pos && beta.pos < lf_read.len() + tr.len()) + .map(|beta| Beta { + pos: beta.pos - lf_read.len(), + value: beta.value, + }) + .collect_vec(); + proj_betas.extend(tr_betas.iter().map(|beta| Beta { + value: beta.value, + pos: beta.pos + lf_read.len(), + })); + + let rf_betas = &read + .betas + .iter() + .filter(|beta| lf_read.len() + tr.len() <= beta.pos) + .map(|beta| Beta { + pos: beta.pos - lf_read.len() - tr.len(), + value: beta.value, + }) + .collect_vec(); + + proj_betas.extend( + project_betas(&rf_bio_align, rf_betas) + .iter() + .map(|beta| Beta { + pos: beta.pos + lf_read.len() + tr.len() + longest_read - read.seq.len(), + value: beta.value, + }), + ); + + // let mut betas = project_betas(&lf_bio_align, lf_betas); + assert_eq!( + read.betas.len(), + lf_betas.len() + tr_betas.len() + rf_betas.len() + ); + + (align, proj_betas) } -fn label_motifs(locus: &Locus, align: &AlignInfo) -> Vec { - let mut motifs = locus.motifs.clone(); - motifs.sort_by_key(|b| std::cmp::Reverse(b.len())); +fn get_flank_align(ref_seq: &[u8], read_seq: &[u8]) -> BioAlign { + let likely_max_len = ref_seq.len() + 50; + let mut aligner = get_aligner(likely_max_len); + aligner.global(read_seq, ref_seq) +} - let mut labels = Vec::new(); - let flank_len = locus.left_flank.len(); +/// Convert a rust-bio alignment into an internal alignments +fn convert(bio_align: &BioAlign, seg_type: SegType) -> Align { + assert!(bio_align.xstart == 0); + assert!(bio_align.ystart == 0); - let tr_start = flank_len; - let tr_end = align.align.y_aln_len() - flank_len; - labels.push(RegionLabel::Flank(0, tr_start)); + let mut align = Vec::new(); + let mut ref_pos = 0; + let mut seq_pos = 0; - let tr_seq = &align.seq[tr_start..tr_end]; - let mut index = 0; - let mut motif_by_base = Vec::new(); - while index < tr_seq.len() { - let mut motif_found = false; - for (motif_index, motif) in motifs.iter().enumerate() { - if index + motif.len() > tr_seq.len() { - break; - } + for (bio_op, group) in &bio_align.operations.iter().group_by(|op| *op) { + let run_len = group.count(); - if &tr_seq[index..index + motif.len()] == motif { - let mut piece = vec![motif_index; motif.len()]; - motif_by_base.append(&mut piece); - motif_found = true; - index += motif.len(); - break; - } - } + let align_seg = match *bio_op { + BioOp::Match => AlignSeg { + width: run_len, + op: AlignOp::Match, + seg_type, + }, + BioOp::Subst => AlignSeg { + width: run_len, + op: AlignOp::Subst, + seg_type, + }, + BioOp::Del => AlignSeg { + width: run_len, + op: AlignOp::Del, + seg_type, + }, + BioOp::Ins => AlignSeg { + width: 0, + op: AlignOp::Ins, + seg_type, + }, + _ => panic!("No logic to handle {bio_op:?}"), + }; - if !motif_found { - motif_by_base.push(motifs.len()); - index += 1; - } - } + align.push(align_seg); - let mut tr_pos = tr_start; - let groups = motif_by_base.iter().group_by(|rec| *rec); - for (motif_index, group) in &groups { - let group = group.collect_vec(); - if *motif_index == motifs.len() { - labels.push(RegionLabel::Other(tr_pos, tr_pos + group.len())); - } else { - let motif = motifs[*motif_index].clone(); - labels.push(RegionLabel::Tr(tr_pos, tr_pos + group.len(), motif)); - } + ref_pos += match *bio_op { + BioOp::Match => run_len, + BioOp::Subst => run_len, + BioOp::Del => run_len, + BioOp::Ins => 0, + _ => panic!("Unhandled operation {:?}", *bio_op), + }; - tr_pos += group.len(); + seq_pos += match *bio_op { + BioOp::Match => run_len, + BioOp::Subst => run_len, + BioOp::Del => 0, + BioOp::Ins => run_len, + _ => panic!("Unhandled operation {:?}", *bio_op), + }; } - assert_eq!(tr_pos, tr_end); + assert_eq!(bio_align.x_aln_len(), seq_pos); + assert_eq!(bio_align.y_aln_len(), ref_pos); - labels.push(RegionLabel::Flank(tr_end, align.align.y_aln_len())); - labels + align } -fn label_flanks(locus: &Locus, align: &AlignInfo) -> Vec { - let flank_len = locus.left_flank.len(); +pub fn plot( + locus: &Locus, + what_to_show: &str, + reads: &[(Align, Vec)], + colors: &ColorMap, +) -> PipePlot { + let height = 4; + let xpos = 0; + let mut ypos = 0; + let mut pipes = Vec::new(); + for (align, betas) in reads { + let (colors, betas) = if what_to_show == "meth" { + (get_meth_colors(&locus.motifs), betas.clone()) + } else { + (colors.clone(), Vec::new()) + }; + let pipe = get_pipe(xpos, ypos, height, align, &betas, &colors); + pipes.push(pipe); + ypos += 5; + } + let mut labels = Vec::new(); - let tr_start = flank_len; - let tr_end = align.align.y_aln_len() - flank_len; - labels.push(RegionLabel::Flank(0, tr_start)); - labels.push(RegionLabel::Other(tr_start, tr_end)); - labels.push(RegionLabel::Flank(tr_end, align.align.y_aln_len())); - labels -} + if what_to_show == "motifs" { + for (index, motif) in locus.motifs.iter().enumerate() { + let color = colors.get(&SegType::Tr(index)).unwrap().to_string(); + labels.push((motif.clone(), color)); + } + } else { + labels = vec![ + ("Methylated".to_string(), Color::Grad(1.0).to_string()), + ("Unmethylated".to_string(), Color::Grad(0.0).to_string()), + ]; + } + ypos += 1; + + let legend = Legend { + xpos, + ypos, + height, + labels, + }; -fn get_tr_len(flank_len: usize, align: &AlignInfo) -> usize { - align.align.ylen - 2 * flank_len + PipePlot { pipes, legend } } -fn add_deletion(max_tr_len: usize, labels: &[RegionLabel], pipe: Pipe) -> Pipe { - let flank_len = match labels.first().unwrap() { - RegionLabel::Flank(start, end) => end - start, - RegionLabel::Other(start, end) => end - start, - RegionLabel::Seq(start, end) => end - start, - RegionLabel::Tr(start, end, _) => end - start, - }; +fn get_pipe( + xpos: u32, + ypos: u32, + height: u32, + align: &Align, + betas: &Betas, + colors: &ColorMap, +) -> Pipe { + let segs = align + .iter() + .map(|align_seg| { + let shape = match align_seg.op { + AlignOp::Del => Shape::HLine, + AlignOp::Ins => Shape::VLine, + AlignOp::Match | AlignOp::Subst => Shape::Rect, + }; + let color = match align_seg.op { + AlignOp::Match => colors.get(&align_seg.seg_type).unwrap(), + AlignOp::Subst => &Color::Gray, + _ => &Color::Black, + }; + Seg { + width: align_seg.width as u32, + color: color.to_string(), + shape, + } + }) + .collect(); - let tr_end = *match labels.last().unwrap() { - RegionLabel::Flank(start, _) => start, - RegionLabel::Other(start, _) => start, - RegionLabel::Seq(start, _) => start, - RegionLabel::Tr(start, _, _) => start, - }; + let bands = betas + .iter() + .map(|beta| Band { + pos: beta.pos as u32, + width: 2, + color: Color::Grad(beta.value).to_string(), + }) + .collect_vec(); - let tr_len = tr_end - flank_len; - let deletion_len = max_tr_len - tr_len; - if deletion_len == 0 { - return pipe; + Pipe { + xpos, + ypos, + height, + segs, + bands, + outline: false, } +} - let mut segs = Vec::new(); - let mut ref_pos = 0; - for seg in pipe.segs { - if ref_pos == tr_end { - segs.push(PipeSeg { - width: deletion_len as u32, - color: Color::Gray, - shape: Shape::HLine, - }); - } - ref_pos += seg.width as usize; +fn label_motifs(motifs: &[String], seq: &str) -> Align { + // Preserve motif order after sorting + let mut motifs = motifs.iter().enumerate().collect_vec(); + motifs.sort_by_key(|(_c, m)| std::cmp::Reverse(m.len())); - segs.push(seg); - } + let mut align = Vec::new(); + let mut pos = 0; + while pos != seq.len() { + let mut motif_found = false; + for (motif_index, motif) in &motifs { + if pos + motif.len() > seq.len() { + continue; + } - let mut betas = Vec::new(); - for beta in pipe.betas { - if beta.pos <= tr_end { - betas.push(beta); - } else { - betas.push(Beta { - value: beta.value, - pos: beta.pos + deletion_len, + if &seq[pos..pos + motif.len()] == *motif { + align.push(AlignSeg { + width: motif.len(), + op: AlignOp::Match, + seg_type: SegType::Tr(*motif_index), + }); + + motif_found = true; + pos += motif.len(); + break; + } + } + + if !motif_found { + align.push(AlignSeg { + width: 1, + op: AlignOp::Match, + seg_type: SegType::Tr(motifs.len()), }); + pos += 1; } } - Pipe { - segs, - betas, - height: pipe.height, - outline: Vec::new(), - scale: Vec::new(), + group(&align) +} + +// TODO: factor out as a generic alignment operation +fn group(align: &Align) -> Align { + let mut grouped_align = Vec::new(); + for ((op, seg_type), group) in &align.iter().group_by(|a| (a.op.clone(), a.seg_type)) { + let width = group.into_iter().map(|seg| seg.width).sum(); + grouped_align.push(AlignSeg { + width, + op, + seg_type, + }); } + grouped_align } -// add tests +type ScoreFunc = fn(u8, u8) -> i32; -#[cfg(test)] -mod tests { - #[test] - fn test_parent_origin_matrix_1() { - let mut motifs = ["AC", "GTG", "AAAAAA", "GGGGG"]; - // motifs.sort_by(|a, b| b.len().cmp(&a.len())); - motifs.sort_by_key(|b| std::cmp::Reverse(b.len())); +fn score(a: u8, b: u8) -> i32 { + if a == b { + 1i32 + } else { + -1i32 } } + +fn get_aligner<'a>(likely_max_len: usize) -> Aligner<&'a ScoreFunc> { + Aligner::with_capacity( + likely_max_len, + likely_max_len, + -5, + -1, + &(score as ScoreFunc), + ) +} diff --git a/src/utils/io_utils.rs b/src/utils/io_utils.rs index b0c8a4d..ac6775a 100644 --- a/src/utils/io_utils.rs +++ b/src/utils/io_utils.rs @@ -1,9 +1,11 @@ use crate::utils::Result; +use std::path::Path; -pub fn create_writer(output_prefix: &str, output_suffix: &str, f: F) -> Result +pub fn create_writer(output_prefix: &Path, output_suffix: &str, f: F) -> Result where F: FnOnce(&str) -> Result, { - let output_path = format!("{}.{}", output_prefix, output_suffix); - f(&output_path) + let mut output_path = output_prefix.to_path_buf().into_os_string(); + output_path.push(format!(".{output_suffix}")); + f(output_path.to_str().unwrap()) }