diff --git a/metagraph/integration_tests/test_transform_anno.py b/metagraph/integration_tests/test_transform_anno.py index 1c1be5493e..3cb530f6f8 100644 --- a/metagraph/integration_tests/test_transform_anno.py +++ b/metagraph/integration_tests/test_transform_anno.py @@ -7,7 +7,7 @@ import glob import os import gzip -from base import TestingBase, METAGRAPH, TEST_DATA_DIR, graph_file_extension, anno_file_extension +from base import TestingBase, METAGRAPH, TEST_DATA_DIR, graph_file_extension, anno_file_extension, DNA_MODE from helpers import get_test_class_name """Test operations on annotation columns""" @@ -546,5 +546,71 @@ def test_round_trip_column_to_brwt_to_flat_to_disk(self): self._assert_query_identical(column_path, converted) +@parameterized_class(('graph_mode',), + # protein graphs have no reverse complement, so only basic mode + input_values=[('basic',)] + ([('canonical',), ('primary',)] if DNA_MODE else []), + class_name_func=get_test_class_name +) +class TestRowDiffCharSSHash(TestingBase): + """row_diff_brwt -> row_diff_char_brwt on an SSHash graph must not change + any query result: only the successor encoding differs, the diffs are the same.""" + + @classmethod + def setUpClass(cls): + cls.tempdir = TemporaryDirectory() + cls.graph_path = cls.tempdir.name + '/graph.sshashdbg' + # SSHash is built from contigs without repeated k-mers, as in the unit tests + # (building it from the raw reads makes it write minimizer files without bound) + contig_mode = 'canonical' if cls.graph_mode == 'primary' else cls.graph_mode + cls._run_command(f'{METAGRAPH} build -p {NUM_THREADS} --mode {contig_mode} -k 21 ' + f'-o {cls.tempdir.name}/contigs {TEST_DATA_DIR}/transcripts_100.fa', + "Build command") + primary_flag = ' --primary-kmers' if cls.graph_mode == 'primary' else '' + cls._run_command(f'{METAGRAPH} transform -p {NUM_THREADS} --to-fasta{primary_flag} ' + f'-o {cls.tempdir.name}/contigs {cls.tempdir.name}/contigs.dbg', + "Transform command") + input_fa = cls.tempdir.name + '/contigs.fasta.gz' + build_mode = cls.graph_mode + cls._run_command(f'{METAGRAPH} build -p {NUM_THREADS} --graph sshash --mode {build_mode} ' + f'-k 21 -o {cls.graph_path} {input_fa}', "Build command") + + cls._annotate_graph(TEST_DATA_DIR + '/transcripts_100.fa', + cls.graph_path, + cls.tempdir.name + '/annotation', + 'row_diff_brwt') + cls.brwt_path = cls.tempdir.name + '/annotation.row_diff_brwt.annodbg' + + out_base = cls.tempdir.name + '/annotation_char' + command = (f'{METAGRAPH} transform_anno --anno-type row_diff_char_brwt ' + f'-i {cls.graph_path} -o {out_base} -p {NUM_THREADS} {cls.brwt_path}') + cls._run_command(command, "Transform_anno to row_diff_char_brwt") + cls.char_path = out_base + '.row_diff_char_brwt.annodbg' + + def test_stats_match(self): + src = self._get_stats(f'-a {self.brwt_path}') + dst = self._get_stats(f'-a {self.char_path}') + self.assertEqual(src['returncode'], 0) + self.assertEqual(dst['returncode'], 0) + self.assertEqual('100', dst['labels']) + self.assertEqual(src['labels'], dst['labels']) + self.assertEqual(src['objects'], dst['objects']) + self.assertEqual(src['density'], dst['density']) + self.assertEqual('row_diff_char_brwt', dst['representation']) + + @parameterized.expand(itertools.product(['labels', 'matches', 'signature'], [1, NUM_THREADS])) + def test_query_identical(self, query_mode, num_threads): + def query(anno_path): + cmd = (f'{METAGRAPH} query -i {self.graph_path} -a {anno_path} ' + f'--query-mode {query_mode} -p {num_threads} ' + f'--min-kmers-fraction-label 0 {TEST_DATA_DIR}/transcripts_1000.fa') + res = subprocess.run(cmd.split(), stdout=PIPE, stderr=PIPE) + self.assertEqual(res.returncode, 0, msg=res.stderr.decode()) + # with several threads, results are printed in no fixed order + return sorted(res.stdout.splitlines()) + expected = query(self.brwt_path) + self.assertGreater(len(expected), 0) + self.assertEqual(expected, query(self.char_path)) + + if __name__ == '__main__': unittest.main() diff --git a/metagraph/src/annotation/annotation_converters.cpp b/metagraph/src/annotation/annotation_converters.cpp index 03aa1d73cb..cae6e09451 100644 --- a/metagraph/src/annotation/annotation_converters.cpp +++ b/metagraph/src/annotation/annotation_converters.cpp @@ -300,6 +300,84 @@ std::pair get_anchors_and_fork_fnames(const std::strin return std::make_pair(anchors_file, fork_succ_file); } +std::pair> +compute_rd_succ_char(const graph::DeBruijnGraph &graph, + const bit_vector &fork_succ) { + using node_index = graph::DeBruijnGraph::node_index; + + const auto &alphabet = graph.alphabet(); + // skip sentinel '$' if present at position 0 + size_t offset = (alphabet[0] == '$') ? 1 : 0; + size_t alpha_size = alphabet.size() - offset; + + // bits per character: ceil(log2(alpha_size)) + uint8_t bits_per_char = 0; + for (size_t s = alpha_size - 1; s > 0; s >>= 1) + ++bits_per_char; + + // build char-to-index map for the non-sentinel alphabet + std::array char_to_idx{}; + for (size_t i = 0; i < alpha_size; ++i) { + char_to_idx[static_cast(alphabet[i + offset])] = i; + } + + // Nodes are processed in blocks in parallel. The block size is a multiple of 64, + // so no two threads write to the same word of |needs_char_bv|. Each block keeps + // its chars packed until all blocks are concatenated in node order. + const uint64_t num_rows = graph.max_index() + 1; + constexpr uint64_t kBlockSize = 64 << 14; + const uint64_t num_blocks = (num_rows + kBlockSize - 1) / kBlockSize; + sdsl::bit_vector needs_char_bv(num_rows, 0); + std::vector> block_chars(num_blocks); + + #pragma omp parallel for num_threads(get_num_threads()) schedule(dynamic) + for (uint64_t b = 0; b < num_blocks; ++b) { + std::vector chars; + const uint64_t end = std::min(num_rows, (b + 1) * kBlockSize); + for (node_index node = std::max(b * kBlockSize, 1); node < end; ++node) { + // If the graph can trivially determine the successor, no need to store + if (!graph.in_graph(node) || graph.has_no_outgoing(node) + || graph.trivial_outgoing(node) != graph::DeBruijnGraph::npos) + continue; + + // Need to store the character for this node + node_index succ = row_diff_successor(graph, node, fork_succ); + int idx = -1; + graph.call_outgoing_kmers(node, [&](node_index adj, char c) { + if (adj == succ) + idx = char_to_idx[static_cast(c)]; + }); + if (idx < 0) { + logger->error("No outgoing edge of node {} leads to its row-diff " + "successor {}", node, succ); + std::exit(1); + } + needs_char_bv[node] = true; + chars.push_back(idx); + } + block_chars[b] = sdsl::int_vector<>(chars.size(), 0, bits_per_char); + std::copy(chars.begin(), chars.end(), block_chars[b].begin()); + } + + size_t count = 0; + for (const auto &chars : block_chars) { + count += chars.size(); + } + sdsl::int_vector<> rd_succ_char(count, 0, bits_per_char); + auto it = rd_succ_char.begin(); + for (auto &chars : block_chars) { + it = std::copy(chars.begin(), chars.end(), it); + chars = sdsl::int_vector<>(); + } + + logger->trace("Computed rd_succ_char: {} nodes total, {} need stored char ({:.1f}%), " + "{} bits/char", + graph.max_index(), count, + 100.0 * count / std::max(graph.max_index(), 1), + bits_per_char); + return { bit_vector_small(std::move(needs_char_bv)), std::move(rd_succ_char) }; +} + std::tuple>, LabelEncoder, uint64_t> load_row_diff_columns(const std::vector &files) { std::vector> columns; diff --git a/metagraph/src/annotation/annotation_converters.hpp b/metagraph/src/annotation/annotation_converters.hpp index ac6a3bc561..798948d743 100644 --- a/metagraph/src/annotation/annotation_converters.hpp +++ b/metagraph/src/annotation/annotation_converters.hpp @@ -139,6 +139,13 @@ void convert_row_diff_to_col_compressed(const std::vector &files, std::pair get_anchors_and_fork_fnames(const std::string &fbase); +// Compute sparse successor character data for CharSucc mode. +// Returns (needs_char bitvector, compact char array). Only nodes where the graph +// can't trivially determine the successor have stored characters. +std::pair> +compute_rd_succ_char(const graph::DeBruijnGraph &graph, + const bit_vector &fork_succ); + template StaticBinRelAnnotator, std::string> load_coords(Annotator&& anno, const std::vector &files); diff --git a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp index dcb44e3496..1995744965 100644 --- a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp +++ b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp @@ -38,6 +38,36 @@ node_index row_diff_successor(const graph::DeBruijnGraph &graph, } } +node_index row_diff_successor_char( + const graph::DeBruijnGraph &graph, + node_index node, + const bit_vector_small &rd_succ_needs_char, + const sdsl::int_vector<> &rd_succ_char) { + // Fast path: if the graph can trivially resolve the successor (e.g., + // succ_is_next_ in SSHash), use it directly without any char lookup. + node_index trivial = graph.trivial_outgoing(node); + if (trivial != graph::DeBruijnGraph::npos) + return trivial; + + // Slow path: look up stored character for this node + assert(node < rd_succ_needs_char.size()); + assert(rd_succ_needs_char[node] && "node must have a stored char"); + uint64_t rank = rd_succ_needs_char.rank1(node); + assert(rank > 0 && rank <= rd_succ_char.size()); + uint8_t idx = rd_succ_char[rank - 1]; + const auto &alphabet = graph.alphabet(); + size_t offset = (alphabet[0] == '$') ? 1 : 0; + assert(idx + offset < alphabet.size()); + char c = alphabet[idx + offset]; + node_index succ = graph.traverse(node, c); + if (succ == graph::DeBruijnGraph::npos) { + // called in parallel regions, where an exception would terminate anyway + common::logger->error("Corrupt rd_succ_char: traverse({}, '{}') returned npos", node, c); + std::exit(1); + } + return succ; +} + namespace matrix { @@ -124,7 +154,16 @@ IRowDiff::get_rd_ids(std::span row_ids, size_t num_thre if (anchor_[row]) break; - node = row_diff_successor(*graph_, node, fork_succ_); + switch (rd_mode_) { + case RowDiffMode::CharSucc: + node = row_diff_successor_char(*graph_, node, + rd_succ_needs_char_, + rd_succ_char_); + break; + case RowDiffMode::Standard: + node = row_diff_successor(*graph_, node, fork_succ_); + break; + } } } diff --git a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp index c559f616a1..fb01740c36 100644 --- a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp +++ b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp @@ -6,6 +6,8 @@ #include #include +#include + #include "annotation/binary_matrix/base/binary_matrix.hpp" #include "annotation/binary_matrix/column_sparse/column_major.hpp" #include "common/batch_accumulator.hpp" @@ -30,8 +32,20 @@ graph::DeBruijnGraph::node_index row_diff_successor(const graph::DeBruijnGraph & graph::DeBruijnGraph::node_index node, const bit_vector &rd_succ); +// CharSucc mode: look up successor from sparse char storage or trivial graph path +graph::DeBruijnGraph::node_index row_diff_successor_char( + const graph::DeBruijnGraph &graph, + graph::DeBruijnGraph::node_index node, + const bit_vector_small &rd_succ_needs_char, + const sdsl::int_vector<> &rd_succ_char); + namespace matrix { +enum class RowDiffMode : uint8_t { + Standard = 0, // original fork_succ bitvector + CharSucc = 1, // successor chars, stored only where trivial_outgoing() fails +}; + const std::string kRowDiffAnchorExt = ".anchors"; const std::string kRowDiffForkSuccExt = ".rd_succ"; @@ -53,6 +67,19 @@ class IRowDiff { const fork_succ_bv_type& fork_succ() const { return fork_succ_; } + const bit_vector_small& rd_succ_needs_char() const { return rd_succ_needs_char_; } + const sdsl::int_vector<>& rd_succ_char() const { return rd_succ_char_; } + + RowDiffMode mode() const { return rd_mode_; } + + // Switch to CharSucc mode: the successor chars replace fork_succ + void set_successor_chars(bit_vector_small&& needs_char, sdsl::int_vector<>&& succ_char) { + rd_succ_needs_char_ = std::move(needs_char); + rd_succ_char_ = std::move(succ_char); + fork_succ_ = fork_succ_bv_type(); + rd_mode_ = RowDiffMode::CharSucc; + } + protected: template Matrix> std::vector @@ -74,8 +101,11 @@ class IRowDiff { H decode_diffs, Callback call_row, size_t num_threads) const; const graph::DeBruijnGraph *graph_ = nullptr; + RowDiffMode rd_mode_ = RowDiffMode::Standard; anchor_bv_type anchor_; fork_succ_bv_type fork_succ_; + bit_vector_small rd_succ_needs_char_; + sdsl::int_vector<> rd_succ_char_; }; /** @@ -151,6 +181,23 @@ class RowDiff : public IRowDiff, public BinaryMatrix { }; +/** + * RowDiffChar: RowDiff variant that stores successor characters only for nodes + * whose successor trivial_outgoing() cannot resolve, instead of the 1-bit + * fork_succ bitvector. Enables single-lookup traversal + * via graph.traverse(node, char) instead of enumerating all outgoing neighbors. + */ +template +class RowDiffChar : public RowDiff { + public: + template + RowDiffChar(const graph::DeBruijnGraph *graph = nullptr, Args&&... args) + : RowDiff(graph, std::forward(args)...) { + this->rd_mode_ = RowDiffMode::CharSucc; + } +}; + + template Matrix> std::vector IRowDiff::get_column_from_rows(const Matrix &matrix, BinaryMatrix::Column column) const { @@ -296,7 +343,21 @@ template bool RowDiff::load(std::istream &f) { auto pos = f.tellg(); std::string version(4, '\0'); - if (f.read(version.data(), 4) && version == "v2.0") { + if (f.read(version.data(), 4) && version == "v3.0") { + // CharSucc mode: anchor + needs_char(bit_vector_small) + rd_succ_char(int_vector<>) + diffs + rd_mode_ = RowDiffMode::CharSucc; + if constexpr(!std::is_same_v) { + auto anchor_start = f.tellg(); + if (!anchor_.load(f) || !rd_succ_needs_char_.load(f)) + return false; + rd_succ_char_.load(f); + if (!f.good()) + return false; + // anchor_ / successor chars are accessed randomly, hint the kernel. + utils::madvise_random_range(f, anchor_start, f.tellg() - anchor_start); + } + } else if (version == "v2.0") { + rd_mode_ = RowDiffMode::Standard; if constexpr(!std::is_same_v) { auto anchor_start = f.tellg(); if (!anchor_.load(f) || !fork_succ_.load(f)) @@ -307,6 +368,7 @@ bool RowDiff::load(std::istream &f) { } else { // backward compatibility f.seekg(pos); + rd_mode_ = RowDiffMode::Standard; if constexpr(!std::is_same_v) { auto anchor_start = f.tellg(); if (!anchor_.load(f)) @@ -325,10 +387,23 @@ bool RowDiff::load(std::istream &f) { template void RowDiff::serialize(std::ostream &f) const { - f.write("v2.0", 4); - if constexpr(!std::is_same_v) { - anchor_.serialize(f); - fork_succ_.serialize(f); + switch (rd_mode_) { + case RowDiffMode::CharSucc: + f.write("v3.0", 4); + if constexpr(!std::is_same_v) { + anchor_.serialize(f); + rd_succ_needs_char_.serialize(f); + rd_succ_char_.serialize(f); + } + break; + case RowDiffMode::Standard: + default: + f.write("v2.0", 4); + if constexpr(!std::is_same_v) { + anchor_.serialize(f); + fork_succ_.serialize(f); + } + break; } diffs_.serialize(f); } diff --git a/metagraph/src/annotation/representation/annotation_matrix/annotation_matrix.cpp b/metagraph/src/annotation/representation/annotation_matrix/annotation_matrix.cpp index 8fb8049a3b..c82f57cd99 100644 --- a/metagraph/src/annotation/representation/annotation_matrix/annotation_matrix.cpp +++ b/metagraph/src/annotation/representation/annotation_matrix/annotation_matrix.cpp @@ -198,6 +198,8 @@ template class StaticBinRelAnnotator, std::string>; template class StaticBinRelAnnotator, std::string>; +template class StaticBinRelAnnotator, std::string>; + template class StaticBinRelAnnotator, std::string>; template class StaticBinRelAnnotator, std::string>; diff --git a/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp b/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp index 151054550d..e2d3c5b0cc 100644 --- a/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp +++ b/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp @@ -49,6 +49,8 @@ typedef StaticBinRelAnnotator, std::string> R typedef StaticBinRelAnnotator, std::string> RowDiffDiskAnnotator; +typedef StaticBinRelAnnotator, std::string> RowDiffCharBRWTAnnotator; + typedef sdsl::dac_vector_dp<> CountsVector; typedef StaticBinRelAnnotator, std::string> IntMultiBRWTAnnotator; @@ -95,6 +97,8 @@ inline const std::string RowDiffRowSparseAnnotator::kExtension = ".row_diff_spar template <> inline const std::string RowDiffDiskAnnotator::kExtension = ".row_diff_disk.annodbg"; template <> +inline const std::string RowDiffCharBRWTAnnotator::kExtension = ".row_diff_char_brwt.annodbg"; +template <> inline const std::string IntMultiBRWTAnnotator::kExtension = ".int_brwt.annodbg"; template <> inline const std::string IntRowDiffBRWTAnnotator::kExtension = ".row_diff_int_brwt.annodbg"; diff --git a/metagraph/src/cli/config/config.cpp b/metagraph/src/cli/config/config.cpp index 2d3346c7b8..87a13b1460 100644 --- a/metagraph/src/cli/config/config.cpp +++ b/metagraph/src/cli/config/config.cpp @@ -648,11 +648,14 @@ Config::Config(int argc, char *argv[]) { || anno_type == RowDiffRowSparse || anno_type == RowDiffDiskCoord || anno_type == RowDiffBRWTCoord - || anno_type == RowDiffCoord; + || anno_type == RowDiffCoord + || anno_type == RowDiffCharBRWT; // For row_diff_brwt and row_diff_flat inputs, the anchor and fork-succ // bitmaps are embedded in the .annodbg, so no graph is needed when // re-formatting to another row_diff representation. + // Computing successor characters for row_diff_char_brwt still needs the graph. const bool input_has_embedded_anchor = fnames.size() + && anno_type != RowDiffCharBRWT && (fnames.front().ends_with(annot::RowDiffBRWTAnnotator::kExtension) || fnames.front().ends_with(annot::RowDiffRowFlatAnnotator::kExtension)); if (input_has_embedded_anchor && fnames.size() != 1) { @@ -793,6 +796,8 @@ std::string Config::annotype_to_string(AnnotationType state) { return "row_sparse"; case RowDiffDisk: return "row_diff_disk"; + case RowDiffCharBRWT: + return "row_diff_char_brwt"; case IntRowDiffDisk: return "row_diff_int_disk"; case RowDiffDiskCoord: @@ -840,6 +845,8 @@ Config::AnnotationType Config::string_to_annotype(const std::string &string) { return AnnotationType::RowSparse; } else if (string == "row_diff_disk") { return AnnotationType::RowDiffDisk; + } else if (string == "row_diff_char_brwt") { + return AnnotationType::RowDiffCharBRWT; } else if (string == "row_diff_int_disk") { return AnnotationType::IntRowDiffDisk; } else if (string == "row_diff_disk_coord") { @@ -959,7 +966,7 @@ void Config::print_usage(const std::string &prog_name, IdentityType identity) { const char annotation_list[] = "\t\t( column, brwt, rb_brwt, int_brwt,\n" "\t\t column_coord, brwt_coord, row_diff_coord, row_diff_brwt_coord,\n" "\t\t row_diff, row_diff_brwt, row_diff_flat, row_diff_sparse, row_diff_int_brwt,\n" - "\t\t row_diff_disk, row_diff_int_disk, row_diff_disk_coord,\n" + "\t\t row_diff_disk, row_diff_int_disk, row_diff_disk_coord, row_diff_char_brwt,\n" "\t\t row, flat, row_sparse, rbfish, bin_rel_wt )"; switch (identity) { diff --git a/metagraph/src/cli/config/config.hpp b/metagraph/src/cli/config/config.hpp index 18cff74623..bbaa97abcd 100644 --- a/metagraph/src/cli/config/config.hpp +++ b/metagraph/src/cli/config/config.hpp @@ -212,6 +212,7 @@ class Config { RowDiffRowFlat, RowDiffRowSparse, RowDiffDisk, + RowDiffCharBRWT, RowFlat, RowSparse, RBFish, diff --git a/metagraph/src/cli/load/load_annotation.cpp b/metagraph/src/cli/load/load_annotation.cpp index 87bf7105ce..8c6b29b0a6 100644 --- a/metagraph/src/cli/load/load_annotation.cpp +++ b/metagraph/src/cli/load/load_annotation.cpp @@ -39,6 +39,9 @@ Config::AnnotationType parse_annotation_type(const std::string &filename) { } else if (filename.ends_with(annot::RowCompressed<>::kExtension)) { return Config::AnnotationType::RowCompressed; + } else if (filename.ends_with(annot::RowDiffCharBRWTAnnotator::kExtension)) { + return Config::AnnotationType::RowDiffCharBRWT; + } else if (filename.ends_with(annot::MultiBRWTAnnotator::kExtension)) { return Config::AnnotationType::BRWT; @@ -142,6 +145,10 @@ initialize_annotation(Config::AnnotationType anno_type, annotation.reset(new annot::RowDiffBRWTAnnotator()); break; } + case Config::RowDiffCharBRWT: { + annotation.reset(new annot::RowDiffCharBRWTAnnotator()); + break; + } case Config::RowDiffRowFlat: { annotation.reset(new annot::RowDiffRowFlatAnnotator()); break; diff --git a/metagraph/src/cli/transform_annotation.cpp b/metagraph/src/cli/transform_annotation.cpp index 343d7f618c..bb3593a72e 100644 --- a/metagraph/src/cli/transform_annotation.cpp +++ b/metagraph/src/cli/transform_annotation.cpp @@ -13,6 +13,7 @@ #include "annotation/annotation_converters.hpp" #include "config/config.hpp" #include "load/load_annotation.hpp" +#include "load/load_graph.hpp" namespace mtg { @@ -772,6 +773,10 @@ int transform_annotation(Config *config) { logger->error("Convert to row_diff first, and then to row_diff_disk"); return 0; } + case Config::RowDiffCharBRWT: { + logger->error("Convert to row_diff_brwt first, and then to row_diff_char_brwt"); + return 0; + } case Config::RowDiff: { auto out_dir = std::filesystem::path(config->outfbase).remove_filename(); convert_to_row_diff(files, config->infbase, config->memory_available * 1e9, @@ -937,6 +942,37 @@ int transform_annotation(Config *config) { } logger->trace("Serialized to {}", config->outfbase); + } else if (input_anno_type == Config::RowDiffBRWT && config->anno_type == Config::RowDiffCharBRWT) { + if (files.size() != 1) { + logger->error("Can only convert row_diff_brwt annotations one at a time"); + exit(1); + } + + logger->trace("Loading annotation..."); + RowDiffBRWTAnnotator annotator; + if (!annotator.load(files[0])) { + logger->error("Cannot load annotations from file '{}'", files[0]); + exit(1); + } + + logger->trace("Loading graph..."); + auto graph = load_critical_dbg(config->infbase); + + logger->trace("Computing rd_succ_char from fork_succ..."); + auto matrix = annotator.release_matrix(); + auto [needs_char, rd_succ_char] = compute_rd_succ_char(*graph, matrix->fork_succ()); + matrix->set_successor_chars(std::move(needs_char), std::move(rd_succ_char)); + + // BRWT diffs are identical to the original, RowDiff::serialize writes the v3.0 layout + const auto &fname = utils::make_suffix(config->outfbase, + RowDiffCharBRWTAnnotator::kExtension); + std::ofstream out = utils::open_new_ofstream(fname); + annotator.get_label_encoder().serialize(out); + matrix->serialize(out); + out.close(); + + logger->trace("Serialized to {}", config->outfbase); + } else { logger->error( "Conversion to other representations is not implemented for {} " diff --git a/metagraph/src/graph/representation/base/sequence_graph.hpp b/metagraph/src/graph/representation/base/sequence_graph.hpp index c0bf4bfb57..3b37e1e2c6 100644 --- a/metagraph/src/graph/representation/base/sequence_graph.hpp +++ b/metagraph/src/graph/representation/base/sequence_graph.hpp @@ -213,6 +213,10 @@ class DeBruijnGraph : public SequenceGraph { virtual bool has_single_outgoing(node_index node) const { return outdegree(node) == 1; } virtual bool has_multiple_outgoing(node_index node) const { return outdegree(node) > 1; } + // Return the unique outgoing neighbor if it can be determined in O(1) + // without expensive lookups (e.g., MPHF). Returns npos otherwise. + virtual node_index trivial_outgoing(node_index) const { return npos; } + virtual size_t indegree(node_index) const = 0; virtual bool has_no_incoming(node_index node) const { return indegree(node) == 0; } virtual bool has_single_incoming(node_index node) const { return indegree(node) == 1; } diff --git a/metagraph/src/graph/representation/hash/dbg_sshash.cpp b/metagraph/src/graph/representation/hash/dbg_sshash.cpp index 32f288a649..7c0cf3332a 100644 --- a/metagraph/src/graph/representation/hash/dbg_sshash.cpp +++ b/metagraph/src/graph/representation/hash/dbg_sshash.cpp @@ -391,6 +391,22 @@ void DBGSSHash::map_to_nodes_sequentially(std::string_view sequence, } } +DBGSSHash::node_index DBGSSHash::trivial_outgoing(node_index node) const { + assert(in_graph(node)); + if (node <= dict_size()) { + if (succ_is_next_[node]) + return node + 1; + } else { + // RC node v: forward counterpart is v - dict_size() + // pred_is_prev_[fwd] means the successor of v is the RC of fwd - 1, + // which is fwd - 1 itself if that k-mer is a palindrome (even k only) + node_index rev_comp = reverse_complement(node); + if (pred_is_prev_[rev_comp]) + return reverse_complement(rev_comp - 1); + } + return npos; +} + DBGSSHash::node_index DBGSSHash::traverse(node_index node, char next_char) const { assert(in_graph(node)); // Fast path: if successor is the next index, just check if the character matches diff --git a/metagraph/src/graph/representation/hash/dbg_sshash.hpp b/metagraph/src/graph/representation/hash/dbg_sshash.hpp index f635fc6574..3ee1d3dea1 100644 --- a/metagraph/src/graph/representation/hash/dbg_sshash.hpp +++ b/metagraph/src/graph/representation/hash/dbg_sshash.hpp @@ -116,6 +116,8 @@ class DBGSSHash : public DeBruijnGraph { const std::string& alphabet() const override final { return alphabet_; } + node_index trivial_outgoing(node_index node) const override; + const dict_t& data() const { return dict_; } node_index reverse_complement(node_index node) const; diff --git a/metagraph/tests/graph/all/test_dbg_node_degree.cpp b/metagraph/tests/graph/all/test_dbg_node_degree.cpp index 7dd7fd1fd1..6ecc8526c6 100644 --- a/metagraph/tests/graph/all/test_dbg_node_degree.cpp +++ b/metagraph/tests/graph/all/test_dbg_node_degree.cpp @@ -319,4 +319,41 @@ TYPED_TEST(DeBruijnGraphTest, is_single_outgoing_for_multiple_valid_edges) { EXPECT_EQ(1u, single_outgoing_counter); } + +// trivial_outgoing() may only shortcut nodes with a single outgoing edge +TYPED_TEST(DeBruijnGraphTest, trivial_outgoing_is_unique_successor) { + // small k give highly branching graphs, so count shortcuts over all k + size_t num_trivial = 0; + for (size_t k = 3; k < 10; k += 2) { + auto graph = build_graph(k, { + "ATGCGATCGATATGCGAGA", + "ATGCGATCGAGACTACGAG", + "GTACGATAGACATGACGAG", + "ACTGACGAGACACAGATGC" + }); + + graph->call_nodes([&](auto node) { + auto succ = graph->trivial_outgoing(node); + if (succ == DeBruijnGraph::npos) + return; + + ++num_trivial; + // look up successors by k-mer, independently of any index shortcuts + std::string suffix = graph->get_node_sequence(node).substr(1); + std::vector outgoing; + for (char c : graph->alphabet()) { + auto next = graph->kmer_to_node(suffix + c); + if (next != DeBruijnGraph::npos) + outgoing.push_back(next); + } + EXPECT_EQ(std::vector{ succ }, outgoing) + << "k=" << k << ", node: " << node; + }); + } + + if constexpr(std::is_same_v) { + EXPECT_LT(0u, num_trivial); + } +} + } // namespace diff --git a/metagraph/tests/graph/test_dbg_canonical.cpp b/metagraph/tests/graph/test_dbg_canonical.cpp index 55e5c38ffe..f38a66b4d5 100644 --- a/metagraph/tests/graph/test_dbg_canonical.cpp +++ b/metagraph/tests/graph/test_dbg_canonical.cpp @@ -1465,6 +1465,51 @@ TYPED_TEST(DeBruijnGraphCanonicalTest, CallUnitigsCross) { } } + +// trivial_outgoing() must also resolve reverse-complement nodes correctly +TYPED_TEST(DeBruijnGraphCanonicalTest, trivial_outgoing_is_unique_successor) { + // small k give highly branching graphs, so count shortcuts over all k; + // even k cover palindromic k-mers (e.g. ACGT, TACGTA, TTACGTAA) + size_t num_trivial = 0; + for (size_t k = 3; k < 10; ++k) { + // canonical SSHash graphs support only odd k, so they have no palindromes + if constexpr(std::is_same_v) { + if (k % 2 == 0) + continue; + } + + auto graph = build_graph(k, { + "ATGCGATCGATATGCGAGA", + "ATGCGATCGAGACTACGAG", + "GTACGATAGACATGACGAG", + "ACTGACGAGACACAGATGC", + "GCAATTACGTAATTGCC" + }, DeBruijnGraph::CANONICAL); + + graph->call_nodes([&](auto node) { + auto succ = graph->trivial_outgoing(node); + if (succ == DeBruijnGraph::npos) + return; + + ++num_trivial; + // look up successors by k-mer, independently of any index shortcuts + std::string suffix = graph->get_node_sequence(node).substr(1); + std::vector outgoing; + for (char c : graph->alphabet()) { + auto next = graph->kmer_to_node(suffix + c); + if (next != DeBruijnGraph::npos) + outgoing.push_back(next); + } + EXPECT_EQ(std::vector{ succ }, outgoing) + << "k=" << k << ", node: " << node; + }); + } + + if constexpr(std::is_same_v) { + EXPECT_LT(0u, num_trivial); + } +} + #endif // ! _PROTEIN_GRAPH } // namespace