From 10b8133f62f9ba970a85a7a56efc58583468a58f Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Mon, 8 Apr 2024 14:48:11 +0200 Subject: [PATCH 01/22] map_to_nodes gets superkmer node idx --- include/buckets.hpp | 21 ++++++++++++++ include/builder/util.hpp | 5 +++- include/dictionary.cpp | 59 ++++++++++++++++++++++++++++++++++++++++ include/dictionary.hpp | 5 ++++ include/minimizers.hpp | 2 +- 5 files changed, 90 insertions(+), 2 deletions(-) diff --git a/include/buckets.hpp b/include/buckets.hpp index 7da1478..34d3805 100644 --- a/include/buckets.hpp +++ b/include/buckets.hpp @@ -103,6 +103,27 @@ struct buckets { } return lookup_result(); } + // superkmer annotation + lookup_result superkmer_id_to_kmer_id(uint64_t super_kmer_id, uint64_t k) const { + uint64_t offset = offsets.access(super_kmer_id); + auto [res, contig_end] = offset_to_id(offset, k); + return res; + } + lookup_result lookup_superkmer_start(uint64_t bucket_id, kmer_t target_kmer, kmer_t target_kmer_rc, + uint64_t k, uint64_t m) const { + auto [begin, end] = locate_bucket(bucket_id); + for (uint64_t super_kmer_id = begin; super_kmer_id != end; ++super_kmer_id) { + if(is_valid(lookup_in_super_kmer(super_kmer_id, target_kmer, k, m)) || + is_valid(lookup_in_super_kmer(super_kmer_id, target_kmer_rc, k, m))){ + uint64_t offset = offsets.access(super_kmer_id); + auto [res, contig_end] = offset_to_id(offset, k); + return res; + } + } + + return lookup_result(); + } + lookup_result lookup_canonical(uint64_t bucket_id, kmer_t target_kmer, kmer_t target_kmer_rc, uint64_t k, uint64_t m) const { diff --git a/include/builder/util.hpp b/include/builder/util.hpp index 54527f9..36f587b 100644 --- a/include/builder/util.hpp +++ b/include/builder/util.hpp @@ -242,7 +242,10 @@ struct minimizers_tuples { } std::string get_minimizers_filename() const { - assert(m_num_files_to_merge > 0); + //assert(m_num_files_to_merge > 0); + if(!(m_num_files_to_merge > 0)){ + throw std::invalid_argument("m_num_files_to_merge > 0"); + } if (m_num_files_to_merge == 1) return get_tmp_output_filename(0); std::stringstream filename; filename << m_tmp_dirname << "/sshash.tmp.run_" << m_run_identifier << ".minimizers.bin"; diff --git a/include/dictionary.cpp b/include/dictionary.cpp index ce5f209..dc0d220 100644 --- a/include/dictionary.cpp +++ b/include/dictionary.cpp @@ -2,6 +2,65 @@ namespace sshash { +uint64_t dictionary::kmer_to_superkmer_idx(std::string_view kmer) const { + // Plan: + // kmer -> get minimizer -> get bucket ->skew_index? -> maybe move the following code to inside m_buckets... + // ->get superkmer range in offsets!!-> + // ->get closest (previous) superkmer coord in offsets + // (+check everything makes sense on contig level) + // reverse complement considerations: primary/canonical graph -> find smaller minimizer between kmer and rc_kmer + // skew index considerations???? + + // from dictionary::lookup_uint_canonical_parsing: + char const* kmer_str_ptr = kmer.data(); + kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str_ptr, m_k); + kmer_t uint_kmer_rc = util::compute_reverse_complement(uint_kmer, m_k); + + uint64_t minimizer = util::compute_minimizer(uint_kmer, m_k, m_m, m_seed); + uint64_t minimizer_rc = util::compute_minimizer(uint_kmer_rc, m_k, m_m, m_seed); + uint64_t bucket_id = m_minimizers.lookup(std::min(minimizer, minimizer_rc)); + + lookup_result res_superkmer; + // no skew index + if (m_skew_index.empty()){ + lookup_result res = m_buckets.lookup_superkmer_start(bucket_id, uint_kmer,uint_kmer_rc, m_k, m_m); + return res.kmer_id; + } + + auto [begin, end] = m_buckets.locate_bucket(bucket_id); + uint64_t num_super_kmers_in_bucket = end - begin; + uint64_t log2_bucket_size = util::ceil_log2_uint32(num_super_kmers_in_bucket); + // superkmer in skew index + if (log2_bucket_size > m_skew_index.min_log2) { + uint64_t pos = m_skew_index.lookup(uint_kmer, log2_bucket_size); + + /* It must hold pos < num_super_kmers_in_bucket for the kmer to exist. */ + if (pos < num_super_kmers_in_bucket) { + uint64_t superkmer_id = begin + pos; + lookup_result res = m_buckets.superkmer_id_to_kmer_id(superkmer_id, m_k); + return res.kmer_id; + } + + uint64_t pos_rc = m_skew_index.lookup(uint_kmer_rc, log2_bucket_size); + + /* It must hold pos < num_super_kmers_in_bucket for the kmer to exist. */ + if (pos_rc < num_super_kmers_in_bucket) { + uint64_t superkmer_id_rc = begin + pos; + lookup_result res = m_buckets.superkmer_id_to_kmer_id(superkmer_id_rc, m_k); + return res.kmer_id; + } + + return constants::invalid_uint64; + + } + + // superkmer not in skew index + lookup_result res = m_buckets.lookup_superkmer_start(bucket_id, uint_kmer,uint_kmer_rc, m_k, m_m); + return res.kmer_id; +} + + + lookup_result dictionary::lookup_uint_regular_parsing(kmer_t uint_kmer) const { uint64_t minimizer = util::compute_minimizer(uint_kmer, m_k, m_m, m_seed); uint64_t bucket_id = m_minimizers.lookup(minimizer); diff --git a/include/dictionary.hpp b/include/dictionary.hpp index 0bf0b6e..60e1a36 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -22,6 +22,7 @@ struct dictionary { uint64_t k() const { return m_k; } uint64_t m() const { return m_m; } uint64_t num_contigs() const { return m_buckets.pieces.size() - 1; } + uint64_t num_superkmers() const { return m_buckets.offsets.size(); } bool canonicalized() const { return m_canonical_parsing; } bool weighted() const { return !m_weights.empty(); } @@ -92,6 +93,8 @@ struct dictionary { void print_space_breakdown() const; void compute_statistics() const; + void superkmer_statistics() const; + template void visit(Visitor& visitor) { visitor.visit(m_size); @@ -105,6 +108,8 @@ struct dictionary { visitor.visit(m_weights); } + uint64_t kmer_to_superkmer_idx(std::string_view kmer) const ; + private: uint64_t m_size; uint64_t m_seed; diff --git a/include/minimizers.hpp b/include/minimizers.hpp index 15e4475..8a9379b 100644 --- a/include/minimizers.hpp +++ b/include/minimizers.hpp @@ -15,7 +15,7 @@ struct minimizers { mphf_config.verbose_output = false; mphf_config.num_threads = 1; uint64_t num_threads = std::thread::hardware_concurrency() >= 8 ? 8 : 1; - if (size >= num_threads) mphf_config.num_threads = num_threads; + if (size >= 10*num_threads) mphf_config.num_threads = num_threads; if (build_config.verbose) { std::cout << "building minimizers MPHF (PTHash) with " << mphf_config.num_threads From 74dd50a087347507dcf12febdf99a4a2d0d0d9a5 Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Mon, 8 Apr 2024 17:06:08 +0200 Subject: [PATCH 02/22] compilation not native --- CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5e92a0d..e3e5601 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -17,8 +17,8 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}) MESSAGE(STATUS "Compiling for processor: " ${CMAKE_HOST_SYSTEM_PROCESSOR}) if (UNIX AND (CMAKE_HOST_SYSTEM_PROCESSOR STREQUAL "x86_64")) - MESSAGE(STATUS "Compiling with flags: -march=native -mbmi2 -msse4.2") - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native") + MESSAGE(STATUS "Compiling with flags: -march=core-avx2 -mbmi2 -msse4.2") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=core-avx2") # Flags for PTHash: set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mbmi2 -msse4.2") # for hardware popcount and pdep endif() From abed673faa1ac73a1a93cedef9dd730445511448 Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Thu, 11 Apr 2024 11:16:02 +0200 Subject: [PATCH 03/22] fixed minimizer bug --- include/buckets.hpp | 15 +++++----- include/dictionary.cpp | 66 +++++++++++++++--------------------------- include/dictionary.hpp | 4 +-- 3 files changed, 33 insertions(+), 52 deletions(-) diff --git a/include/buckets.hpp b/include/buckets.hpp index 34d3805..3976d50 100644 --- a/include/buckets.hpp +++ b/include/buckets.hpp @@ -109,18 +109,17 @@ struct buckets { auto [res, contig_end] = offset_to_id(offset, k); return res; } - lookup_result lookup_superkmer_start(uint64_t bucket_id, kmer_t target_kmer, kmer_t target_kmer_rc, + + + lookup_result lookup_superkmer_start(uint64_t begin, uint64_t end, kmer_t target_kmer, uint64_t k, uint64_t m) const { - auto [begin, end] = locate_bucket(bucket_id); for (uint64_t super_kmer_id = begin; super_kmer_id != end; ++super_kmer_id) { - if(is_valid(lookup_in_super_kmer(super_kmer_id, target_kmer, k, m)) || - is_valid(lookup_in_super_kmer(super_kmer_id, target_kmer_rc, k, m))){ - uint64_t offset = offsets.access(super_kmer_id); - auto [res, contig_end] = offset_to_id(offset, k); - return res; + auto res = lookup_in_super_kmer(super_kmer_id, target_kmer, k, m); + if(res.kmer_id != constants::invalid_uint64){ + assert(is_valid(res)); + return superkmer_id_to_kmer_id(super_kmer_id, k); } } - return lookup_result(); } diff --git a/include/dictionary.cpp b/include/dictionary.cpp index dc0d220..1825d8d 100644 --- a/include/dictionary.cpp +++ b/include/dictionary.cpp @@ -1,65 +1,47 @@ #include "dictionary.hpp" namespace sshash { - -uint64_t dictionary::kmer_to_superkmer_idx(std::string_view kmer) const { - // Plan: - // kmer -> get minimizer -> get bucket ->skew_index? -> maybe move the following code to inside m_buckets... - // ->get superkmer range in offsets!!-> - // ->get closest (previous) superkmer coord in offsets - // (+check everything makes sense on contig level) - // reverse complement considerations: primary/canonical graph -> find smaller minimizer between kmer and rc_kmer - // skew index considerations???? - - // from dictionary::lookup_uint_canonical_parsing: - char const* kmer_str_ptr = kmer.data(); - kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str_ptr, m_k); - kmer_t uint_kmer_rc = util::compute_reverse_complement(uint_kmer, m_k); - +lookup_result dictionary::kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const { uint64_t minimizer = util::compute_minimizer(uint_kmer, m_k, m_m, m_seed); - uint64_t minimizer_rc = util::compute_minimizer(uint_kmer_rc, m_k, m_m, m_seed); - uint64_t bucket_id = m_minimizers.lookup(std::min(minimizer, minimizer_rc)); + uint64_t bucket_id = m_minimizers.lookup(minimizer); + auto [begin, end] = m_buckets.locate_bucket(bucket_id); - lookup_result res_superkmer; - // no skew index - if (m_skew_index.empty()){ - lookup_result res = m_buckets.lookup_superkmer_start(bucket_id, uint_kmer,uint_kmer_rc, m_k, m_m); - return res.kmer_id; - } + if (m_skew_index.empty()) return m_buckets.lookup_superkmer_start(begin, end, uint_kmer, m_k, m_m); - auto [begin, end] = m_buckets.locate_bucket(bucket_id); uint64_t num_super_kmers_in_bucket = end - begin; uint64_t log2_bucket_size = util::ceil_log2_uint32(num_super_kmers_in_bucket); // superkmer in skew index if (log2_bucket_size > m_skew_index.min_log2) { uint64_t pos = m_skew_index.lookup(uint_kmer, log2_bucket_size); - /* It must hold pos < num_super_kmers_in_bucket for the kmer to exist. */ if (pos < num_super_kmers_in_bucket) { uint64_t superkmer_id = begin + pos; lookup_result res = m_buckets.superkmer_id_to_kmer_id(superkmer_id, m_k); - return res.kmer_id; - } - - uint64_t pos_rc = m_skew_index.lookup(uint_kmer_rc, log2_bucket_size); - - /* It must hold pos < num_super_kmers_in_bucket for the kmer to exist. */ - if (pos_rc < num_super_kmers_in_bucket) { - uint64_t superkmer_id_rc = begin + pos; - lookup_result res = m_buckets.superkmer_id_to_kmer_id(superkmer_id_rc, m_k); - return res.kmer_id; + return res; } + return lookup_result(); + } + // superkmer not in skew index + return m_buckets.lookup_superkmer_start(begin, end, uint_kmer, m_k, m_m); - return constants::invalid_uint64; - - } +} +uint64_t dictionary::kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const { + // Plan: + // kmer -> get minimizer -> get bucket ->skew_index? -> + // ->get superkmer range in offsets!!-> + // ->get superkmer coords in offsets + // (+check everything makes sense on contig level) + // reverse complement considerations: primary/canonical graph -> find smaller minimizer between kmer and rc_kmer - // superkmer not in skew index - lookup_result res = m_buckets.lookup_superkmer_start(bucket_id, uint_kmer,uint_kmer_rc, m_k, m_m); + kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); + lookup_result res = kmer_to_superkmer_idx_helper(uint_kmer); + if(res.kmer_id == constants::invalid_uint64 && check_reverse_complement){ + kmer_t uint_kmer_rc = util::compute_reverse_complement(uint_kmer, m_k); + res = kmer_to_superkmer_idx_helper(uint_kmer_rc); + } return res.kmer_id; -} - +} lookup_result dictionary::lookup_uint_regular_parsing(kmer_t uint_kmer) const { uint64_t minimizer = util::compute_minimizer(uint_kmer, m_k, m_m, m_seed); diff --git a/include/dictionary.hpp b/include/dictionary.hpp index 60e1a36..1caf64e 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -107,8 +107,8 @@ struct dictionary { visitor.visit(m_skew_index); visitor.visit(m_weights); } - - uint64_t kmer_to_superkmer_idx(std::string_view kmer) const ; + lookup_result kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const ; + uint64_t kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const ; private: uint64_t m_size; From ee9c468ce08e918c6f50fafde51a221a7d73b7b9 Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Fri, 12 Apr 2024 16:43:18 +0200 Subject: [PATCH 04/22] fixed bug in skew index case --- include/dictionary.cpp | 46 +++++++++++++++++++++++------------------- 1 file changed, 25 insertions(+), 21 deletions(-) diff --git a/include/dictionary.cpp b/include/dictionary.cpp index 1825d8d..4b4691a 100644 --- a/include/dictionary.cpp +++ b/include/dictionary.cpp @@ -1,6 +1,8 @@ #include "dictionary.hpp" namespace sshash { + +////////////////////////////////////////// lookup_result dictionary::kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const { uint64_t minimizer = util::compute_minimizer(uint_kmer, m_k, m_m, m_seed); uint64_t bucket_id = m_minimizers.lookup(minimizer); @@ -15,9 +17,10 @@ lookup_result dictionary::kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const { uint64_t pos = m_skew_index.lookup(uint_kmer, log2_bucket_size); /* It must hold pos < num_super_kmers_in_bucket for the kmer to exist. */ if (pos < num_super_kmers_in_bucket) { - uint64_t superkmer_id = begin + pos; - lookup_result res = m_buckets.superkmer_id_to_kmer_id(superkmer_id, m_k); - return res; + if(m_buckets.lookup_in_super_kmer(begin + pos, uint_kmer, m_k, m_m).kmer_id != sshash::constants::invalid_uint64){ + uint64_t superkmer_id = begin + pos; + return m_buckets.superkmer_id_to_kmer_id(superkmer_id, m_k); + } } return lookup_result(); } @@ -25,24 +28,6 @@ lookup_result dictionary::kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const { return m_buckets.lookup_superkmer_start(begin, end, uint_kmer, m_k, m_m); } -uint64_t dictionary::kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const { - // Plan: - // kmer -> get minimizer -> get bucket ->skew_index? -> - // ->get superkmer range in offsets!!-> - // ->get superkmer coords in offsets - // (+check everything makes sense on contig level) - // reverse complement considerations: primary/canonical graph -> find smaller minimizer between kmer and rc_kmer - - kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); - lookup_result res = kmer_to_superkmer_idx_helper(uint_kmer); - if(res.kmer_id == constants::invalid_uint64 && check_reverse_complement){ - kmer_t uint_kmer_rc = util::compute_reverse_complement(uint_kmer, m_k); - res = kmer_to_superkmer_idx_helper(uint_kmer_rc); - } - return res.kmer_id; - -} - lookup_result dictionary::lookup_uint_regular_parsing(kmer_t uint_kmer) const { uint64_t minimizer = util::compute_minimizer(uint_kmer, m_k, m_m, m_seed); uint64_t bucket_id = m_minimizers.lookup(minimizer); @@ -110,6 +95,25 @@ lookup_result dictionary::lookup_advanced(char const* string_kmer, kmer_t uint_kmer = util::string_to_uint_kmer(string_kmer, m_k); return lookup_advanced_uint(uint_kmer, check_reverse_complement); } +////////////////////////////////////////// +uint64_t dictionary::kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const { + // Plan: + // kmer -> get minimizer -> get bucket ->skew_index? -> + // ->get superkmer range in offsets!!-> + // ->get superkmer coords in offsets + // (+check everything makes sense on contig level) + // reverse complement considerations: primary graph -> look for both minimizers of kmer and rc_kmer + + kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); + lookup_result res = kmer_to_superkmer_idx_helper(uint_kmer); + assert(res.kmer_orientation == constants::forward_orientation); + if(res.kmer_id == constants::invalid_uint64 && check_reverse_complement){ + kmer_t uint_kmer_rc = util::compute_reverse_complement(uint_kmer, m_k); + res = kmer_to_superkmer_idx_helper(uint_kmer_rc); + res.kmer_orientation = constants::backward_orientation; + } + return res.kmer_id; +} lookup_result dictionary::lookup_advanced_uint(kmer_t uint_kmer, bool check_reverse_complement) const { if (m_canonical_parsing) return lookup_uint_canonical_parsing(uint_kmer); From a7d4a79b1b86dcb23afc0c3e3a32712b7ab02b75 Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Thu, 18 Apr 2024 11:02:50 +0200 Subject: [PATCH 05/22] changes for loss-less annotation: added superkmer mask and return superkmer id with superkmer start index --- include/buckets.hpp | 6 ++-- include/dictionary.cpp | 16 +++++---- include/dictionary.hpp | 6 ++-- include/statistics.cpp | 76 +++++++++++++++++++++++++++++++++++++++++- 4 files changed, 91 insertions(+), 13 deletions(-) diff --git a/include/buckets.hpp b/include/buckets.hpp index 3976d50..3359e21 100644 --- a/include/buckets.hpp +++ b/include/buckets.hpp @@ -111,16 +111,16 @@ struct buckets { } - lookup_result lookup_superkmer_start(uint64_t begin, uint64_t end, kmer_t target_kmer, + std::pair lookup_superkmer_start(uint64_t begin, uint64_t end, kmer_t target_kmer, uint64_t k, uint64_t m) const { for (uint64_t super_kmer_id = begin; super_kmer_id != end; ++super_kmer_id) { auto res = lookup_in_super_kmer(super_kmer_id, target_kmer, k, m); if(res.kmer_id != constants::invalid_uint64){ assert(is_valid(res)); - return superkmer_id_to_kmer_id(super_kmer_id, k); + return std::pair(superkmer_id_to_kmer_id(super_kmer_id, k), super_kmer_id); // CHECK THIS IS CORRECT!?!?! YES } } - return lookup_result(); + return std::pair(lookup_result(), constants::invalid_uint64); } diff --git a/include/dictionary.cpp b/include/dictionary.cpp index 4b4691a..5a8bbe2 100644 --- a/include/dictionary.cpp +++ b/include/dictionary.cpp @@ -3,7 +3,7 @@ namespace sshash { ////////////////////////////////////////// -lookup_result dictionary::kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const { +std::pair dictionary::kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const { uint64_t minimizer = util::compute_minimizer(uint_kmer, m_k, m_m, m_seed); uint64_t bucket_id = m_minimizers.lookup(minimizer); auto [begin, end] = m_buckets.locate_bucket(bucket_id); @@ -19,10 +19,10 @@ lookup_result dictionary::kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const { if (pos < num_super_kmers_in_bucket) { if(m_buckets.lookup_in_super_kmer(begin + pos, uint_kmer, m_k, m_m).kmer_id != sshash::constants::invalid_uint64){ uint64_t superkmer_id = begin + pos; - return m_buckets.superkmer_id_to_kmer_id(superkmer_id, m_k); + return std::pair(m_buckets.superkmer_id_to_kmer_id(superkmer_id, m_k), superkmer_id); // CHECK THIS IS CORRECT!?!?! YES } } - return lookup_result(); + return std::pair(lookup_result(), constants::invalid_uint64); } // superkmer not in skew index return m_buckets.lookup_superkmer_start(begin, end, uint_kmer, m_k, m_m); @@ -96,7 +96,7 @@ lookup_result dictionary::lookup_advanced(char const* string_kmer, return lookup_advanced_uint(uint_kmer, check_reverse_complement); } ////////////////////////////////////////// -uint64_t dictionary::kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const { +std::pair dictionary::kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const { // Plan: // kmer -> get minimizer -> get bucket ->skew_index? -> // ->get superkmer range in offsets!!-> @@ -105,14 +105,16 @@ uint64_t dictionary::kmer_to_superkmer_idx(char const* kmer_str, bool check_reve // reverse complement considerations: primary graph -> look for both minimizers of kmer and rc_kmer kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); - lookup_result res = kmer_to_superkmer_idx_helper(uint_kmer); + auto[res, s_idx] = kmer_to_superkmer_idx_helper(uint_kmer); assert(res.kmer_orientation == constants::forward_orientation); if(res.kmer_id == constants::invalid_uint64 && check_reverse_complement){ kmer_t uint_kmer_rc = util::compute_reverse_complement(uint_kmer, m_k); - res = kmer_to_superkmer_idx_helper(uint_kmer_rc); + auto pair = kmer_to_superkmer_idx_helper(uint_kmer_rc); + res = pair.first; + s_idx = pair.second; res.kmer_orientation = constants::backward_orientation; } - return res.kmer_id; + return std::pair(res.kmer_id, s_idx); } lookup_result dictionary::lookup_advanced_uint(kmer_t uint_kmer, bool check_reverse_complement) const { diff --git a/include/dictionary.hpp b/include/dictionary.hpp index 1caf64e..b1c3dda 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -95,6 +95,8 @@ struct dictionary { void superkmer_statistics() const; + std::vector build_superkmer_bv(const std::function (std::string)> &get_annotation_labels) const; + template void visit(Visitor& visitor) { visitor.visit(m_size); @@ -107,8 +109,8 @@ struct dictionary { visitor.visit(m_skew_index); visitor.visit(m_weights); } - lookup_result kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const ; - uint64_t kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const ; + std::pair kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const ; + std::pair kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const ; private: uint64_t m_size; diff --git a/include/statistics.cpp b/include/statistics.cpp index 3b77d67..73eac05 100644 --- a/include/statistics.cpp +++ b/include/statistics.cpp @@ -1,6 +1,5 @@ #include "dictionary.hpp" #include "buckets_statistics.hpp" - namespace sshash { void dictionary::compute_statistics() const { @@ -47,5 +46,80 @@ void dictionary::compute_statistics() const { buckets_stats.print_full(); std::cout << "DONE" << std::endl; } +bool equal(const std::vector& input1, const std::vector& input2) { + if(input1.size() != input2.size()){ + return false; + } + for(size_t i = 0; i < input1.size(); i++){ + if(input1.at(i) != input2.at(i)){ + return false; + } + } + return true; +} + +std::vector dictionary::build_superkmer_bv(const std::function (std::string)> &get_annotation_labels) const { + uint64_t num_kmers = size(); + uint64_t num_minimizers = m_minimizers.size(); + uint64_t num_super_kmers = m_buckets.offsets.size(); + + buckets_statistics buckets_stats(num_minimizers, num_kmers, num_super_kmers); + + std::cout << "building super kmer mask..." << std::endl; + std::vector non_mono_superkmer (num_super_kmers, false); + size_t superkmer_idx = 0; + + std::cout<<"iterating through buckets :\n"; + for (uint64_t bucket_id = 0; bucket_id != num_minimizers; ++bucket_id) { + std::cout << int(100*bucket_id/num_minimizers) <<"%" << '\r'<< std::flush; + auto [begin, end] = m_buckets.locate_bucket(bucket_id); + uint64_t num_super_kmers_in_bucket = end - begin; + buckets_stats.add_num_super_kmers_in_bucket(num_super_kmers_in_bucket); + for (uint64_t super_kmer_id = begin; super_kmer_id != end; ++super_kmer_id, ++superkmer_idx) { + //std::cout << super_kmer_id <<" "; + uint64_t offset = m_buckets.offsets.access(super_kmer_id); + auto [_, contig_end] = m_buckets.offset_to_id(offset, m_k); + (void)_; + bit_vector_iterator bv_it(m_buckets.strings, 2 * offset); + uint64_t window_size = std::min(m_k - m_m + 1, contig_end - offset - m_k + 1); + uint64_t prev_minimizer = constants::invalid_uint64; + std::vector prev_labels = {}; + uint64_t w = 0; + for (; w != window_size; ++w) { + uint64_t kmer = bv_it.read_and_advance_by_two(2 * m_k); + auto [minimizer, pos] = util::compute_minimizer_pos(kmer, m_k, m_m, m_seed); + if (m_canonical_parsing) { + uint64_t kmer_rc = util::compute_reverse_complement(kmer, m_k); + auto [minimizer_rc, pos_rc] = + util::compute_minimizer_pos(kmer_rc, m_k, m_m, m_seed); + if (minimizer_rc < minimizer) { + minimizer = minimizer_rc; + pos = pos_rc; + } + } + + std::string kmer_str = ""; + util::uint_kmer_to_string(kmer, &kmer_str[0], m_k); + std::vector labels = get_annotation_labels(kmer_str); + if(prev_labels.size() == 0){ + prev_labels = labels; + } else if(non_mono_superkmer[superkmer_idx] == false && !equal(labels, prev_labels)){ + non_mono_superkmer[superkmer_idx] = 1; + prev_labels = labels; + } + if (prev_minimizer != constants::invalid_uint64 and minimizer != prev_minimizer) { + break; + } + prev_minimizer = minimizer; + } + buckets_stats.add_num_kmers_in_super_kmer(num_super_kmers_in_bucket, w); + } + } + std::cout << "DONE" << std::endl; + std::cout<< "numbers of superkmers: " << superkmer_idx << std::endl; + + buckets_stats.print_full(); + return non_mono_superkmer; +} } // namespace sshash From c008c20d9ab1fc513f498af731dcb7dc91eba38e Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Mon, 22 Apr 2024 09:59:18 +0200 Subject: [PATCH 06/22] fixed empty bitvec bug --- include/dictionary.hpp | 2 +- include/statistics.cpp | 38 ++++++++++++++++++-------------------- 2 files changed, 19 insertions(+), 21 deletions(-) diff --git a/include/dictionary.hpp b/include/dictionary.hpp index b1c3dda..3a79159 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -95,7 +95,7 @@ struct dictionary { void superkmer_statistics() const; - std::vector build_superkmer_bv(const std::function (std::string)> &get_annotation_labels) const; + std::vector build_superkmer_bv(const std::function (std::string_view)> &get_annotation_labels) const; template void visit(Visitor& visitor) { diff --git a/include/statistics.cpp b/include/statistics.cpp index 73eac05..becb672 100644 --- a/include/statistics.cpp +++ b/include/statistics.cpp @@ -58,25 +58,22 @@ bool equal(const std::vector& input1, const std::vector dictionary::build_superkmer_bv(const std::function (std::string)> &get_annotation_labels) const { +std::vector dictionary::build_superkmer_bv(const std::function (std::string_view)> &get_annotation_labels) const { uint64_t num_kmers = size(); uint64_t num_minimizers = m_minimizers.size(); uint64_t num_super_kmers = m_buckets.offsets.size(); - buckets_statistics buckets_stats(num_minimizers, num_kmers, num_super_kmers); - std::cout << "building super kmer mask..." << std::endl; std::vector non_mono_superkmer (num_super_kmers, false); - size_t superkmer_idx = 0; + //size_t superkmer_idx = 0; + uint64_t one_pc_buckets = num_minimizers/100; std::cout<<"iterating through buckets :\n"; for (uint64_t bucket_id = 0; bucket_id != num_minimizers; ++bucket_id) { - std::cout << int(100*bucket_id/num_minimizers) <<"%" << '\r'<< std::flush; + if(bucket_id%one_pc_buckets==0)std::cout << bucket_id/one_pc_buckets <<"%" << '\r'<< std::flush; auto [begin, end] = m_buckets.locate_bucket(bucket_id); uint64_t num_super_kmers_in_bucket = end - begin; - buckets_stats.add_num_super_kmers_in_bucket(num_super_kmers_in_bucket); - for (uint64_t super_kmer_id = begin; super_kmer_id != end; ++super_kmer_id, ++superkmer_idx) { - //std::cout << super_kmer_id <<" "; + for (uint64_t super_kmer_id = begin; super_kmer_id != end; ++super_kmer_id) { uint64_t offset = m_buckets.offsets.access(super_kmer_id); auto [_, contig_end] = m_buckets.offset_to_id(offset, m_k); (void)_; @@ -86,7 +83,7 @@ std::vector dictionary::build_superkmer_bv(const std::function prev_labels = {}; uint64_t w = 0; for (; w != window_size; ++w) { - uint64_t kmer = bv_it.read_and_advance_by_two(2 * m_k); + kmer_t kmer = bv_it.read_and_advance_by_two(2 * m_k); auto [minimizer, pos] = util::compute_minimizer_pos(kmer, m_k, m_m, m_seed); if (m_canonical_parsing) { uint64_t kmer_rc = util::compute_reverse_complement(kmer, m_k); @@ -98,27 +95,28 @@ std::vector dictionary::build_superkmer_bv(const std::function labels = get_annotation_labels(kmer_str); - if(prev_labels.size() == 0){ + std::vector labels = get_annotation_labels(kmer_str); + if(prev_labels.size() == 0){ prev_labels = labels; } else if(non_mono_superkmer[superkmer_idx] == false && !equal(labels, prev_labels)){ - non_mono_superkmer[superkmer_idx] = 1; + non_mono_superkmer[superkmer_idx] = true; prev_labels = labels; } - if (prev_minimizer != constants::invalid_uint64 and minimizer != prev_minimizer) { - break; - } - prev_minimizer = minimizer; } - buckets_stats.add_num_kmers_in_super_kmer(num_super_kmers_in_bucket, w); } } std::cout << "DONE" << std::endl; - std::cout<< "numbers of superkmers: " << superkmer_idx << std::endl; + //std::cout<< "numbers of superkmers: " << superkmer_idx << std::endl; - buckets_stats.print_full(); return non_mono_superkmer; } From 74dc101663be5fd13a01419ab1e5fd32728d2151 Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Mon, 22 Apr 2024 11:52:39 +0200 Subject: [PATCH 07/22] small fixes --- include/statistics.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/statistics.cpp b/include/statistics.cpp index becb672..20c0773 100644 --- a/include/statistics.cpp +++ b/include/statistics.cpp @@ -59,7 +59,7 @@ bool equal(const std::vector& input1, const std::vector dictionary::build_superkmer_bv(const std::function (std::string_view)> &get_annotation_labels) const { - uint64_t num_kmers = size(); + //uint64_t num_kmers = size(); uint64_t num_minimizers = m_minimizers.size(); uint64_t num_super_kmers = m_buckets.offsets.size(); @@ -72,7 +72,7 @@ std::vector dictionary::build_superkmer_bv(const std::function dictionary::build_superkmer_bv(const std::function labels = get_annotation_labels(kmer_str); if(prev_labels.size() == 0){ prev_labels = labels; - } else if(non_mono_superkmer[superkmer_idx] == false && !equal(labels, prev_labels)){ - non_mono_superkmer[superkmer_idx] = true; + } else if(non_mono_superkmer[super_kmer_id] == false && !equal(labels, prev_labels)){ + non_mono_superkmer[super_kmer_id] = true; prev_labels = labels; } } From 3e9f4943753f3898f527e0401a72eba891728896 Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Wed, 24 Apr 2024 15:41:12 +0200 Subject: [PATCH 08/22] changes to decrease runtime --- include/statistics.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/include/statistics.cpp b/include/statistics.cpp index 20c0773..987694e 100644 --- a/include/statistics.cpp +++ b/include/statistics.cpp @@ -46,12 +46,12 @@ void dictionary::compute_statistics() const { buckets_stats.print_full(); std::cout << "DONE" << std::endl; } -bool equal(const std::vector& input1, const std::vector& input2) { +inline bool equal(const std::vector& input1, const std::vector& input2) { if(input1.size() != input2.size()){ return false; } for(size_t i = 0; i < input1.size(); i++){ - if(input1.at(i) != input2.at(i)){ + if(input1[i] != input2[i]){ return false; } } @@ -66,15 +66,14 @@ std::vector dictionary::build_superkmer_bv(const std::function non_mono_superkmer (num_super_kmers, false); //size_t superkmer_idx = 0; - uint64_t one_pc_buckets = num_minimizers/100; - + uint64_t one_pc_buckets = std::ceil(num_minimizers/100.0); std::cout<<"iterating through buckets :\n"; for (uint64_t bucket_id = 0; bucket_id != num_minimizers; ++bucket_id) { if(bucket_id%one_pc_buckets==0)std::cout << bucket_id/one_pc_buckets <<"%" << '\r'<< std::flush; - auto [begin, end] = m_buckets.locate_bucket(bucket_id); + auto [begin, end] = m_buckets.locate_bucket(bucket_id); //uint64_t num_super_kmers_in_bucket = end - begin; for (uint64_t super_kmer_id = begin; super_kmer_id != end; ++super_kmer_id) { - uint64_t offset = m_buckets.offsets.access(super_kmer_id); + uint64_t offset = m_buckets.offsets.access(super_kmer_id); auto [_, contig_end] = m_buckets.offset_to_id(offset, m_k); (void)_; bit_vector_iterator bv_it(m_buckets.strings, 2 * offset); @@ -105,11 +104,13 @@ std::vector dictionary::build_superkmer_bv(const std::function labels = get_annotation_labels(kmer_str); + if(prev_labels.size() == 0){ prev_labels = labels; } else if(non_mono_superkmer[super_kmer_id] == false && !equal(labels, prev_labels)){ - non_mono_superkmer[super_kmer_id] = true; + non_mono_superkmer[super_kmer_id] = true; prev_labels = labels; + break; } } } From 1254a7d61fac45b684dc118cd3284528db2f13a5 Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Wed, 24 Apr 2024 16:04:45 +0200 Subject: [PATCH 09/22] removed unnecessary line --- include/statistics.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/include/statistics.cpp b/include/statistics.cpp index 987694e..a629d45 100644 --- a/include/statistics.cpp +++ b/include/statistics.cpp @@ -107,9 +107,8 @@ std::vector dictionary::build_superkmer_bv(const std::function Date: Fri, 26 Apr 2024 17:08:34 +0200 Subject: [PATCH 10/22] getting superkmer label batches --- include/dictionary.hpp | 2 +- include/statistics.cpp | 49 +++++++++++++++++++++++++----------------- 2 files changed, 30 insertions(+), 21 deletions(-) diff --git a/include/dictionary.hpp b/include/dictionary.hpp index 3a79159..de31540 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -95,7 +95,7 @@ struct dictionary { void superkmer_statistics() const; - std::vector build_superkmer_bv(const std::function (std::string_view)> &get_annotation_labels) const; + std::vector build_superkmer_bv(const std::function &monochromatic_labels) const; template void visit(Visitor& visitor) { diff --git a/include/statistics.cpp b/include/statistics.cpp index a629d45..f561764 100644 --- a/include/statistics.cpp +++ b/include/statistics.cpp @@ -58,7 +58,12 @@ inline bool equal(const std::vector& input1, const std::vector dictionary::build_superkmer_bv(const std::function (std::string_view)> &get_annotation_labels) const { +void uint_kmer_to_last_char(kmer_t x, char* str, uint64_t k){ + x >>= 2*(k-1); + str[0] = util::uint64_to_char(x & 3); +} + +std::vector dictionary::build_superkmer_bv(const std::function &monochromatic_labels) const { //uint64_t num_kmers = size(); uint64_t num_minimizers = m_minimizers.size(); uint64_t num_super_kmers = m_buckets.offsets.size(); @@ -70,48 +75,52 @@ std::vector dictionary::build_superkmer_bv(const std::function(m_k - m_m + 1, contig_end - offset - m_k + 1); uint64_t prev_minimizer = constants::invalid_uint64; std::vector prev_labels = {}; - uint64_t w = 0; - for (; w != window_size; ++w) { + std::string superkmer = ""; + size_t num_chars_to_add = m_k; + for (uint64_t w = 0; w != window_size; ++w) { kmer_t kmer = bv_it.read_and_advance_by_two(2 * m_k); auto [minimizer, pos] = util::compute_minimizer_pos(kmer, m_k, m_m, m_seed); if (m_canonical_parsing) { uint64_t kmer_rc = util::compute_reverse_complement(kmer, m_k); - auto [minimizer_rc, pos_rc] = - util::compute_minimizer_pos(kmer_rc, m_k, m_m, m_seed); + auto [minimizer_rc, pos_rc] = util::compute_minimizer_pos(kmer_rc, m_k, m_m, m_seed); if (minimizer_rc < minimizer) { minimizer = minimizer_rc; pos = pos_rc; } } - //check if superkmer has ended + //check if superkmer has ended if (prev_minimizer != constants::invalid_uint64 and minimizer != prev_minimizer) { break; } - prev_minimizer = minimizer; - - // get labels and compare them to previous ones - std::string kmer_str(m_k,'_'); - util::uint_kmer_to_string(kmer, &kmer_str[0], m_k); - std::vector labels = get_annotation_labels(kmer_str); - - if(prev_labels.size() == 0){ - prev_labels = labels; - } else if(!equal(labels, prev_labels)){ - non_mono_superkmer[super_kmer_id] = true; - break; + prev_minimizer = minimizer; + + //get kmers into superkmer + std::string chars_to_add(num_chars_to_add, '_'); + if(num_chars_to_add > 1){ + util::uint_kmer_to_string(kmer, &chars_to_add[0], m_k); + }else{ + uint_kmer_to_last_char(kmer, &chars_to_add[0], m_k); + num_chars_to_add = 1; } + superkmer += chars_to_add; + + } + // get labels and compare them to previous ones + if(!monochromatic_labels(superkmer)){ + non_mono_superkmer[super_kmer_id] = true; } + } } std::cout << "DONE" << std::endl; From d650c0dff4c743d8cc93c8041c65a0917ce6b59f Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Sat, 27 Apr 2024 08:42:40 +0200 Subject: [PATCH 11/22] tiny bug fix --- include/statistics.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/include/statistics.cpp b/include/statistics.cpp index f561764..eb16f8c 100644 --- a/include/statistics.cpp +++ b/include/statistics.cpp @@ -78,7 +78,7 @@ std::vector dictionary::build_superkmer_bv(const std::function dictionary::build_superkmer_bv(const std::function 1){ util::uint_kmer_to_string(kmer, &chars_to_add[0], m_k); + num_chars_to_add = 1; }else{ uint_kmer_to_last_char(kmer, &chars_to_add[0], m_k); - num_chars_to_add = 1; } superkmer += chars_to_add; - } + } + //std::cout<<"superkmer "< Date: Mon, 29 Apr 2024 00:08:59 +0200 Subject: [PATCH 12/22] paralellized bit vector construction --- include/dictionary.hpp | 1 + include/statistics.cpp | 18 +++++++----------- 2 files changed, 8 insertions(+), 11 deletions(-) diff --git a/include/dictionary.hpp b/include/dictionary.hpp index de31540..0d63ca4 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -5,6 +5,7 @@ #include "buckets.hpp" #include "skew_index.hpp" #include "weights.hpp" +#include namespace sshash { diff --git a/include/statistics.cpp b/include/statistics.cpp index eb16f8c..8b46977 100644 --- a/include/statistics.cpp +++ b/include/statistics.cpp @@ -73,6 +73,10 @@ std::vector dictionary::build_superkmer_bv(const std::function dictionary::build_superkmer_bv(const std::function dictionary::build_superkmer_bv(const std::function Date: Mon, 29 Apr 2024 16:32:47 +0200 Subject: [PATCH 13/22] splitting and merging of superkmer vector for parallelization --- include/dictionary.hpp | 2 +- include/statistics.cpp | 43 +++++++++++++++++++++++++++++++----------- 2 files changed, 33 insertions(+), 12 deletions(-) diff --git a/include/dictionary.hpp b/include/dictionary.hpp index 0d63ca4..9d69cc7 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -96,7 +96,7 @@ struct dictionary { void superkmer_statistics() const; - std::vector build_superkmer_bv(const std::function &monochromatic_labels) const; + std::vector build_superkmer_bv(const std::function &monochromatic_labels) const; template void visit(Visitor& visitor) { diff --git a/include/statistics.cpp b/include/statistics.cpp index 8b46977..10db83b 100644 --- a/include/statistics.cpp +++ b/include/statistics.cpp @@ -1,5 +1,9 @@ +#include +#include +#include #include "dictionary.hpp" #include "buckets_statistics.hpp" +#include namespace sshash { void dictionary::compute_statistics() const { @@ -63,22 +67,28 @@ void uint_kmer_to_last_char(kmer_t x, char* str, uint64_t k){ str[0] = util::uint64_to_char(x & 3); } -std::vector dictionary::build_superkmer_bv(const std::function &monochromatic_labels) const { +std::vector dictionary::build_superkmer_bv(const std::function &monochromatic_labels) const { //uint64_t num_kmers = size(); uint64_t num_minimizers = m_minimizers.size(); - uint64_t num_super_kmers = m_buckets.offsets.size(); + //uint64_t num_super_kmers = m_buckets.offsets.size(); std::cout << "building super kmer mask..." << std::endl; - std::vector non_mono_superkmer (num_super_kmers, false); + //std::vector non_mono_superkmer (num_super_kmers, false); //size_t superkmer_idx = 0; uint64_t one_pc_buckets = std::ceil(num_minimizers/100.0); std::cout<<"iterating through buckets :\n"; uint64_t num_threads = std::thread::hardware_concurrency(); if (num_minimizers < num_threads) num_threads = num_minimizers; - std::mutex vec_mutex; + std::vector> indeces (num_threads); + time_t my_time = time(NULL); + printf("%s", ctime(&my_time)); #pragma omp parallel for num_threads(num_threads) schedule(dynamic) for (uint64_t bucket_id = 0; bucket_id != num_minimizers; ++bucket_id) { - if(bucket_id%one_pc_buckets==0)std::cout << bucket_id/one_pc_buckets <<"%" << '\r'<< std::flush; + if(bucket_id%one_pc_buckets==0){ + my_time = time(NULL); + printf("%s", ctime(&my_time)); + std::cout <<" "<< bucket_id/one_pc_buckets <<"%" << std::endl; + } auto [begin, end] = m_buckets.locate_bucket(bucket_id); //uint64_t num_super_kmers_in_bucket = end - begin; for (uint64_t super_kmer_id = begin; super_kmer_id != end; ++super_kmer_id) { @@ -114,16 +124,27 @@ std::vector dictionary::build_superkmer_bv(const std::function & first, const std::vector& second) {return first.size() < second.size(); }); + std::vector merged = std::move(indeces[0]); + for(uint64_t i = 1; i < num_threads; i++){ + std::vector dest_aux; + auto& src = indeces[i]; + std::merge(merged.begin(),merged.end(),src.begin(), src.end(),std::back_inserter(dest_aux)); + merged = std::move(dest_aux); + } + my_time = time(NULL); + printf("%s", ctime(&my_time)); + std::cout<< " merged! \n"; + return merged; } } // namespace sshash From f2eaf607a827dfeea3fcc32df0337df4206297f6 Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Thu, 2 May 2024 00:52:29 +0200 Subject: [PATCH 14/22] more efficient query for kmers in non-monochromatic superkmers --- include/dictionary.cpp | 4 ++++ include/dictionary.hpp | 2 ++ 2 files changed, 6 insertions(+) diff --git a/include/dictionary.cpp b/include/dictionary.cpp index 5a8bbe2..39989f0 100644 --- a/include/dictionary.cpp +++ b/include/dictionary.cpp @@ -28,6 +28,10 @@ std::pair dictionary::kmer_to_superkmer_idx_helper(kmer return m_buckets.lookup_superkmer_start(begin, end, uint_kmer, m_k, m_m); } +uint64_t dictionary::look_up_from_superkmer_id(uint64_t superkmer_id, char const* kmer_str){ + kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); + return m_buckets.lookup_in_super_kmer(superkmer_id, uint_kmer, m_k, m_m).kmer_id; +} lookup_result dictionary::lookup_uint_regular_parsing(kmer_t uint_kmer) const { uint64_t minimizer = util::compute_minimizer(uint_kmer, m_k, m_m, m_seed); uint64_t bucket_id = m_minimizers.lookup(minimizer); diff --git a/include/dictionary.hpp b/include/dictionary.hpp index 3a79159..da56713 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -111,6 +111,8 @@ struct dictionary { } std::pair kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const ; std::pair kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const ; + uint64_t look_up_from_superkmer_id(uint64_t superkmer_id, char const* kmer_str); + private: uint64_t m_size; From 9e48f313701583866ad44f55968de5cb88e5d7d9 Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Tue, 7 May 2024 14:15:28 +0200 Subject: [PATCH 15/22] reverse complement bug fix --- include/dictionary.cpp | 12 ++++++++---- include/dictionary.hpp | 2 +- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/include/dictionary.cpp b/include/dictionary.cpp index 39989f0..1fe28b9 100644 --- a/include/dictionary.cpp +++ b/include/dictionary.cpp @@ -28,9 +28,15 @@ std::pair dictionary::kmer_to_superkmer_idx_helper(kmer return m_buckets.lookup_superkmer_start(begin, end, uint_kmer, m_k, m_m); } -uint64_t dictionary::look_up_from_superkmer_id(uint64_t superkmer_id, char const* kmer_str){ +uint64_t dictionary::look_up_from_superkmer_id(uint64_t superkmer_id, char const* kmer_str, bool check_reverse_complement){ kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); - return m_buckets.lookup_in_super_kmer(superkmer_id, uint_kmer, m_k, m_m).kmer_id; + + auto res = m_buckets.lookup_in_super_kmer(superkmer_id, uint_kmer, m_k, m_m).kmer_id; + if (check_reverse_complement and res == constants::invalid_uint64) { + kmer_t uint_kmer_rc = util::compute_reverse_complement(uint_kmer, m_k); + res = m_buckets.lookup_in_super_kmer(superkmer_id, uint_kmer_rc, m_k, m_m).kmer_id; + } + return res; } lookup_result dictionary::lookup_uint_regular_parsing(kmer_t uint_kmer) const { uint64_t minimizer = util::compute_minimizer(uint_kmer, m_k, m_m, m_seed); @@ -110,13 +116,11 @@ std::pair dictionary::kmer_to_superkmer_idx(char const* kmer kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); auto[res, s_idx] = kmer_to_superkmer_idx_helper(uint_kmer); - assert(res.kmer_orientation == constants::forward_orientation); if(res.kmer_id == constants::invalid_uint64 && check_reverse_complement){ kmer_t uint_kmer_rc = util::compute_reverse_complement(uint_kmer, m_k); auto pair = kmer_to_superkmer_idx_helper(uint_kmer_rc); res = pair.first; s_idx = pair.second; - res.kmer_orientation = constants::backward_orientation; } return std::pair(res.kmer_id, s_idx); } diff --git a/include/dictionary.hpp b/include/dictionary.hpp index 9335c88..2a21367 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -112,7 +112,7 @@ struct dictionary { } std::pair kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const ; std::pair kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const ; - uint64_t look_up_from_superkmer_id(uint64_t superkmer_id, char const* kmer_str); + uint64_t look_up_from_superkmer_id(uint64_t superkmer_id, char const* kmer_str, bool check_reverse_complement); private: From 9be36ab842cebf07000b4a8b6a7fa78cb131bb37 Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Sun, 12 May 2024 13:36:08 +0200 Subject: [PATCH 16/22] returning kmer index during superkmer lookup --- include/buckets.hpp | 12 ++++++++---- include/dictionary.cpp | 29 +++++++++++------------------ include/dictionary.hpp | 7 ++++--- 3 files changed, 23 insertions(+), 25 deletions(-) diff --git a/include/buckets.hpp b/include/buckets.hpp index 3359e21..dbe1c5b 100644 --- a/include/buckets.hpp +++ b/include/buckets.hpp @@ -5,7 +5,11 @@ #include "ef_sequence.hpp" namespace sshash { - +struct superkmer_result { + uint64_t kmer_idx; + uint64_t superkmer_idx; + uint64_t superkmer_id; + }; struct buckets { std::pair offset_to_id(uint64_t offset, uint64_t k) const { auto [pos, contig_begin, contig_end] = pieces.locate(offset); @@ -111,16 +115,16 @@ struct buckets { } - std::pair lookup_superkmer_start(uint64_t begin, uint64_t end, kmer_t target_kmer, + superkmer_result lookup_superkmer_start(uint64_t begin, uint64_t end, kmer_t target_kmer, uint64_t k, uint64_t m) const { for (uint64_t super_kmer_id = begin; super_kmer_id != end; ++super_kmer_id) { auto res = lookup_in_super_kmer(super_kmer_id, target_kmer, k, m); if(res.kmer_id != constants::invalid_uint64){ assert(is_valid(res)); - return std::pair(superkmer_id_to_kmer_id(super_kmer_id, k), super_kmer_id); // CHECK THIS IS CORRECT!?!?! YES + return {res.kmer_id, superkmer_id_to_kmer_id(super_kmer_id, k).kmer_id, super_kmer_id}; } } - return std::pair(lookup_result(), constants::invalid_uint64); + return {constants::invalid_uint64, constants::invalid_uint64, constants::invalid_uint64}; } diff --git a/include/dictionary.cpp b/include/dictionary.cpp index 1fe28b9..28ee623 100644 --- a/include/dictionary.cpp +++ b/include/dictionary.cpp @@ -3,7 +3,7 @@ namespace sshash { ////////////////////////////////////////// -std::pair dictionary::kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const { +superkmer_result dictionary::kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const { uint64_t minimizer = util::compute_minimizer(uint_kmer, m_k, m_m, m_seed); uint64_t bucket_id = m_minimizers.lookup(minimizer); auto [begin, end] = m_buckets.locate_bucket(bucket_id); @@ -17,12 +17,13 @@ std::pair dictionary::kmer_to_superkmer_idx_helper(kmer uint64_t pos = m_skew_index.lookup(uint_kmer, log2_bucket_size); /* It must hold pos < num_super_kmers_in_bucket for the kmer to exist. */ if (pos < num_super_kmers_in_bucket) { - if(m_buckets.lookup_in_super_kmer(begin + pos, uint_kmer, m_k, m_m).kmer_id != sshash::constants::invalid_uint64){ + auto res_kmer = m_buckets.lookup_in_super_kmer(begin + pos, uint_kmer, m_k, m_m).kmer_id; + if(res_kmer != sshash::constants::invalid_uint64){ uint64_t superkmer_id = begin + pos; - return std::pair(m_buckets.superkmer_id_to_kmer_id(superkmer_id, m_k), superkmer_id); // CHECK THIS IS CORRECT!?!?! YES + return {res_kmer, m_buckets.superkmer_id_to_kmer_id(superkmer_id, m_k).kmer_id, superkmer_id}; } } - return std::pair(lookup_result(), constants::invalid_uint64); + return {constants::invalid_uint64, constants::invalid_uint64, constants::invalid_uint64}; } // superkmer not in skew index return m_buckets.lookup_superkmer_start(begin, end, uint_kmer, m_k, m_m); @@ -106,23 +107,15 @@ lookup_result dictionary::lookup_advanced(char const* string_kmer, return lookup_advanced_uint(uint_kmer, check_reverse_complement); } ////////////////////////////////////////// -std::pair dictionary::kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const { - // Plan: - // kmer -> get minimizer -> get bucket ->skew_index? -> - // ->get superkmer range in offsets!!-> - // ->get superkmer coords in offsets - // (+check everything makes sense on contig level) - // reverse complement considerations: primary graph -> look for both minimizers of kmer and rc_kmer - +superkmer_result dictionary::kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const { + kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); - auto[res, s_idx] = kmer_to_superkmer_idx_helper(uint_kmer); - if(res.kmer_id == constants::invalid_uint64 && check_reverse_complement){ + superkmer_result sk_res = kmer_to_superkmer_idx_helper(uint_kmer); + if(sk_res.kmer_idx == constants::invalid_uint64 && check_reverse_complement){ kmer_t uint_kmer_rc = util::compute_reverse_complement(uint_kmer, m_k); - auto pair = kmer_to_superkmer_idx_helper(uint_kmer_rc); - res = pair.first; - s_idx = pair.second; + sk_res = kmer_to_superkmer_idx_helper(uint_kmer_rc); } - return std::pair(res.kmer_id, s_idx); + return sk_res; } lookup_result dictionary::lookup_advanced_uint(kmer_t uint_kmer, bool check_reverse_complement) const { diff --git a/include/dictionary.hpp b/include/dictionary.hpp index 2a21367..8c57122 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -110,10 +110,11 @@ struct dictionary { visitor.visit(m_skew_index); visitor.visit(m_weights); } - std::pair kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const ; - std::pair kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const ; + superkmer_result kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const ; + superkmer_result kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const ; uint64_t look_up_from_superkmer_id(uint64_t superkmer_id, char const* kmer_str, bool check_reverse_complement); - + + private: uint64_t m_size; From dd9684fd8d9d9c8d95bf991ce82a48edd82cf109 Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Fri, 24 May 2024 11:30:52 +0200 Subject: [PATCH 17/22] added separate stats function --- include/dictionary.hpp | 5 +++ include/statistics.cpp | 89 ++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 90 insertions(+), 4 deletions(-) diff --git a/include/dictionary.hpp b/include/dictionary.hpp index 8c57122..48c4647 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -96,7 +96,12 @@ struct dictionary { void superkmer_statistics() const; + // pre: monochromatic_labels must return true iff same labels are found for all k-mers in given string + // post: returns indeces of superkmers that are not monochromatic std::vector build_superkmer_bv(const std::function &monochromatic_labels) const; + //pre: get_annotation_labels must return vector of labels for given k-mer + //post: number of k-mers and color changes for every super-k-mer are written to files + void make_superkmer_stats(const std::function (std::string_view)> &get_annotation_labels) const; template void visit(Visitor& visitor) { diff --git a/include/statistics.cpp b/include/statistics.cpp index 10db83b..224dd92 100644 --- a/include/statistics.cpp +++ b/include/statistics.cpp @@ -115,7 +115,7 @@ std::vector dictionary::build_superkmer_bv(const std::function 1){ util::uint_kmer_to_string(kmer, &chars_to_add[0], m_k); - num_chars_to_add = 1; + num_chars_to_add = 1; }else{ uint_kmer_to_last_char(kmer, &chars_to_add[0], m_k); } @@ -137,14 +137,95 @@ std::vector dictionary::build_superkmer_bv(const std::function merged = std::move(indeces[0]); for(uint64_t i = 1; i < num_threads; i++){ std::vector dest_aux; - auto& src = indeces[i]; - std::merge(merged.begin(),merged.end(),src.begin(), src.end(),std::back_inserter(dest_aux)); - merged = std::move(dest_aux); + auto& src = indeces[i]; + std::merge(merged.begin(),merged.end(),src.begin(), src.end(),std::back_inserter(dest_aux)); + merged = std::move(dest_aux); } my_time = time(NULL); printf("%s", ctime(&my_time)); std::cout<< " merged! \n"; return merged; } +void write_vec(const std::vector& vec, std::string output_path){ + std::ofstream outputFile(output_path); + if (!outputFile.is_open()) { + std::cerr << "Error opening the file." << std::endl; + } + for (const auto& element : vec) { + outputFile << element << " "; + } + outputFile.close(); +} +void dictionary::make_superkmer_stats(const std::function (std::string_view)> &get_annotation_labels) const { + uint64_t num_minimizers = m_minimizers.size(); + uint64_t num_super_kmers = m_buckets.offsets.size(); + + std::cout << "gathering super-k-mer stats..." << std::endl; + std::vector kmer_per_superkmer (num_super_kmers, 0); + std::vector color_changes_superkmer (num_super_kmers, 0); + + uint64_t num_threads = std::thread::hardware_concurrency(); + if (num_minimizers < num_threads) num_threads = num_minimizers; + std::mutex color_lock; + std::mutex kmer_lock; + + uint64_t one_pc_buckets = std::ceil(num_minimizers/100.0); + time_t my_time = time(NULL); + std::cout<<"iterating through buckets :\n"; + #pragma omp parallel for num_threads(num_threads) schedule(dynamic) + for (uint64_t bucket_id = 0; bucket_id != num_minimizers; ++bucket_id) { + if(bucket_id%one_pc_buckets==0){ + my_time = time(NULL); + printf("%s", ctime(&my_time)); + std::cout <<" "<< bucket_id/one_pc_buckets <<"%" << std::endl; + } + auto [begin, end] = m_buckets.locate_bucket(bucket_id); + for (uint64_t super_kmer_id = begin; super_kmer_id != end; ++super_kmer_id) { + uint64_t offset = m_buckets.offsets.access(super_kmer_id); + auto [_, contig_end] = m_buckets.offset_to_id(offset, m_k); + (void)_; + bit_vector_iterator bv_it(m_buckets.strings, 2 * offset); + uint64_t window_size = std::min(m_k - m_m + 1, contig_end - offset - m_k + 1); + uint64_t prev_minimizer = constants::invalid_uint64; + std::vector prev_labels = {}; + + unsigned num_kmers = 0; + unsigned num_col = 0; + for (uint64_t w = 0; w != window_size; ++w) { + kmer_t kmer = bv_it.read_and_advance_by_two(2 * m_k); + auto [minimizer, pos] = util::compute_minimizer_pos(kmer, m_k, m_m, m_seed); + + //check if superkmer has ended + if (prev_minimizer != constants::invalid_uint64 and minimizer != prev_minimizer) { + break; + } + prev_minimizer = minimizer; + num_kmers++; + + // get labels and compare them to previous ones + std::string kmer_str(m_k,'_'); + util::uint_kmer_to_string(kmer, &kmer_str[0], m_k); + std::vector labels = get_annotation_labels(kmer_str); + + if(prev_labels.size() == 0){ + prev_labels = labels; + } else if(!equal(labels, prev_labels)){ + num_col++; + } + } + color_lock.lock(); + color_changes_superkmer[super_kmer_id] = num_col; + color_lock.unlock(); + kmer_lock.lock(); + kmer_per_superkmer[super_kmer_id] = num_kmers; + kmer_lock.unlock(); + } + } + std::cout << "DONE" << std::endl; + + std::cout<<"writing kmer and color stats to files...\n"; + write_vec(color_changes_superkmer,"final_color_stats.txt"); + write_vec(kmer_per_superkmer,"final_kmer_stats.txt"); +} } // namespace sshash From e156b27fc6f604c466b90e7f26a7c4e9279cc18d Mon Sep 17 00:00:00 2001 From: Marianna Marzetta Date: Fri, 24 May 2024 13:11:23 +0200 Subject: [PATCH 18/22] clean-up --- include/dictionary.cpp | 5 ++--- include/dictionary.hpp | 17 +++++++---------- 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/include/dictionary.cpp b/include/dictionary.cpp index 28ee623..1effe4e 100644 --- a/include/dictionary.cpp +++ b/include/dictionary.cpp @@ -2,7 +2,6 @@ namespace sshash { -////////////////////////////////////////// superkmer_result dictionary::kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const { uint64_t minimizer = util::compute_minimizer(uint_kmer, m_k, m_m, m_seed); uint64_t bucket_id = m_minimizers.lookup(minimizer); @@ -29,7 +28,7 @@ superkmer_result dictionary::kmer_to_superkmer_idx_helper(kmer_t uint_kmer) cons return m_buckets.lookup_superkmer_start(begin, end, uint_kmer, m_k, m_m); } -uint64_t dictionary::look_up_from_superkmer_id(uint64_t superkmer_id, char const* kmer_str, bool check_reverse_complement){ +uint64_t dictionary::look_up_from_superkmer_id(uint64_t superkmer_id, char const* kmer_str, bool check_reverse_complement) const { kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); auto res = m_buckets.lookup_in_super_kmer(superkmer_id, uint_kmer, m_k, m_m).kmer_id; @@ -106,7 +105,7 @@ lookup_result dictionary::lookup_advanced(char const* string_kmer, kmer_t uint_kmer = util::string_to_uint_kmer(string_kmer, m_k); return lookup_advanced_uint(uint_kmer, check_reverse_complement); } -////////////////////////////////////////// + superkmer_result dictionary::kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const { kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); diff --git a/include/dictionary.hpp b/include/dictionary.hpp index 48c4647..1454db1 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -94,15 +94,17 @@ struct dictionary { void print_space_breakdown() const; void compute_statistics() const; - void superkmer_statistics() const; - - // pre: monochromatic_labels must return true iff same labels are found for all k-mers in given string + // pre: monochromatic_labels returns true iff same labels are found for all k-mers in given string // post: returns indeces of superkmers that are not monochromatic std::vector build_superkmer_bv(const std::function &monochromatic_labels) const; - //pre: get_annotation_labels must return vector of labels for given k-mer - //post: number of k-mers and color changes for every super-k-mer are written to files + // pre: get_annotation_labels returns vector of labels for given k-mer + // post: number of k-mers and color changes for every super-k-mer are written to files void make_superkmer_stats(const std::function (std::string_view)> &get_annotation_labels) const; + superkmer_result kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const ; + superkmer_result kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const ; + uint64_t look_up_from_superkmer_id(uint64_t superkmer_id, char const* kmer_str, bool check_reverse_complement) const ; + template void visit(Visitor& visitor) { visitor.visit(m_size); @@ -115,11 +117,6 @@ struct dictionary { visitor.visit(m_skew_index); visitor.visit(m_weights); } - superkmer_result kmer_to_superkmer_idx_helper(kmer_t uint_kmer) const ; - superkmer_result kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const ; - uint64_t look_up_from_superkmer_id(uint64_t superkmer_id, char const* kmer_str, bool check_reverse_complement); - - private: uint64_t m_size; From aa66464c9dd9f4ea681548fd04f4ce2eaf5c457a Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Thu, 25 Jul 2024 15:56:59 +0200 Subject: [PATCH 19/22] Fix compilation errors --- include/dictionary.hpp | 2 ++ include/dictionary.impl | 10 ++++++---- include/statistics.impl | 13 ++++++++----- 3 files changed, 16 insertions(+), 9 deletions(-) diff --git a/include/dictionary.hpp b/include/dictionary.hpp index 97437e8..e1f5863 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -1,5 +1,7 @@ #pragma once +#include + #include "util.hpp" #include "minimizers.hpp" #include "buckets.hpp" diff --git a/include/dictionary.impl b/include/dictionary.impl index e6ae767..e3ddabf 100644 --- a/include/dictionary.impl +++ b/include/dictionary.impl @@ -30,11 +30,12 @@ superkmer_result dictionary::kmer_to_superkmer_idx_helper(kmer_t uint_km template uint64_t dictionary::look_up_from_superkmer_id(uint64_t superkmer_id, char const* kmer_str, bool check_reverse_complement) const { - kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); + kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); auto res = m_buckets.lookup_in_super_kmer(superkmer_id, uint_kmer, m_k, m_m).kmer_id; if (check_reverse_complement and res == constants::invalid_uint64) { - kmer_t uint_kmer_rc = util::compute_reverse_complement(uint_kmer, m_k); + kmer_t uint_kmer_rc = uint_kmer; + uint_kmer_rc.reverse_complement_inplace(m_k); res = m_buckets.lookup_in_super_kmer(superkmer_id, uint_kmer_rc, m_k, m_m).kmer_id; } return res; @@ -110,10 +111,11 @@ uint64_t dictionary::lookup_uint(kmer_t uint_kmer, bool check_reverse_co template superkmer_result dictionary::kmer_to_superkmer_idx(char const* kmer_str, bool check_reverse_complement) const { - kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); + kmer_t uint_kmer = util::string_to_uint_kmer(kmer_str, m_k); superkmer_result sk_res = kmer_to_superkmer_idx_helper(uint_kmer); if(sk_res.kmer_idx == constants::invalid_uint64 && check_reverse_complement){ - kmer_t uint_kmer_rc = util::compute_reverse_complement(uint_kmer, m_k); + kmer_t uint_kmer_rc = uint_kmer; + uint_kmer_rc.reverse_complement_inplace(m_k); sk_res = kmer_to_superkmer_idx_helper(uint_kmer_rc); } return sk_res; diff --git a/include/statistics.impl b/include/statistics.impl index c9f6824..ae33275 100644 --- a/include/statistics.impl +++ b/include/statistics.impl @@ -64,12 +64,14 @@ inline bool equal(const std::vector& input1, const std::vector>= 2*(k-1); - str[0] = util::uint64_to_char(x & 3); +template +void uint_kmer_to_last_char(kmer_t x, char* str, uint64_t k) { + x.drop_chars(k - 1); + str[0] = kmer_t::uint64_to_char(x.pop_char(1)); } -std::vector dictionary::build_superkmer_bv(const std::function &monochromatic_labels) const { +template +std::vector dictionary::build_superkmer_bv(const std::function &monochromatic_labels) const { //uint64_t num_kmers = size(); uint64_t num_minimizers = m_minimizers.size(); //uint64_t num_super_kmers = m_buckets.offsets.size(); @@ -158,7 +160,8 @@ void write_vec(const std::vector& vec, std::string output_path){ } outputFile.close(); } -void dictionary::make_superkmer_stats(const std::function (std::string_view)> &get_annotation_labels) const { +template +void dictionary::make_superkmer_stats(const std::function (std::string_view)> &get_annotation_labels) const { uint64_t num_minimizers = m_minimizers.size(); uint64_t num_super_kmers = m_buckets.offsets.size(); From 25cf8e443d8db6b0b8178c16bbf5c2fe41f951a6 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Thu, 25 Jul 2024 16:36:17 +0200 Subject: [PATCH 20/22] Fixes --- include/statistics.impl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/include/statistics.impl b/include/statistics.impl index ae33275..b36ddb3 100644 --- a/include/statistics.impl +++ b/include/statistics.impl @@ -67,7 +67,7 @@ inline bool equal(const std::vector& input1, const std::vector void uint_kmer_to_last_char(kmer_t x, char* str, uint64_t k) { x.drop_chars(k - 1); - str[0] = kmer_t::uint64_to_char(x.pop_char(1)); + str[0] = kmer_t::uint64_to_char(x.pop_char()); } template @@ -99,14 +99,14 @@ std::vector dictionary::build_superkmer_bv(const std::function uint64_t offset = m_buckets.offsets.access(super_kmer_id); auto [_, contig_end] = m_buckets.offset_to_id(offset, m_k); (void)_; - bit_vector_iterator bv_it(m_buckets.strings, 2 * offset); + bit_vector_iterator bv_it(m_buckets.strings, kmer_t::bits_per_char * offset); uint64_t window_size = std::min(m_k - m_m + 1, contig_end - offset - m_k + 1); - uint64_t prev_minimizer = constants::invalid_uint64; + kmer_t prev_minimizer = constants::invalid_uint64; std::vector prev_labels = {}; std::string superkmer = ""; size_t num_chars_to_add = m_k; for (uint64_t w = 0; w != window_size; ++w) { - kmer_t kmer = bv_it.read_and_advance_by_two(2 * m_k); + kmer_t kmer = bv_it.read_and_advance_by_char(m_k); auto [minimizer, pos] = util::compute_minimizer_pos(kmer, m_k, m_m, m_seed); //check if superkmer has ended @@ -150,7 +150,7 @@ std::vector dictionary::build_superkmer_bv(const std::function std::cout<< " merged! \n"; return merged; } -void write_vec(const std::vector& vec, std::string output_path){ +static void write_vec(const std::vector& vec, std::string output_path){ std::ofstream outputFile(output_path); if (!outputFile.is_open()) { std::cerr << "Error opening the file." << std::endl; @@ -189,15 +189,15 @@ void dictionary::make_superkmer_stats(const std::function bv_it(m_buckets.strings, kmer_t::bits_per_char * offset); uint64_t window_size = std::min(m_k - m_m + 1, contig_end - offset - m_k + 1); - uint64_t prev_minimizer = constants::invalid_uint64; + kmer_t prev_minimizer = constants::invalid_uint64; std::vector prev_labels = {}; unsigned num_kmers = 0; unsigned num_col = 0; for (uint64_t w = 0; w != window_size; ++w) { - kmer_t kmer = bv_it.read_and_advance_by_two(2 * m_k); + kmer_t kmer = bv_it.read_and_advance_by_char(m_k); auto [minimizer, pos] = util::compute_minimizer_pos(kmer, m_k, m_m, m_seed); //check if superkmer has ended From 4246537ae422c332882b2449c1933d2d966df9ef Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Thu, 25 Jul 2024 17:30:14 +0200 Subject: [PATCH 21/22] Fix --- include/statistics.impl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/statistics.impl b/include/statistics.impl index b36ddb3..86462f5 100644 --- a/include/statistics.impl +++ b/include/statistics.impl @@ -15,8 +15,8 @@ void dictionary::compute_statistics() const { buckets_statistics buckets_stats(num_minimizers, num_kmers, num_super_kmers); std::cout << "computing buckets statistics..." << std::endl; - - for (uint64_t bucket_id = 0; bucket_id != num_minimizers; ++bucket_id) { + #pragma omp parallel for num_threads(num_threads) schedule(dynamic) + for (uint64_t bucket_id = 0; bucket_id < num_minimizers; ++bucket_id) { auto [begin, end] = m_buckets.locate_bucket(bucket_id); uint64_t num_super_kmers_in_bucket = end - begin; buckets_stats.add_num_super_kmers_in_bucket(num_super_kmers_in_bucket); @@ -87,7 +87,7 @@ std::vector dictionary::build_superkmer_bv(const std::function time_t my_time = time(NULL); printf("%s", ctime(&my_time)); #pragma omp parallel for num_threads(num_threads) schedule(dynamic) - for (uint64_t bucket_id = 0; bucket_id != num_minimizers; ++bucket_id) { + for (uint64_t bucket_id = 0; bucket_id < num_minimizers; ++bucket_id) { if(bucket_id%one_pc_buckets==0){ my_time = time(NULL); printf("%s", ctime(&my_time)); @@ -178,7 +178,7 @@ void dictionary::make_superkmer_stats(const std::function Date: Thu, 25 Jul 2024 18:58:16 +0200 Subject: [PATCH 22/22] Fix --- include/statistics.impl | 1 - 1 file changed, 1 deletion(-) diff --git a/include/statistics.impl b/include/statistics.impl index 86462f5..e5b7518 100644 --- a/include/statistics.impl +++ b/include/statistics.impl @@ -15,7 +15,6 @@ void dictionary::compute_statistics() const { buckets_statistics buckets_stats(num_minimizers, num_kmers, num_super_kmers); std::cout << "computing buckets statistics..." << std::endl; - #pragma omp parallel for num_threads(num_threads) schedule(dynamic) for (uint64_t bucket_id = 0; bucket_id < num_minimizers; ++bucket_id) { auto [begin, end] = m_buckets.locate_bucket(bucket_id); uint64_t num_super_kmers_in_bucket = end - begin;