From e9ccb1469138fbfda5c87f0a3ae79d44e360ff9d Mon Sep 17 00:00:00 2001 From: Harun Mustafa Date: Sun, 13 Oct 2024 15:33:12 +0200 Subject: [PATCH 001/308] initial commit of alignment redo --- .../graph/alignment_redone/aligner_config.cpp | 255 ++++++++++ .../graph/alignment_redone/aligner_config.hpp | 99 ++++ .../graph/alignment_redone/aln_chainer.hpp | 252 ++++++++++ .../src/graph/alignment_redone/aln_cigar.cpp | 240 ++++++++++ .../src/graph/alignment_redone/aln_cigar.hpp | 114 +++++ .../src/graph/alignment_redone/aln_match.cpp | 123 +++++ .../src/graph/alignment_redone/aln_match.hpp | 152 ++++++ .../src/graph/alignment_redone/aln_query.hpp | 34 ++ .../src/graph/alignment_redone/aln_seeder.cpp | 434 ++++++++++++++++++ .../src/graph/alignment_redone/aln_seeder.hpp | 32 ++ .../alignment_redone/annotation_buffer.cpp | 221 +++++++++ .../alignment_redone/annotation_buffer.hpp | 85 ++++ 12 files changed, 2041 insertions(+) create mode 100644 metagraph/src/graph/alignment_redone/aligner_config.cpp create mode 100644 metagraph/src/graph/alignment_redone/aligner_config.hpp create mode 100644 metagraph/src/graph/alignment_redone/aln_chainer.hpp create mode 100644 metagraph/src/graph/alignment_redone/aln_cigar.cpp create mode 100644 metagraph/src/graph/alignment_redone/aln_cigar.hpp create mode 100644 metagraph/src/graph/alignment_redone/aln_match.cpp create mode 100644 metagraph/src/graph/alignment_redone/aln_match.hpp create mode 100644 metagraph/src/graph/alignment_redone/aln_query.hpp create mode 100644 metagraph/src/graph/alignment_redone/aln_seeder.cpp create mode 100644 metagraph/src/graph/alignment_redone/aln_seeder.hpp create mode 100644 metagraph/src/graph/alignment_redone/annotation_buffer.cpp create mode 100644 metagraph/src/graph/alignment_redone/annotation_buffer.hpp diff --git a/metagraph/src/graph/alignment_redone/aligner_config.cpp b/metagraph/src/graph/alignment_redone/aligner_config.cpp new file mode 100644 index 0000000000..ec6d04e857 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aligner_config.cpp @@ -0,0 +1,255 @@ +#include "aligner_config.hpp" + +#include "aln_cigar.hpp" +#include "kmer/alphabets.hpp" +#include "common/logger.hpp" + + +namespace mtg::graph::align { + +using mtg::common::logger; + +void DBGAlignerConfig::print_summary() const { + logger->trace("Alignment settings:"); + logger->trace("\t Alignments to report: {}", num_alternative_paths); + logger->trace("\t Min seed length: {}", min_seed_length); + logger->trace("\t Max seed length: {}", max_seed_length == std::numeric_limits::max() + ? "inf" : std::to_string(max_seed_length)); + logger->trace("\t Max num seeds per locus: {}", max_num_seeds_per_locus); + logger->trace("\t Max num nodes per sequence char: {}", max_nodes_per_seq_char); + logger->trace("\t Max RAM per alignment: {}", max_ram_per_alignment); + logger->trace("\t Gap opening penalty: {}", (int64_t)gap_opening_penalty); + logger->trace("\t Gap extension penalty: {}", (int64_t)gap_extension_penalty); + logger->trace("\t Min alignment score: {}", min_path_score); + logger->trace("\t X drop-off: {}", xdrop); + logger->trace("\t Exact nucleotide match threshold: {}", min_exact_match); + logger->trace("\t Chain alignments: {}", chain_alignments); + + logger->trace("\t Scoring matrix: {}", alignment_edit_distance ? "unit costs" : "matrix"); + if (!alignment_edit_distance) { + logger->trace("\t\t Match score: {}", (int64_t)alignment_match_score); + logger->trace("\t\t (DNA) Transition score: {}", (int64_t)alignment_mm_transition_score); + logger->trace("\t\t (DNA) Transversion score: {}", (int64_t)alignment_mm_transversion_score); + } +} + +// check to make sure the current scoring system won't underflow +bool DBGAlignerConfig::check_config_scores() const { + int8_t min_penalty_score = std::numeric_limits::max(); + for (const auto &row : score_matrix) { + min_penalty_score = std::min(min_penalty_score, + *std::min_element(row.begin(), row.end())); + } + + if (gap_opening_penalty * 2 >= min_penalty_score) { + common::logger->error( + "|gap_opening_penalty| * 2 should be greater than greatest mismatch penalty: " + "{} >= {}", + gap_opening_penalty * 2, (score_t)min_penalty_score + ); + return false; + } + + min_penalty_score = std::min({ min_penalty_score, + gap_opening_penalty, + gap_extension_penalty }); + + assert(min_penalty_score < 0 && "min scores must be negative"); + + if (min_cell_score >= std::numeric_limits::min() - min_penalty_score) + return true; + + std::cerr << "min_cell_score is too small: " << min_cell_score << std::endl; + return false; +} + +DBGAlignerConfig::score_t DBGAlignerConfig +::score_cigar(std::string_view reference, + std::string_view query, + const Cigar &cigar) const { + assert(cigar.is_valid(reference, query)); + + if (cigar.empty()) + return 0; + + score_t score = (!cigar.get_clipping() ? left_end_bonus : 0) + + (!cigar.get_end_clipping() ? right_end_bonus : 0); + + auto ref_it = reference.begin(); + auto alt_it = query.begin(); + auto it = cigar.data().begin(); + if (it->first == Cigar::CLIPPED) + ++it; + + for ( ; it != cigar.data().end(); ++it) { + const auto &op = *it; + switch (op.first) { + case Cigar::CLIPPED: { + if (it + 1 != cigar.data().end()) + alt_it += op.second; + } break; + case Cigar::MATCH: { + score += match_score(std::string_view(ref_it, op.second)); + ref_it += op.second; + alt_it += op.second; + } break; + case Cigar::MISMATCH: { + score += score_sequences(std::string_view(ref_it, op.second), + std::string_view(alt_it, op.second)); + ref_it += op.second; + alt_it += op.second; + } break; + case Cigar::INSERTION: { + score += gap_opening_penalty + (op.second - 1) * gap_extension_penalty; + alt_it += op.second; + } break; + case Cigar::DELETION: { + score += gap_opening_penalty + (op.second - 1) * gap_extension_penalty; + ref_it += op.second; + if (it >= cigar.data().begin() + 2 && (it - 2)->first == Cigar::DELETION + && (it - 1)->first == Cigar::NODE_INSERTION) { + score -= gap_opening_penalty - gap_extension_penalty; + } + } break; + case Cigar::NODE_INSERTION: { + score += gap_opening_penalty + (op.second - 1) * gap_extension_penalty; + } break; + } + } + + assert(ref_it == reference.end()); + assert(alt_it == query.end()); + + return score; +} + +void DBGAlignerConfig::set_scoring_matrix() { + if (alignment_edit_distance) { + // TODO: REPLACE THIS + #if _PROTEIN_GRAPH + const auto *alphabet = kmer::alphabets::kAlphabetProtein; + const auto *alphabet_encoding = kmer::alphabets::kCharToProtein; + #elif _DNA_GRAPH || _DNA5_GRAPH || _DNA_CASE_SENSITIVE_GRAPH + const auto *alphabet = kmer::alphabets::kAlphabetDNA; + const auto *alphabet_encoding = kmer::alphabets::kCharToDNA; + #else + static_assert(false, + "Define an alphabet: either " + "_DNA_GRAPH, _DNA5_GRAPH, _PROTEIN_GRAPH, or _DNA_CASE_SENSITIVE_GRAPH." + ); + #endif + + score_matrix = unit_scoring_matrix(1, alphabet, alphabet_encoding); + left_end_bonus = 0; + right_end_bonus = 0; + + } else { + #if _PROTEIN_GRAPH + score_matrix = score_matrix_blosum62; + #elif _DNA_GRAPH || _DNA5_GRAPH || _DNA_CASE_SENSITIVE_GRAPH + score_matrix = dna_scoring_matrix(alignment_match_score, + -alignment_mm_transition_score, + -alignment_mm_transversion_score); + #else + static_assert(false, + "Define an alphabet: either " + "_DNA_GRAPH, _DNA5_GRAPH, _PROTEIN_GRAPH, or _DNA_CASE_SENSITIVE_GRAPH." + ); + #endif + } +} + +DBGAlignerConfig::ScoreMatrix DBGAlignerConfig +::dna_scoring_matrix(int8_t match_score, + int8_t mm_transition_score, + int8_t mm_transversion_score) { + ScoreMatrix score_matrix; + for (auto& row : score_matrix) { + row.fill(mm_transversion_score); + } + + score_matrix['A']['G'] = mm_transition_score; + score_matrix['G']['A'] = mm_transition_score; + score_matrix['C']['T'] = mm_transition_score; + score_matrix['T']['C'] = mm_transition_score; + score_matrix['A']['A'] = match_score; + score_matrix['C']['C'] = match_score; + score_matrix['G']['G'] = match_score; + score_matrix['T']['T'] = match_score; + + return score_matrix; +} + +DBGAlignerConfig::ScoreMatrix DBGAlignerConfig +::unit_scoring_matrix(int8_t match_score, + const std::string &alphabet, + const uint8_t *encoding) { + ScoreMatrix score_matrix; + for (auto& row : score_matrix) { + row.fill(-match_score); + } + + for (uint8_t c : alphabet) { + // if a character is invalid, don't count matches of that character + if (encoding[c] == encoding[0]) + continue; + + char upper = toupper(c); + + score_matrix[upper][upper] = match_score; + } + + return score_matrix; +} + +DBGAlignerConfig::ScoreMatrix blosum62_scoring_matrix() { + std::string alphabet = "ARNDCQEGHILKMFPSTWYVBZX"; + + std::vector> scores = { + { 4, -1, -2, -2, 0, -1, -1, 0, -2, -1, -1, -1, -1, -2, -1, 1, 0, -3, -2, 0, -2, -1, 0 }, + { -1, 5, 0, -2, -3, 1, 0, -2, 0, -3, -2, 2, -1, -3, -2, -1, -1, -3, -2, -3, -1, 0, -1 }, + { -2, 0, 6, 1, -3, 0, 0, 0, 1, -3, -3, 0, -2, -3, -2, 1, 0, -4, -2, -3, 3, 0, -1 }, + { -2, -2, 1, 6, -3, 0, 2, -1, -1, -3, -4, -1, -3, -3, -1, 0, -1, -4, -3, -3, 4, 1, -1 }, + { 0, -3, -3, -3, 9, -3, -4, -3, -3, -1, -1, -3, -1, -2, -3, -1, -1, -2, -2, -1, -3, -3, -2 }, + { -1, 1, 0, 0, -3, 5, 2, -2, 0, -3, -2, 1, 0, -3, -1, 0, -1, -2, -1, -2, 0, 3, -1 }, + { -1, 0, 0, 2, -4, 2, 5, -2, 0, -3, -3, 1, -2, -3, -1, 0, -1, -3, -2, -2, 1, 4, -1 }, + { 0, -2, 0, -1, -3, -2, -2, 6, -2, -4, -4, -2, -3, -3, -2, 0, -2, -2, -3, -3, -1, -2, -1 }, + { -2, 0, 1, -1, -3, 0, 0, -2, 8, -3, -3, -1, -2, -1, -2, -1, -2, -2, 2, -3, 0, 0, -1 }, + { -1, -3, -3, -3, -1, -3, -3, -4, -3, 4, 2, -3, 1, 0, -3, -2, -1, -3, -1, 3, -3, -3, -1 }, + { -1, -2, -3, -4, -1, -2, -3, -4, -3, 2, 4, -2, 2, 0, -3, -2, -1, -2, -1, 1, -4, -3, -1 }, + { -1, 2, 0, -1, -3, 1, 1, -2, -1, -3, -2, 5, -1, -3, -1, 0, -1, -3, -2, -2, 0, 1, -1 }, + { -1, -1, -2, -3, -1, 0, -2, -3, -2, 1, 2, -1, 5, 0, -2, -1, -1, -1, -1, 1, -3, -1, -1 }, + { -2, -3, -3, -3, -2, -3, -3, -3, -1, 0, 0, -3, 0, 6, -4, -2, -2, 1, 3, -1, -3, -3, -1 }, + { -1, -2, -2, -1, -3, -1, -1, -2, -2, -3, -3, -1, -2, -4, 7, -1, -1, -4, -3, -2, -2, -1, -2 }, + { 1, -1, 1, 0, -1, 0, 0, 0, -1, -2, -2, 0, -1, -2, -1, 4, 1, -3, -2, -2, 0, 0, 0 }, + { 0, -1, 0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1, 1, 5, -2, -2, 0, -1, -1, 0 }, + { -3, -3, -4, -4, -2, -2, -3, -2, -2, -3, -2, -3, -1, 1, -4, -3, -2, 11, 2, -3, -4, -3, -2 }, + { -2, -2, -2, -3, -2, -1, -2, -3, 2, -1, -1, -2, -1, 3, -3, -2, -2, 2, 7, -1, -3, -2, -1 }, + { 0, -3, -3, -3, -1, -2, -2, -3, -3, 3, 1, -2, 1, -1, -2, -2, 0, -3, -1, 4, -3, -2, -1 }, + { -2, -1, 3, 4, -3, 0, 1, -1, 0, -3, -4, 0, -3, -3, -2, 0, -1, -4, -3, -3, 4, 1, -1 }, + { -1, 0, 0, 1, -3, 3, 4, -2, 0, -3, -3, 1, -1, -3, -1, 0, -1, -3, -2, -2, 1, 4, -1 }, + { 0, -1, -1, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, 0, 0, -2, -1, -1, -1, -1, -1 } + }; + + DBGAlignerConfig::ScoreMatrix score_matrix; + + for (size_t i = 0; i < score_matrix.size(); ++i) { + score_matrix[i].fill(-4); + + // meant to handle the letters J, O, U + score_matrix[i][i] = 1; + } + + for (size_t i = 0; i < alphabet.size(); ++i) { + for (size_t j = 0; j < alphabet.size(); ++j) { + score_matrix[alphabet[i]][alphabet[j]] = scores[i][j]; + } + } + + return score_matrix; +} + +const DBGAlignerConfig::ScoreMatrix DBGAlignerConfig::score_matrix_blosum62 + = blosum62_scoring_matrix(); + +} // namespace mtg::graph::align diff --git a/metagraph/src/graph/alignment_redone/aligner_config.hpp b/metagraph/src/graph/alignment_redone/aligner_config.hpp new file mode 100644 index 0000000000..f372b3bb38 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aligner_config.hpp @@ -0,0 +1,99 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + + +namespace mtg::graph::align { + +class Cigar; + +struct DBGAlignerConfig { + typedef int32_t score_t; + typedef std::array ScoreMatrixRow; + typedef std::array ScoreMatrix; + + size_t num_alternative_paths = 1; + size_t min_seed_length = 0; + size_t max_seed_length = 0; + size_t max_num_seeds_per_locus = std::numeric_limits::max(); + + size_t max_dist_between_seeds = 400; + size_t gap_shrinking_factor = 4; + double chaining_algorithm_switch_cutoff = 1.0; + double chains_per_char = 0; + size_t chaining_bandwidth = 65; + + // Lowest possible score. 100 is added to prevent underflow during operations. + // For this to work, all penalties should be less than 100. + // This is checked whenever an aligner is initialized. + constexpr static score_t ninf = std::numeric_limits::min() + 100; + + // thresholds for scores + score_t min_cell_score = ninf; + score_t min_path_score = 0; + score_t xdrop = std::numeric_limits::max(); + + double min_exact_match = 0.0; + double max_nodes_per_seq_char = std::numeric_limits::max(); + double max_ram_per_alignment = std::numeric_limits::max(); + double rel_score_cutoff = 0.0; + + int8_t gap_opening_penalty = -5; + int8_t gap_extension_penalty = -2; + int8_t left_end_bonus = 0; + int8_t right_end_bonus = 0; + + bool forward_and_reverse_complement = true; + bool chain_alignments = false; + bool post_chain_alignments = false; + bool global_xdrop = true; + bool allow_left_trim = true; + bool no_backtrack = false; + bool seed_complexity_filter = true; + + bool alignment_edit_distance; + int8_t alignment_match_score; + int8_t alignment_mm_transition_score; + int8_t alignment_mm_transversion_score; + + ScoreMatrix score_matrix; + + void print_summary() const; + + score_t score_sequences(std::string_view a, std::string_view b) const { + return std::inner_product( + a.begin(), a.end(), b.begin(), score_t(0), std::plus(), + [&](char a, char b) -> score_t { return score_matrix[a][b]; } + ); + } + + score_t match_score(std::string_view query) const { + return score_sequences(query, query); + } + + score_t score_cigar(std::string_view reference, + std::string_view query, + const Cigar &cigar) const; + + bool check_config_scores() const; + + void set_scoring_matrix(); + + // Protein matrices + static const ScoreMatrix score_matrix_blosum62; + + static ScoreMatrix dna_scoring_matrix(int8_t match_score, + int8_t mm_transition_score, + int8_t mm_transversion_score); + + static ScoreMatrix unit_scoring_matrix(int8_t match_score, + const std::string &alphabet, + const uint8_t *encoding); +}; + +} // namespace mtg::graph::align diff --git a/metagraph/src/graph/alignment_redone/aln_chainer.hpp b/metagraph/src/graph/alignment_redone/aln_chainer.hpp new file mode 100644 index 0000000000..cd22bfb356 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_chainer.hpp @@ -0,0 +1,252 @@ +#pragma once + +#include "aligner_config.hpp" +#include "aln_query.hpp" +#include "aln_match.hpp" + +namespace mtg::graph::align { + +template +using ChainScores = std::vector>; + +template +using ScoreUpdater = std::function; + +template +using AnchorChain = std::vector>; + +template +using ExtensionStarter = std::function&, // chain + const std::vector& // score traceback + )>; + + +template +using AnchorConnector = std::function::iterator, + const ScoreUpdater& + )>; + +using AlignmentCallback = std::function; + +template +using AnchorExtender = std::function; + +template +void extend_chain(const AnchorChain &chain, + const std::vector &score_traceback, + const AnchorExtender &anchor_extender, + const AlignmentCallback &callback, + const std::function &terminate = []() { return false; }) { + std::vector alns; + alns.emplace_back(*chain.back().first); + + for (auto it = chain.rbegin() + 1; it != chain.rend(); ++it) { + std::vector next_alns; + for (auto&& aln : alns) { + size_t dist = (it - 1)->second; + anchor_extender((it - 1)->first, it->first, std::move(aln), dist, 0, + [&](Alignment&& next_aln) { next_alns.emplace_back(std::move(next_aln)); } + ); + } + std::swap(next_alns, alns); + } + // auto jt = score_traceback.begin(); + // std::vector alns; + // alns.emplace_back(*chain[0].first); + // ++jt; + + // for (auto it = chain.begin() + 1; it != chain.end(); ++it, ++jt) { + // assert(jt != score_traceback.end()); + // std::vector next_alns; + // for (auto&& aln : alns) { + // anchor_extender(it->first, (it - 1)->first, std::move(aln), it->second, *jt, + // [&](Alignment&& next_aln) { + // next_alns.emplace_back(std::move(next_aln)); + // } + // ); + // } + // std::swap(next_alns, alns); + // } + + assert(jt == score_traceback.end()); + + for (auto&& aln : alns) { + if (terminate()) + return; + + callback(std::move(aln)); + } +} + +template +void chain_anchors(const Query &query, + const DBGAlignerConfig &config, + AnchorIt begin, + AnchorIt end, + const AnchorConnector &anchor_connector, + const ExtensionStarter &extension_starter, + const AnchorExtender &anchor_extender, + const AlignmentCallback &callback = [](Alignment&&) {}, + const std::function &terminate = []() { return false; }) { + if (terminate() || begin == end) + return; + + ssize_t query_size = query.get_query().size(); + + std::sort(begin, end, [&](const Anchor &a, const Anchor &b) { + return std::make_tuple(a.get_orientation(), a.get_label_class(), b.get_end_clipping()) + < std::make_tuple(b.get_orientation(), b.get_label_class(), a.get_end_clipping()); + }); + + ChainScores chain_scores; + chain_scores.reserve(end - begin); + std::transform(begin, end, std::back_inserter(chain_scores), [&](const Anchor &a) { + return std::make_tuple(a.get_score(), end, 0); + }); + + // forward pass + ssize_t max_gap_between_anchors = std::min( + static_cast(config.max_dist_between_seeds), + query_size + ); + + auto forward_pass = [&](AnchorIt begin, + AnchorIt end, + typename ChainScores::iterator chain_scores) { + if (begin == end) + return; + + auto make_anchor_connector = [&](const AnchorIt i) { + return [&,i](DBGAlignerConfig::score_t score, const AnchorIt last, size_t dist) { + assert(last != i); + auto &[max_score, best_last, best_dist] = *(chain_scores + (i - begin)); + if (score > max_score) { + max_score = score; + best_last = last; + best_dist = dist; + return true; + } else { + return false; + } + }; + }; + + if (static_cast(end - begin) / query_size + <= config.chaining_algorithm_switch_cutoff) { + // if there are fewer seeds, this algorithm is faster + ssize_t b = max_gap_between_anchors; + ssize_t b_last; + DBGAlignerConfig::score_t best_score = std::get<0>(chain_scores[0]); + ssize_t best_cost = std::numeric_limits::max(); + DBGAlignerConfig::score_t match_score = config.match_score("A"); + do { + AnchorIt j = begin; + for (AnchorIt i = begin; i != end; ++i) { + j = std::find_if(j, end, [&](const Anchor &a) { + return a.get_end_clipping() <= b + i->get_end_clipping(); + }); + + // align anchors [j, i) to i + anchor_connector(*i, b, j, i, chain_scores + (j - begin), + make_anchor_connector(i)); + best_score = std::max(best_score, std::get<0>(chain_scores[i - begin])); + } + b_last = b; + b *= config.gap_shrinking_factor; + + // best_score = 1/2(match_score(query_size + match_spelling_size - best_cost)) + // => 2*best_score/match_score - query_size - match_spelling_size = -best_cost + // => best_cost = match_spelling_size + query_size - 2*best_score/match_score + + // if we assume that match_spelling_size == query_size + best_cost = 2*(query_size - best_score/match_score); + } while (best_cost > b_last); + } else { + // otherwise, use this algorithm + for (AnchorIt i = begin + 1; i != end; ++i) { + AnchorIt j = i - std::min(static_cast(i - begin), config.chaining_bandwidth); + anchor_connector(*i, std::numeric_limits::max(), j, i, + chain_scores + (j - begin), + make_anchor_connector(i)); + } + } + }; + + auto cur_begin = begin; + auto jt = chain_scores.begin(); + for (auto it = begin + 1; it != end; ++it) { + if (it->get_orientation() != (it - 1)->get_orientation() + || it->get_label_class() != (it - 1)->get_label_class()) { + forward_pass(cur_begin, it, jt); + jt += it - cur_begin; + cur_begin = it; + } + } + + forward_pass(cur_begin, end, jt); + + // backtracking + std::vector> best_chains; + best_chains.reserve(chain_scores.size()); + for (size_t i = 0; i < chain_scores.size(); ++i) { + const auto &[score, last, dist] = chain_scores[i]; + + if (score > 0) + best_chains.emplace_back(-score, (begin + i)->get_orientation(), i); + } + + std::sort(best_chains.begin(), best_chains.end()); + + sdsl::bit_vector used(chain_scores.size()); + for (auto [nscore, orientation, i] : best_chains) { + if (terminate()) + return; + + if (used[i]) + continue; + + AnchorChain chain; + std::vector scores; + AnchorIt last_anchor = begin + i; + auto [score, last, dist] = chain_scores[i]; + chain.emplace_back(last_anchor, dist); + assert(score == -nscore); + scores.emplace_back(score); + while (last != end) { + last_anchor = last; + std::tie(score, last, dist) = chain_scores[last - begin]; + chain.emplace_back(last_anchor, dist); + scores.emplace_back(score); + if (used[last_anchor - begin]) { + chain.clear(); + scores.clear(); + break; + } + } + + std::reverse(chain.begin(), chain.end()); + std::reverse(scores.begin(), scores.end()); + + if (chain.size() && extension_starter(chain, scores)) { + for (const auto &[a_ptr, dist] : chain) { + used[a_ptr - begin] = true; + } + + extend_chain(chain, scores, anchor_extender, callback, terminate); + } + } +} + +} // namespace mtg::graph::align \ No newline at end of file diff --git a/metagraph/src/graph/alignment_redone/aln_cigar.cpp b/metagraph/src/graph/alignment_redone/aln_cigar.cpp new file mode 100644 index 0000000000..a330ab41d8 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_cigar.cpp @@ -0,0 +1,240 @@ +#include "aln_cigar.hpp" + +#include "kmer/alphabets.hpp" + +namespace mtg::graph::align { + + +OperatorTable initialize_opt_table() { + OperatorTable char_to_op; + + // TODO: fix this when alphabets are no longer set at compile time + #if _PROTEIN_GRAPH + const auto *alphabet = kmer::alphabets::kAlphabetProtein; + const auto *alphabet_encoding = kmer::alphabets::kCharToProtein; + #elif _DNA_CASE_SENSITIVE_GRAPH + const auto *alphabet = kmer::alphabets::kAlphabetDNA; + const auto *alphabet_encoding = kmer::alphabets::kCharToDNA; + #elif _DNA5_GRAPH + const auto *alphabet = kmer::alphabets::kAlphabetDNA; + const auto *alphabet_encoding = kmer::alphabets::kCharToDNA; + #elif _DNA_GRAPH + const auto *alphabet = kmer::alphabets::kAlphabetDNA; + const auto *alphabet_encoding = kmer::alphabets::kCharToDNA; + #else + static_assert(false, + "Define an alphabet: either " + "_DNA_GRAPH, _DNA5_GRAPH, _PROTEIN_GRAPH, or _DNA_CASE_SENSITIVE_GRAPH." + ); + #endif + + for (auto &row : char_to_op) { + row.fill(Cigar::MISMATCH); + } + + for (uint8_t c : std::string(alphabet)) { + if (alphabet_encoding[c] == alphabet_encoding[0]) + continue; + + char upper = toupper(c); + char lower = tolower(c); + + char_to_op[upper][upper] + = char_to_op[upper][lower] + = char_to_op[lower][upper] + = char_to_op[lower][lower] = Cigar::MATCH; + } + + return char_to_op; +} + +Cigar::Cigar(std::string_view cigar_str) { + std::string op_count; + for (char c : cigar_str) { + switch (c) { + case '=': + cigar_.emplace_back(Cigar::MATCH, std::stol(op_count)); + op_count.clear(); + break; + case 'X': + cigar_.emplace_back(Cigar::MISMATCH, std::stol(op_count)); + op_count.clear(); + break; + case 'I': + cigar_.emplace_back(Cigar::INSERTION, std::stol(op_count)); + op_count.clear(); + break; + case 'D': + cigar_.emplace_back(Cigar::DELETION, std::stol(op_count)); + op_count.clear(); + break; + case 'S': + cigar_.emplace_back(Cigar::CLIPPED, std::stol(op_count)); + op_count.clear(); + break; + case 'G': + cigar_.emplace_back(Cigar::NODE_INSERTION, std::stol(op_count)); + break; + default: + op_count += c; + } + } +} + +std::string Cigar::to_string() const { + std::string cigar_string; + + for (const auto &pair : cigar_) { + cigar_string += std::to_string(pair.second) + opt_to_char(pair.first); + } + + return cigar_string; +} + +void Cigar::append(Operator op, LengthType num) { + if (!num) + return; + + if (cigar_.empty() || cigar_.back().first != op) { + cigar_.emplace_back(op, num); + } else { + cigar_.back().second += num; + } +} + +void Cigar::append(Cigar&& other) { + if (other.empty()) + return; + + append(other.cigar_.front().first, other.cigar_.front().second); + cigar_.insert(cigar_.end(), std::next(other.cigar_.begin()), other.cigar_.end()); +} + +bool Cigar::is_valid(std::string_view reference, std::string_view query) const { + auto ref_it = reference.begin(); + auto alt_it = query.begin(); + + for (size_t i = 0; i < cigar_.size(); ++i) { + const auto &op = cigar_[i]; + if (!op.second) { + std::cerr << "Empty operation found in CIGAR" << std::endl + << to_string() << std::endl + << reference << std::endl + << query << std::endl; + return false; + } + + switch (op.first) { + case CLIPPED: { + if ((ref_it != reference.begin() || alt_it != query.begin()) + && (ref_it != reference.end() || alt_it != query.end())) { + if (alt_it > query.end() - op.second) { + std::cerr << "Query too short after " + << Cigar::opt_to_char(op.first) << std::endl + << to_string() << std::endl + << reference << std::endl + << query << std::endl; + return false; + } + + alt_it += op.second; + } + } break; + case MATCH: + case MISMATCH: { + if (ref_it > reference.end() - op.second) { + std::cerr << "Reference too short after " + << Cigar::opt_to_char(op.first) << std::endl + << to_string() << std::endl + << reference << std::endl + << query << std::endl; + return false; + } + + if (alt_it > query.end() - op.second) { + std::cerr << "Query too short after " + << Cigar::opt_to_char(op.first) << std::endl + << to_string() << std::endl + << reference << std::endl + << query << std::endl; + return false; + } + + if (std::equal(ref_it, ref_it + op.second, alt_it) + == (op.first != Cigar::MATCH)) { + std::cerr << "Mismatch despite MATCH in CIGAR" << std::endl + << to_string() << std::endl + << reference << std::endl + << query << std::endl; + return false; + } + + ref_it += op.second; + alt_it += op.second; + } break; + case INSERTION: { + if (i && cigar_[i - 1].first == DELETION) { + std::cerr << "INSERTION after DELETION" << std::endl + << to_string() << std::endl + << reference << std::endl + << query << std::endl; + return false; + } + + if (alt_it > query.end() - op.second) { + std::cerr << "Query too short after " + << Cigar::opt_to_char(op.first) << std::endl + << to_string() << std::endl + << reference << std::endl + << query << std::endl; + return false; + } + + alt_it += op.second; + } break; + case DELETION: { + if (i && cigar_[i - 1].first == INSERTION) { + std::cerr << "DELETION after INSERTION" << std::endl + << to_string() << std::endl + << reference << std::endl + << query << std::endl; + return false; + } + + if (ref_it > reference.end() - op.second) { + std::cerr << "Reference too short after " + << Cigar::opt_to_char(op.first) << std::endl + << to_string() << std::endl + << reference << std::endl + << query << std::endl; + return false; + } + + ref_it += op.second; + } break; + case NODE_INSERTION: { + // do nothing + } break; + } + } + + if (ref_it != reference.end()) { + std::cerr << "Reference end not reached" << std::endl + << to_string() << std::endl + << reference << std::endl + << query << std::endl; + return false; + } + + if (alt_it != query.end()) { + std::cerr << "Query end not reached" << std::endl + << to_string() << std::endl + << reference << std::endl + << query << std::endl; + return false; + } + + return true; +} + +} // namespace mtg::graph::align diff --git a/metagraph/src/graph/alignment_redone/aln_cigar.hpp b/metagraph/src/graph/alignment_redone/aln_cigar.hpp new file mode 100644 index 0000000000..2c3a3941ae --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_cigar.hpp @@ -0,0 +1,114 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +namespace mtg::graph::align { + +class Cigar { + public: + enum Operator : int8_t { + CLIPPED, + MISMATCH, + MATCH, + DELETION, + INSERTION, + NODE_INSERTION + }; + + typedef uint32_t LengthType; + typedef std::pair value_type; + + Cigar(Operator op = CLIPPED, LengthType num = 0) + : cigar_(num ? 1 : 0, std::make_pair(op, num)) { } + + // See section 1.4 in https://samtools.github.io/hts-specs/SAMv1.pdf for + // a specification of the CIGAR string format. + // e.g., 3=1X2I3D for 3 matches, 1 mismatch, 2 insertions, 3 deletions + // The symbol 'G' is introduced to indicate the insertion of a graph node. + Cigar(std::string_view cigar_str); + + size_t size() const { return cigar_.size(); } + bool empty() const { return cigar_.empty(); } + + std::string to_string() const; + + void append(Operator op, LengthType num = 1); + void append(Cigar&& other); + + LengthType trim_clipping() { + if (cigar_.size() && cigar_.front().first == CLIPPED) { + LengthType ret_val = cigar_.front().second; + cigar_.erase(cigar_.begin(), cigar_.begin() + 1); + return ret_val; + } else { + return 0; + } + } + + LengthType trim_end_clipping() { + if (cigar_.size() && cigar_.back().first == CLIPPED) { + LengthType ret_val = cigar_.back().second; + cigar_.pop_back(); + return ret_val; + } else { + return 0; + } + } + + LengthType get_clipping() const { + return cigar_.size() && cigar_.front().first == CLIPPED ? cigar_.front().second : 0; + } + + LengthType get_end_clipping() const { + return cigar_.size() && cigar_.back().first == CLIPPED ? cigar_.back().second : 0; + } + + void extend_clipping(LengthType n) { + assert(cigar_.size()); + if (cigar_.front().first != CLIPPED) { + cigar_.insert(cigar_.begin(), value_type(CLIPPED, n)); + } else { + cigar_.front().second += n; + } + } + + std::vector& data() { return cigar_; } + const std::vector& data() const { return cigar_; } + + bool operator==(const Cigar &other) const { return cigar_ == other.cigar_; } + bool operator!=(const Cigar &other) const { return !(*this == other); } + + size_t get_num_matches() const { + return std::accumulate(cigar_.begin(), cigar_.end(), 0, + [&](size_t old, const value_type &op) { + return old + (op.first == MATCH) * op.second; + }); + } + + // Return true if the cigar is valid + bool is_valid(std::string_view reference, std::string_view query) const; + + bool is_exact_match(size_t query_size) const { + return cigar_.size() == 1 && cigar_.front() == value_type{ MATCH, query_size }; + } + + static constexpr char opt_to_char(Cigar::Operator op) { return op_str_[op]; } + + private: + static constexpr char op_str_[] = "SX=DIG"; + std::vector cigar_; +}; + + +typedef std::array, 128> OperatorTable; +OperatorTable initialize_opt_table(); +static const OperatorTable kCharToOp = initialize_opt_table(); + + +} // namespace mtg::graph::align diff --git a/metagraph/src/graph/alignment_redone/aln_match.cpp b/metagraph/src/graph/alignment_redone/aln_match.cpp new file mode 100644 index 0000000000..d5ef048dfc --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_match.cpp @@ -0,0 +1,123 @@ +#include "aln_match.hpp" + +#include "common/logger.hpp" + +namespace mtg::graph::align { + +using mtg::common::logger; + +// Return the string spelled by the path. This path may have disconnects (if it came) +// from a chain alignment), so this method handles that case. If there is an invalid +// edge, or if there is too long of a stretch of npos nodes, this throws a runtime error. +std::string spell_path(const DeBruijnGraph &graph, + const std::vector &path) { + std::string seq; + + if (path.empty()) + return seq; + + seq.reserve(path.size() + graph.get_k() - 1); + + size_t num_dummy = 0; + size_t num_unknown = 0; + if (path.front()) { + seq += graph.get_node_sequence(path.front()); + } else { + num_unknown = graph.get_k(); + seq += std::string(num_unknown, '$'); + num_dummy = 1; + } + + for (size_t i = 1; i < path.size(); ++i) { + if (path[i]) { + if (num_dummy) { + seq += '$'; + ++num_unknown; + std::string next_seq = graph.get_node_sequence(path[i]); + auto it = seq.end() - next_seq.size(); + for (char c : next_seq) { + if (*it == '$' && c != '$') { + --num_unknown; + *it = c; + } + + ++it; + } + num_dummy = 0; + } else { + char next = '\0'; + graph.call_outgoing_kmers(path[i - 1], [&](auto next_node, char c) { + if (next_node == path[i]) + next = c; + }); + + if (!next) { + logger->error("Invalid edge {} {}\t{} {}", + path[i - 1], path[i], + graph.get_node_sequence(path[i - 1]), + graph.get_node_sequence(path[i])); + throw std::runtime_error(""); + } + + seq += next; + if (num_unknown) { + std::string next_seq = graph.get_node_sequence(path[i]); + auto it = seq.end() - next_seq.size(); + for (char c : next_seq) { + if (*it == '$' && c != '$') { + --num_unknown; + *it = c; + } + + ++it; + } + } + } + } else { + seq += '$'; + ++num_dummy; + ++num_unknown; + } + } + + assert(seq.size() == path.size() + graph.get_k() - 1); + return seq; +} + +void Anchor::append(const Anchor &other, const DeBruijnGraph *graph) { + assert(query_.data() == other.query_.data()); + assert(query_.data() + query_.size() == other.query_.data() + other.query_.size()); + assert(label_id_ == other.label_id_); + + if (empty()) { + *this = other; + return; + } + + assert(other.seed_.data() >= seed_.data()); + + if (seed_.data() + seed_.size() >= other.seed_.data() + other.seed_.size()) + return; + + if (seed_.data() == other.seed_.data()) { + assert(coord_ == other.coord_); + *this = other; + return; + } + + auto next = path_.begin() + (other.seed_.data() - seed_.data()); + assert(next <= path_.end()); + + seed_ = std::string_view( + seed_.data(), + other.seed_.data() + other.seed_.size() - seed_.data() + ); + + path_.erase(next, path_.end()); + path_.insert(path_.end(), other.path_.begin(), other.path_.end()); + suffix_ = other.suffix_; + + assert(is_valid()); +} + +} // namespace mtg::graph::align \ No newline at end of file diff --git a/metagraph/src/graph/alignment_redone/aln_match.hpp b/metagraph/src/graph/alignment_redone/aln_match.hpp new file mode 100644 index 0000000000..4120ab51b1 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_match.hpp @@ -0,0 +1,152 @@ +#pragma once + +#include + +#include "aln_cigar.hpp" +#include "aligner_config.hpp" + +#include "graph/annotated_dbg.hpp" +#include "annotation/int_matrix/base/int_matrix.hpp" + +namespace mtg::graph::align { + +std::string spell_path(const DeBruijnGraph &graph, + const std::vector &path); + +class Match { + public: + using node_index = SequenceGraph::node_index; + using score_t = DBGAlignerConfig::score_t; + + static constexpr score_t ninf = DBGAlignerConfig::ninf; + + virtual ~Match() {} + + std::string_view get_query() const { return query_; } + std::string_view get_seed() const { return seed_; } + + const std::vector& get_path() const { return path_; } + + virtual std::string_view get_spelling() const = 0; + virtual std::string get_path_spelling() const = 0; + virtual std::string_view get_trim_spelling() const = 0; + + score_t get_score() const { return score_; } + bool get_orientation() const { return orientation_; } + + virtual size_t get_end_trim() const = 0; + + bool empty() const { return path_.empty(); } + + size_t get_clipping() const { return seed_.data() - query_.data(); } + size_t get_end_clipping() const { + return query_.data() + query_.size() - (seed_.data() + seed_.size()); + } + + bool is_spelling_valid(const DeBruijnGraph &graph) const { + return get_path_spelling() == spell_path(graph, path_); + } + + protected: + Match(std::string_view query, + size_t begin, + size_t end, + bool orientation, + std::vector&& path, + score_t score) + : query_(query), orientation_(orientation), seed_(query_.substr(begin, end - begin)), + path_(std::move(path)), score_(score) { + assert(is_spelling_valid()); + } + + std::string_view query_; + bool orientation_; + std::string_view seed_; + std::vector path_; + score_t score_; +}; + +class Anchor : public Match { + public: + using label_class_t = uint64_t; + static constexpr label_class_t nlabel = std::numeric_limits::max(); + static constexpr int64_t ncoord = std::numeric_limits::max(); + + Anchor(Anchor&&) = default; + Anchor(const Anchor&) = default; + Anchor& operator=(const Anchor&) = default; + Anchor& operator=(Anchor&&) = default; + + Anchor(std::string_view query, + size_t begin, + size_t end, + bool orientation, + std::vector&& path, + const DBGAlignerConfig &config, + std::string&& suffix = "", + label_class_t label_class = nlabel, + int64_t coord = ncoord) + : Match(query, begin, end, orientation, + std::move(path), config.match_score(query)), + label_class_(label_class), coord_(coord), suffix_(std::move(suffix)) {} + + std::string_view get_spelling() const override final { return seed_; } + std::string get_path_spelling() const override final { return std::string(seed_) + suffix_; } + std::string_view get_trim_spelling() const override final { return suffix_; } + + size_t get_end_trim() const override final { return suffix_.size(); } + + void append(const Anchor &other, const DeBruijnGraph *graph = nullptr); + + label_class_t get_label_class() const { return label_class_; } + + private: + label_class_t label_class_; + int64_t coord_; + std::string suffix_; +}; + +class Alignment : public Match { + public: + Alignment(const Alignment&) = default; + Alignment(Alignment&&) = default; + Alignment& operator=(const Alignment&) = default; + Alignment& operator=(Alignment&&) = default; + + Alignment(const Anchor &anchor) + : Match(anchor.get_query(), anchor.get_seed().data() - anchor.get_query().data(), + anchor.get_query().data() + anchor.get_query().size() - (anchor.get_seed().data() + anchor.get_seed().size()), + anchor.get_orientation(), + std::vector(anchor.get_path()), + anchor.get_score()), + end_trim_(anchor.get_end_trim()), + path_spelling_(anchor.get_path_spelling()), + cigar_(std::string(get_clipping(), 'S') + + std::string(seed_.size(), '=') + + std::string(get_end_clipping(), 'S')), + label_classes_(path_.size(), anchor.get_label_class()) { + assert(path_spelling_.size() == query_.size() + end_trim_); + } + + const Cigar& get_cigar() const { return cigar_; } + + size_t get_end_trim() const override final { return end_trim_; } + + std::string_view get_spelling() const override final { + return std::string_view(path_spelling_.c_str(), path_spelling_.size() - end_trim_); + } + + std::string get_path_spelling() const override final { return path_spelling_; } + std::string_view get_trim_spelling() const override final { + return std::string_view(path_spelling_.c_str() + path_spelling_.size() - end_trim_, + end_trim_); + } + + private: + size_t end_trim_; + std::string path_spelling_; + Cigar cigar_; + std::vector label_classes_; +}; + +} // namespace mtg::graph::align \ No newline at end of file diff --git a/metagraph/src/graph/alignment_redone/aln_query.hpp b/metagraph/src/graph/alignment_redone/aln_query.hpp new file mode 100644 index 0000000000..d029099007 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_query.hpp @@ -0,0 +1,34 @@ +#pragma once + +#include + +#include "graph/representation/base/sequence_graph.hpp" +#include "common/seq_tools/reverse_complement.hpp" + +namespace mtg::graph::align { + +class Query { + public: + template + Query(const DeBruijnGraph &graph, Args&&... args) + : graph_(graph), fwd_(std::forward(args)...), rc_(generate_rc()) {} + + const std::string& get_query(bool orientation = false) const { + return !orientation ? fwd_ : rc_; + } + + const DeBruijnGraph& get_graph() const { return graph_; } + + private: + const DeBruijnGraph &graph_; + std::string fwd_; + std::string rc_; + + std::string generate_rc() { + std::string rc = fwd_; + reverse_complement(rc_.begin(), rc_.end()); + return rc; + } +}; + +} // namespace mtg::graph::align \ No newline at end of file diff --git a/metagraph/src/graph/alignment_redone/aln_seeder.cpp b/metagraph/src/graph/alignment_redone/aln_seeder.cpp new file mode 100644 index 0000000000..08870caacd --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_seeder.cpp @@ -0,0 +1,434 @@ +#include "aln_seeder.hpp" + +#include "aln_chainer.hpp" +#include "common/logger.hpp" +#include "common/vector_map.hpp" + +namespace mtg::graph::align { + +std::vector Seeder::get_inexact_anchors() const { + auto anchors = get_anchors(); + std::vector inexact_anchors; + inexact_anchors.reserve(anchors.size()); + std::transform(anchors.begin(), anchors.end(), std::back_inserter(inexact_anchors), + [&](const Anchor &a) { return Alignment(a); }); + return inexact_anchors; +} + +std::vector ExactSeeder::get_anchors() const { + const DeBruijnGraph &graph = query_.get_graph(); + std::vector anchors; + + if (query_.get_query().size() < config_.min_seed_length) + return anchors; + + if (config_.min_seed_length == graph.get_k()) { + for (bool orientation : { false, true }) { + size_t begin = 0; + graph.map_to_nodes_sequentially(query_.get_query(orientation), + [&](Match::node_index node) { + if (node != DeBruijnGraph::npos) { + anchors.emplace_back(query_.get_query(), + begin, begin + graph.get_k(), + orientation, + std::vector{ node }, + config_); + } + + ++begin; + } + ); + } + + return anchors; + } + + return anchors; +} + +void call_distances(const DeBruijnGraph &graph, + DeBruijnGraph::node_index start, + DeBruijnGraph::node_index target, + const std::function &callback, + const std::function &terminate_branch = [](ssize_t) { return false; }, + const std::function &terminate = []() { return false; }) { + // TODO: this is a slow DFS-based method, replace with something smarter + std::vector> traversal_stack; + traversal_stack.emplace_back(start, 0); + while (traversal_stack.size() && !terminate()) { + auto [node, coord_dist] = traversal_stack.back(); + traversal_stack.pop_back(); + + if (node == target) + callback(coord_dist); + + if (terminate_branch(coord_dist)) + break; + + ++coord_dist; + graph.adjacent_outgoing_nodes(node, [&](DeBruijnGraph::node_index next) { + traversal_stack.emplace_back(next, coord_dist); + }); + } +} + +void global_align(const DeBruijnGraph &graph, + const DBGAlignerConfig &config, + const Alignment &aln, + const Anchor &target, + size_t dist, + const AlignmentCallback &callback) { + assert(dist > 0); + + std::string_view query_window = aln.get_query(); + query_window.remove_suffix(aln.get_end_clipping() + aln.get_seed().size()); + query_window.remove_prefix(target.get_clipping()); + if (query_window.empty()) + return; + + // score = 1/2(match_score(query_size + match_spelling_size - cost)) + // => 2*score/match_score - query_size - match_spelling_size = -cost + // => cost = query_size + match_spelling_size - 2*score/match_score + + DBGAlignerConfig::score_t match_score = config.match_score("A"); + auto score_to_cost = [&,match_score](DBGAlignerConfig::score_t score, ssize_t query=0, ssize_t match=0) -> ssize_t { + query = query_window.size(); + match = dist; + return query + match - 2 * (score - aln.get_score()) / match_score; + }; + + auto cost_to_score = [&,match_score](ssize_t cost, ssize_t query=0, ssize_t match=0) -> DBGAlignerConfig::score_t { + query = query_window.size(); + match = dist; + return match_score * (query + match - cost) / 2 + aln.get_score(); + }; + + size_t mismatch_cost = score_to_cost(-match_score); + size_t gap_opn = score_to_cost(config.gap_opening_penalty); + size_t gap_ext = score_to_cost(config.gap_extension_penalty); + + size_t query_dist = 0; + auto node = aln.get_path()[0]; + size_t best_dist = 0; + for (auto it = query_window.rbegin(); it != query_window.rend(); ++it) { + if (best_dist == dist || !graph.has_single_incoming(node)) + break; + + if (auto prev = graph.traverse_back(node, *it)) { + ++best_dist; + ++query_dist; + node = prev; + } else { + break; + } + } + + // S(cost, query_dist, node) = best_dist + using NodeToBestDists = VectorMap; // S(cost, query_dist) + using WaveFront = std::vector; // S(cost) + using ScoreTable = std::vector; + + ScoreTable S; // best cost + ScoreTable E; // best cost (last operation is insertion) + ScoreTable F; // best cost (last operation is deletion) + + auto fill_table = [&]() -> size_t { + auto &wf = S.emplace_back(); + wf.resize(query_dist + 1); + wf[query_dist][node] = best_dist; + if (best_dist == dist) { + if (node == target.get_path().back()) + return 0; + + return std::numeric_limits::max(); + } + + for (size_t cost = 0; cost < S.size(); ++cost) { + if (S[cost].empty()) + continue; + + for (size_t query_dist = 0; query_dist < S[cost].size(); ++query_dist) { + if (S[cost][query_dist].empty()) + continue; + + size_t it_dist = 0; + for (auto it = S[cost][query_dist].begin(); it != S[cost][query_dist].end(); ++it, ++it_dist) { + assert(query_dist <= query_window.size()); + + auto set_value = [&](auto &table, size_t cost, size_t query_dist, auto node, size_t dist) -> size_t { + if (cost >= table.size()) { + table.resize(cost + 1); + if (&table == &S) + it = S[cost][query_dist].begin() + it_dist; + } + + if (query_dist >= table[cost].size()) { + table[cost].resize(query_dist + 1); + if (&table == &S) + it = S[cost][query_dist].begin() + it_dist; + } + + auto &bucket = table[cost][query_dist][node]; + if (&table == &S) + it = S[cost][query_dist].begin() + it_dist; + + bucket = std::max(bucket, dist); + + return bucket; + }; + + DeBruijnGraph::node_index node = it->first; + size_t best_dist = it->second; + + if (best_dist == dist) { + if (node == target.get_path().back() && query_dist == query_window.size()) + return cost; + + continue; + } + + // DBGAlignerConfig::score_t prev_score = cost_to_score(cost, query_dist, best_dist); + + // forward creation of insertions + if (query_dist + 1 <= query_window.size()) { + // extend a previous insertion + if (cost < E.size() && E[cost].size() && query_dist < E[cost].size()) { + auto find = E[cost][query_dist].find(node); + if (find != E[cost][query_dist].end()) { + size_t last_dist = find->second; + ssize_t next_ext_cost = cost + gap_ext; + // ssize_t next_ext_cost = score_to_cost(prev_score + config.gap_extension_penalty, query_dist + 1, last_dist); + size_t stored_dist = set_value(E, next_ext_cost, query_dist + 1, node, last_dist); + set_value(S, next_ext_cost, query_dist + 1, node, stored_dist); + } + } + + // open an insertion + ssize_t next_opn_cost = cost + gap_opn; + // ssize_t next_opn_cost = score_to_cost(prev_score + config.gap_opening_penalty, query_dist + 1, best_dist); + size_t stored_dist = set_value(E, next_opn_cost, query_dist + 1, node, best_dist); + set_value(S, next_opn_cost, query_dist + 1, node, stored_dist); + } + + if (best_dist < dist) { + if (query_dist < query_window.size()) { + // match + std::string_view local_query_window = query_window; + local_query_window.remove_suffix(query_dist); + + // const auto &score_row = config.score_matrix[local_query_window.back()]; + graph.call_incoming_kmers(node, [&](DeBruijnGraph::node_index prev, char c) { + size_t cur_best = best_dist + 1; + size_t cur_query_dist = query_dist + 1; + size_t cur_cost = c == local_query_window.back() + ? cost + : cost + mismatch_cost; + // : score_to_cost(prev_score + score_row[c], cur_query_dist, cur_best); + for (auto jt = local_query_window.rbegin() + 1; jt != local_query_window.rend(); ++jt) { + if (cur_best == dist || !graph.has_single_incoming(prev)) + break; + + if (auto pprev = graph.traverse_back(prev, *jt)) { + ++cur_best; + ++cur_query_dist; + prev = pprev; + } else { + break; + } + } + + assert(cur_query_dist > query_dist || cur_cost > cost); + set_value(S, cur_cost, cur_query_dist, prev, cur_best); + }); + + } else { + // deletion + std::vector prevs; + graph.adjacent_incoming_nodes(node, [&](auto prev) { + prevs.emplace_back(prev); + }); + if (prevs.size()) { + if (cost < F.size() && F[cost].size() && query_dist < F[cost].size()) { + // extension + auto find = F[cost][query_dist].find(node); + if (find != F[cost][query_dist].end()) { + size_t last_dist = find->second; + ssize_t next_ext_cost = cost + gap_ext; + // ssize_t next_ext_cost = score_to_cost(prev_score + config.gap_extension_penalty, query_dist, last_dist + 1); + for (auto prev : prevs) { + size_t stored_value = set_value(F, next_ext_cost, query_dist, prev, last_dist + 1); + set_value(S, next_ext_cost, query_dist, prev, stored_value); + } + } + } + + ssize_t next_opn_cost = cost + gap_opn; + // ssize_t next_opn_cost = score_to_cost(prev_score + config.gap_opening_penalty, query_dist, best_dist + 1); + for (auto prev : prevs) { + size_t stored_value = set_value(F, next_opn_cost, query_dist, prev, best_dist + 1); + set_value(S, next_opn_cost, query_dist, prev, stored_value); + } + } + } + } + } + } + } + + return std::numeric_limits::max(); + }; + + size_t cost = fill_table(); + if (cost < std::numeric_limits::max()) { + common::logger->trace("FOUND ALIGNMENT AT COST {}", cost); + } else { + throw std::runtime_error("No alignment found"); + } +} + +std::vector ExactSeeder::get_inexact_anchors() const { + std::vector anchors = get_anchors(); + std::vector alignments; + + using AnchorIt = std::vector::iterator; + + chain_anchors(query_, config_, anchors.begin(), anchors.end(), + [this](const Anchor &a_j, + ssize_t max_dist, + AnchorIt begin, + AnchorIt end, + typename ChainScores::iterator chain_scores, + const ScoreUpdater &update_score) { + assert(a_i.get_label_class() == a_j.get_label_class()); + assert(a_i.get_orientation() == a_j.get_orientation()); + + const DeBruijnGraph &graph = query_.get_graph(); + std::string_view query_j = a_j.get_seed(); + const DBGAlignerConfig::score_t &score_j = std::get<0>(*(chain_scores + (end - begin))); + + for (auto it = begin; it != end; ++it) { + const Anchor &a_i = *it; + std::string_view query_i = a_i.get_seed(); + + DBGAlignerConfig::score_t dist = query_j.end() - query_i.end(); + + if (dist <= 0 || query_i.begin() >= query_j.begin()) { + ++chain_scores; + return; + } + + DBGAlignerConfig::score_t base_score = std::get<0>(*chain_scores); + + if (query_i.end() <= query_j.begin()) { + base_score += a_j.get_score(); + } else { + std::string_view ext(query_i.data() + query_i.size(), + query_j.end() - query_i.end()); + base_score += config_.match_score(ext); + + if (!a_j.get_end_clipping()) + base_score += config_.right_end_bonus; + } + + if (base_score <= score_j) { + ++chain_scores; + return; + } + + DBGAlignerConfig::score_t score = base_score; + DBGAlignerConfig::score_t coord_dist = dist; + size_t front_dist = a_j.get_clipping() - a_i.get_clipping(); + if (front_dist < a_i.get_path().size()) { + if (a_i.get_path()[front_dist] != a_j.get_path()[0]) { + ++chain_scores; + return; + } + } else if (front_dist < a_i.get_path().size() + graph.get_k() - 1) { + std::string path_spelling_i = a_i.get_path_spelling(); + std::string path_spelling_j = a_j.get_path_spelling(); + auto [i_mismatch, j_mismatch] = std::mismatch( + path_spelling_i.begin() + front_dist, path_spelling_i.end(), + path_spelling_j.c_str(), path_spelling_j.c_str() + path_spelling_j.size() + ); + + if (i_mismatch != path_spelling_i.end()) { + ++chain_scores; + return; + } + + if (path_spelling_i.size() - front_dist < graph.get_k() - 1) { + // not enough overlap, we need to see if a path exists + auto j_check_end = path_spelling_j.c_str() + graph.get_k(); + assert(j_mismatch < j_check_end); + Match::node_index last_node = DeBruijnGraph::npos; + graph.traverse(a_i.get_path().back(), + j_mismatch, j_check_end, + [&](Match::node_index node) { + last_node = node; + ++j_mismatch; + }); + if (j_mismatch != j_check_end || last_node != a_j.get_path()[0]) { + ++chain_scores; + return; + } + + assert(last_node != DeBruijnGraph::npos); + } + } else { + // completely disjoint, we need to explore freely + coord_dist = -DBGAlignerConfig::ninf; + ssize_t min_diff = std::numeric_limits::max(); + call_distances(graph, a_i.get_path().back(), a_j.get_path()[0], + [&](ssize_t path_dist) { + DBGAlignerConfig::score_t cur_coord_dist + = a_i.get_end_trim() + path_dist + a_j.get_path().size() - 1 - a_j.get_end_trim(); + ssize_t cur_diff = abs(cur_coord_dist - dist); + if (cur_diff < min_diff) { + min_diff = cur_diff; + coord_dist = cur_coord_dist; + } + }, + [&](ssize_t cur_dist) { + if (cur_dist > dist) { + return cur_dist - dist > std::min(min_diff, max_dist); + } else { + return dist - cur_dist > max_dist; + } + }, + [&]() { return min_diff == 0; } + ); + + if (coord_dist == -DBGAlignerConfig::ninf) { + ++chain_scores; + return; + } + } + + DBGAlignerConfig::score_t diff = std::abs(coord_dist - dist); + score += config_.gap_opening_penalty + + (diff - 1) * config_.gap_extension_penalty; + + update_score(score, it, coord_dist); + + ++chain_scores; + } + }, + [&](const AnchorChain &chain, const std::vector &score_traceback) { + return true; + }, + [this](AnchorIt next, + AnchorIt last, + Alignment&& aln, + size_t last_to_next_dist, + DBGAlignerConfig::score_t score_up_to_now, + const AlignmentCallback &callback) { + size_t next_to_last_dist = last_to_next_dist - next->get_seed().size() + last->get_seed().size() - last->get_path().size() + 1; + global_align(query_.get_graph(), config_, aln, *next, next_to_last_dist, callback); + }, + [&alignments](Alignment&& alignment) { alignments.emplace_back(std::move(alignment)); } + ); + + return alignments; +} + +} // namespace mtg::graph::align \ No newline at end of file diff --git a/metagraph/src/graph/alignment_redone/aln_seeder.hpp b/metagraph/src/graph/alignment_redone/aln_seeder.hpp new file mode 100644 index 0000000000..9719f4030e --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_seeder.hpp @@ -0,0 +1,32 @@ +#pragma once + +#include "aln_match.hpp" +#include "aln_query.hpp" +#include "aligner_config.hpp" + +namespace mtg::graph::align { + +class Seeder { + public: + Seeder(const Query &query, const DBGAlignerConfig &config) + : query_(query), config_(config) {} + + virtual ~Seeder() {} + + virtual std::vector get_anchors() const = 0; + virtual std::vector get_inexact_anchors() const; + + const Query& get_query() const { return query_; } + + protected: + const Query &query_; + const DBGAlignerConfig &config_; +}; + +class ExactSeeder : public Seeder { + public: + std::vector get_anchors() const override; + std::vector get_inexact_anchors() const override; +}; + +} // namespace mtg::graph::align \ No newline at end of file diff --git a/metagraph/src/graph/alignment_redone/annotation_buffer.cpp b/metagraph/src/graph/alignment_redone/annotation_buffer.cpp new file mode 100644 index 0000000000..4020f312a7 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/annotation_buffer.cpp @@ -0,0 +1,221 @@ +#include "annotation_buffer.hpp" + +#include "graph/representation/rc_dbg.hpp" +#include "graph/representation/succinct/dbg_succinct.hpp" +#include "graph/representation/canonical_dbg.hpp" +#include "annotation/binary_matrix/base/binary_matrix.hpp" +#include "common/utils/template_utils.hpp" + +namespace mtg { +namespace graph { +namespace align { + +using mtg::common::logger; + +typedef annot::matrix::BinaryMatrix::Row Row; +typedef annot::matrix::BinaryMatrix::Column Column; + +// dummy index for an unfetched annotations +static constexpr size_t nannot = std::numeric_limits::max(); + +AnnotationBuffer::AnnotationBuffer(const DeBruijnGraph &graph, const Annotator &annotator) + : graph_(graph), + annotator_(annotator), + multi_int_(dynamic_cast(&annotator_.get_matrix())), + canonical_(dynamic_cast(&graph_)), + column_sets_({ {} }) { + if (multi_int_ && graph_.get_mode() != DeBruijnGraph::BASIC) { + multi_int_ = nullptr; + logger->warn("Coordinates not supported when aligning to CANONICAL " + "or PRIMARY mode graphs"); + } +} + +void AnnotationBuffer::fetch_queued_annotations() { + assert(graph_.get_mode() != DeBruijnGraph::PRIMARY + && "PRIMARY graphs must be wrapped into CANONICAL"); + + std::vector queued_nodes; + std::vector queued_rows; + + const DeBruijnGraph *base_graph = &graph_; + + if (canonical_) + base_graph = &canonical_->get_graph(); + + const auto *dbg_succ = dynamic_cast(base_graph); + const boss::BOSS *boss = dbg_succ ? &dbg_succ->get_boss() : nullptr; + + for (const auto &path : queued_paths_) { + std::vector base_path; + if (base_graph->get_mode() == DeBruijnGraph::CANONICAL) { + // TODO: avoid this call of spell_path + std::string query = spell_path(graph_, path); + base_path = map_to_nodes(*base_graph, query); + + } else if (canonical_) { + base_path.reserve(path.size()); + for (node_index node : path) { + base_path.emplace_back(canonical_->get_base_node(node)); + } + + } else { + assert(graph_.get_mode() == DeBruijnGraph::BASIC); + base_path = path; + if (dynamic_cast(&graph_)) + std::reverse(base_path.begin(), base_path.end()); + } + + assert(base_path.size() == path.size()); + + for (size_t i = 0; i < path.size(); ++i) { + if (base_path[i] == DeBruijnGraph::npos) { + // this can happen when the base graph is CANONICAL and path[i] is a + // dummy node + if (node_to_cols_.try_emplace(path[i], 0).second && has_coordinates()) + label_coords_.emplace_back(); + + continue; + } + + if (boss && !boss->get_W(dbg_succ->kmer_to_boss_index(base_path[i]))) { + // skip dummy nodes + if (node_to_cols_.try_emplace(base_path[i], 0).second && has_coordinates()) + label_coords_.emplace_back(); + + if (graph_.get_mode() == DeBruijnGraph::CANONICAL + && base_path[i] != path[i] + && node_to_cols_.emplace(path[i], 0).second && has_coordinates()) { + label_coords_.emplace_back(); + } + + continue; + } + + Row row = AnnotatedDBG::graph_to_anno_index(base_path[i]); + if (canonical_ || graph_.get_mode() == DeBruijnGraph::BASIC) { + if (node_to_cols_.try_emplace(base_path[i], nannot).second) { + queued_rows.push_back(row); + queued_nodes.push_back(base_path[i]); + } + + continue; + } + + assert(graph_.get_mode() == DeBruijnGraph::CANONICAL); + + auto find_a = node_to_cols_.find(path[i]); + auto find_b = node_to_cols_.find(base_path[i]); + + if (find_a == node_to_cols_.end() && find_b == node_to_cols_.end()) { + node_to_cols_.try_emplace(path[i], nannot); + queued_rows.push_back(row); + queued_nodes.push_back(path[i]); + + if (path[i] != base_path[i]) { + node_to_cols_.emplace(base_path[i], nannot); + queued_rows.push_back(row); + queued_nodes.push_back(base_path[i]); + } + } else if (find_a == node_to_cols_.end() && find_b != node_to_cols_.end()) { + node_to_cols_.try_emplace(path[i], find_b->second); + if (find_b->second == nannot) { + queued_rows.push_back(row); + queued_nodes.push_back(path[i]); + } + } else if (find_a != node_to_cols_.end() && find_b == node_to_cols_.end()) { + node_to_cols_.try_emplace(base_path[i], find_a->second); + } else { + size_t label_i = std::min(find_a->second, find_b->second); + if (label_i != nannot) { + find_a.value() = label_i; + find_b.value() = label_i; + } + } + } + } + + queued_paths_.clear(); + + if (queued_nodes.empty()) + return; + + auto push_node_labels = [&](auto node_it, auto row_it, auto&& labels) { + assert(node_it != queued_nodes.end()); + assert(node_to_cols_.count(*node_it)); + assert(node_to_cols_.count(AnnotatedDBG::anno_to_graph_index(*row_it))); + + size_t label_i = cache_column_set(std::move(labels)); + node_index base_node = AnnotatedDBG::anno_to_graph_index(*row_it); + if (graph_.get_mode() == DeBruijnGraph::BASIC) { + assert(base_node == *node_it); + node_to_cols_[*node_it] = label_i; + } else if (canonical_) { + node_to_cols_[base_node] = label_i; + } else { + node_to_cols_[*node_it] = label_i; + if (base_node != *node_it && node_to_cols_.try_emplace(base_node, label_i).second + && has_coordinates()) { + label_coords_.emplace_back(label_coords_.back()); + } + } + }; + + auto node_it = queued_nodes.begin(); + auto row_it = queued_rows.begin(); + if (has_coordinates()) { + assert(multi_int_); + // extract both labels and coordinates, then store them separately + for (auto&& row_tuples : multi_int_->get_row_tuples(queued_rows)) { + std::sort(row_tuples.begin(), row_tuples.end(), utils::LessFirst()); + Columns labels; + labels.reserve(row_tuples.size()); + label_coords_.emplace_back(); + label_coords_.back().reserve(row_tuples.size()); + for (auto&& [label, coords] : row_tuples) { + labels.push_back(label); + label_coords_.back().emplace_back(coords.begin(), coords.end()); + } + push_node_labels(node_it++, row_it++, std::move(labels)); + } + } else { + for (auto&& labels : annotator_.get_matrix().get_rows(queued_rows)) { + std::sort(labels.begin(), labels.end()); + push_node_labels(node_it++, row_it++, std::move(labels)); + } + } + +#ifndef NDEBUG + for (const auto &[node, val] : node_to_cols_) { + assert(val != nannot); + } +#endif +} + +auto AnnotationBuffer::get_labels_and_coords(node_index node) const + -> std::pair { + std::pair ret_val { nullptr, nullptr }; + + if (canonical_) + node = canonical_->get_base_node(node); + + auto it = node_to_cols_.find(node); + + // if the node hasn't been seen before, or if its annotations haven't + // been fetched, return nothing + if (it == node_to_cols_.end() || it->second == nannot) + return ret_val; + + ret_val.first = &column_sets_.data()[it->second]; + + if (has_coordinates()) { + assert(static_cast(it - node_to_cols_.begin()) < label_coords_.size()); + ret_val.second = &label_coords_[it - node_to_cols_.begin()]; + } + + return ret_val; +} + +} // namespace align +} // namespace graph +} // namespace mtg diff --git a/metagraph/src/graph/alignment_redone/annotation_buffer.hpp b/metagraph/src/graph/alignment_redone/annotation_buffer.hpp new file mode 100644 index 0000000000..128d56797f --- /dev/null +++ b/metagraph/src/graph/alignment_redone/annotation_buffer.hpp @@ -0,0 +1,85 @@ +#pragma once + +#include "aln_match.hpp" +#include "graph/annotated_dbg.hpp" +#include "annotation/int_matrix/base/int_matrix.hpp" +#include "common/vector_set.hpp" +#include "common/vector_map.hpp" +#include "common/hashers/hash.hpp" + +namespace mtg::graph { + +class CanonicalDBG; + +namespace align { + +// caches queried annotations to speed up next queries (labels with or w/o coordinates) +class AnnotationBuffer { + public: + typedef AnnotatedDBG::Annotator Annotator; + typedef DeBruijnGraph::node_index node_index; + using Column = annot::matrix::BinaryMatrix::Column; + using Tuple = SmallVector; + using Columns = Vector; + using CoordinateSet = Vector; + + AnnotationBuffer(const DeBruijnGraph &graph, const Annotator &annotator); + + void queue_path(std::vector&& path) { + queued_paths_.push_back(std::move(path)); + } + + // fetch annotations for the queued nodes from the buffer and reset the buffer + void fetch_queued_annotations(); + + bool has_coordinates() const { return multi_int_; } + + // Get the annotations and coordinates of a node if they have been fetched. + // The returned pointers are valid until next fetch_queued_annotations(). + std::pair + get_labels_and_coords(node_index node) const; + + // get the labels of a node if they have been fetched + inline const Columns* get_labels(node_index node) const { + return get_labels_and_coords(node).first; + } + + const Annotator& get_annotator() const { return annotator_; } + + size_t num_nodes_buffered() const { return node_to_cols_.size(); } + size_t num_column_sets() const { return column_sets_.size(); } + + // This method lets the caller push additional column sets to dictionary + // `column_sets_`. These column sets can later be fetched with `get_cached_column_set()` + template + inline size_t cache_column_set(Args&&... args) { + auto it = column_sets_.emplace(std::forward(args)...).first; + assert(std::is_sorted(it->begin(), it->end())); + return it - column_sets_.begin(); + } + + // Fetch a label set given its index returned by `cache_column_set()` + inline const Columns& get_cached_column_set(size_t i) const { + assert(i < column_sets_.size()); + return column_sets_.data()[i]; + } + + private: + const DeBruijnGraph &graph_; + const Annotator &annotator_; + const annot::matrix::MultiIntMatrix *multi_int_; + const CanonicalDBG *canonical_; + + // keep a unique set of annotation rows + // the first element is the empty label set + VectorSet column_sets_; + // map node to index in |column_sets_| + VectorMap node_to_cols_; + // coordinate sets for all nodes in |node_to_cols_| in the same order + std::vector label_coords_; + // buffer of paths to later querying with fetch_queued_annotations() + std::vector> queued_paths_; +}; + +} // namespace align +} // namespace mtg::graph From 18fca444dfb6de55b789c5a9973df596aaf2ea5a Mon Sep 17 00:00:00 2001 From: Harun Mustafa Date: Sun, 13 Oct 2024 18:21:57 +0200 Subject: [PATCH 002/308] update --- .../src/graph/alignment_redone/aln_seeder.cpp | 113 ++++++++++++++---- 1 file changed, 93 insertions(+), 20 deletions(-) diff --git a/metagraph/src/graph/alignment_redone/aln_seeder.cpp b/metagraph/src/graph/alignment_redone/aln_seeder.cpp index 08870caacd..54d454d134 100644 --- a/metagraph/src/graph/alignment_redone/aln_seeder.cpp +++ b/metagraph/src/graph/alignment_redone/aln_seeder.cpp @@ -107,22 +107,6 @@ void global_align(const DeBruijnGraph &graph, size_t gap_opn = score_to_cost(config.gap_opening_penalty); size_t gap_ext = score_to_cost(config.gap_extension_penalty); - size_t query_dist = 0; - auto node = aln.get_path()[0]; - size_t best_dist = 0; - for (auto it = query_window.rbegin(); it != query_window.rend(); ++it) { - if (best_dist == dist || !graph.has_single_incoming(node)) - break; - - if (auto prev = graph.traverse_back(node, *it)) { - ++best_dist; - ++query_dist; - node = prev; - } else { - break; - } - } - // S(cost, query_dist, node) = best_dist using NodeToBestDists = VectorMap; // S(cost, query_dist) using WaveFront = std::vector; // S(cost) @@ -133,6 +117,21 @@ void global_align(const DeBruijnGraph &graph, ScoreTable F; // best cost (last operation is deletion) auto fill_table = [&]() -> size_t { + size_t query_dist = 0; + auto node = aln.get_path()[0]; + size_t best_dist = 0; + for (auto it = query_window.rbegin(); it != query_window.rend(); ++it) { + if (best_dist == dist || !graph.has_single_incoming(node)) + break; + + if (auto prev = graph.traverse_back(node, *it)) { + ++best_dist; + ++query_dist; + node = prev; + } else { + break; + } + } auto &wf = S.emplace_back(); wf.resize(query_dist + 1); wf[query_dist][node] = best_dist; @@ -279,11 +278,85 @@ void global_align(const DeBruijnGraph &graph, }; size_t cost = fill_table(); - if (cost < std::numeric_limits::max()) { - common::logger->trace("FOUND ALIGNMENT AT COST {}", cost); - } else { + if (cost == std::numeric_limits::max()) throw std::runtime_error("No alignment found"); - } + + common::logger->trace("FOUND ALIGNMENT AT COST {}", cost); + + // backtrack + DeBruijnGraph::node_index node = target.get_path().back(); + size_t query_dist = query_window.size(); + + assert(cost < S.size() && query_dist < S[cost].size()); + auto bt = S[cost][query_dist].find(node); + assert(bt != S[cost][query_dist].end()); + assert(bt->second == dist); + +// while (bt != S[cost][query_dist].end()) { +// assert(bt->second == dist); +// if (cost < E.size() && query_dist < E[cost].size()) { +// // check insertion +// auto et = E[cost][query_dist].find(node); +// if (et != E[cost][query_dist].end() && et->second == dist) { +// // S == E, so we have an insertion. check if it's an extension or an insertion + +// // backtrack through all +// bool found = true; +// while (found) { +// found = false; +// size_t prev_cost = cost - gap_ext; +// if (query_dist - 1 < E[prev_cost].size()) { +// auto prev_et = E[prev_cost][query_dist - 1].find(node); +// if (prev_et != E[prev_cost][query_dist - 1].end() && prev_et->second == dist) { +// cost = prev_cost; +// --query_dist; +// et = prev_et; +// found = true; +// } +// } +// } + +// size_t prev_cost = cost - gap_opn; +// if (query_dist - 1 < E[prev_cost].size()) { + +// } else { + +// } +// } +// } + +// if (cost < F.size() && query_dist < F[cost].size()) { +// // check deletion +// auto ft = F[cost][query_dist].find(node); +// if (ft != F[cost][query_dist].end() && ft->second == dist) { +// // S == F, so we may have a deletion +// // TODO +// } +// } + +// if (query_dist > 0 && dist > 0) { +// // check match +// DeBruijnGraph::node_index found_prev = DeBruijnGraph::npos; +// size_t cur_query_dist = query_dist; +// for (auto it = S[cost].rbegin(); it != S[cost].rend(); ++it, --cur_query_dist) { +// for (const auto &[prev_node, prev_dist] : *it) { +// if (query_dist > cur_query_dist && dist > prev_dist && dist - prev_dist == query_dist - cur_query_dist) { +// assert(found_prev == DeBruijnGraph::npos); +// found_prev = prev_node; +// #ifdef NDEBUG +// break; +// #endif +// } +// } +// #ifdef NDEBUG +// if (found_prev != DeBruijnGraph::npos) +// break; +// #endif +// } + +// if (found_prev) +// } +// } } std::vector ExactSeeder::get_inexact_anchors() const { From 60b886855a504e07d61a9066d5bb504bdbf73c6c Mon Sep 17 00:00:00 2001 From: Harun Mustafa Date: Tue, 15 Oct 2024 19:41:24 +0200 Subject: [PATCH 003/308] fixes --- .../graph/alignment_redone/aligner_config.cpp | 9 +- .../graph/alignment_redone/aligner_config.hpp | 2 +- .../graph/alignment_redone/aln_chainer.hpp | 45 ++- .../src/graph/alignment_redone/aln_cigar.cpp | 26 +- .../src/graph/alignment_redone/aln_cigar.hpp | 2 +- .../src/graph/alignment_redone/aln_match.cpp | 6 +- .../src/graph/alignment_redone/aln_match.hpp | 71 +++- .../src/graph/alignment_redone/aln_query.hpp | 2 +- .../src/graph/alignment_redone/aln_seeder.cpp | 335 ++++++++++++------ .../src/graph/alignment_redone/aln_seeder.hpp | 5 +- .../alignment_redone/annotation_buffer.cpp | 2 +- .../alignment_redone/annotation_buffer.hpp | 3 +- 12 files changed, 340 insertions(+), 168 deletions(-) diff --git a/metagraph/src/graph/alignment_redone/aligner_config.cpp b/metagraph/src/graph/alignment_redone/aligner_config.cpp index ec6d04e857..bc0fecf0aa 100644 --- a/metagraph/src/graph/alignment_redone/aligner_config.cpp +++ b/metagraph/src/graph/alignment_redone/aligner_config.cpp @@ -5,7 +5,7 @@ #include "common/logger.hpp" -namespace mtg::graph::align { +namespace mtg::graph::align_redone { using mtg::common::logger; @@ -78,15 +78,16 @@ ::score_cigar(std::string_view reference, auto ref_it = reference.begin(); auto alt_it = query.begin(); auto it = cigar.data().begin(); - if (it->first == Cigar::CLIPPED) + if (it->first == Cigar::CLIPPED) { + alt_it += it->second; ++it; + } for ( ; it != cigar.data().end(); ++it) { const auto &op = *it; switch (op.first) { case Cigar::CLIPPED: { - if (it + 1 != cigar.data().end()) - alt_it += op.second; + alt_it += op.second; } break; case Cigar::MATCH: { score += match_score(std::string_view(ref_it, op.second)); diff --git a/metagraph/src/graph/alignment_redone/aligner_config.hpp b/metagraph/src/graph/alignment_redone/aligner_config.hpp index f372b3bb38..a169f74604 100644 --- a/metagraph/src/graph/alignment_redone/aligner_config.hpp +++ b/metagraph/src/graph/alignment_redone/aligner_config.hpp @@ -8,7 +8,7 @@ #include -namespace mtg::graph::align { +namespace mtg::graph::align_redone { class Cigar; diff --git a/metagraph/src/graph/alignment_redone/aln_chainer.hpp b/metagraph/src/graph/alignment_redone/aln_chainer.hpp index cd22bfb356..f5d1792d6a 100644 --- a/metagraph/src/graph/alignment_redone/aln_chainer.hpp +++ b/metagraph/src/graph/alignment_redone/aln_chainer.hpp @@ -1,10 +1,14 @@ #pragma once +#include + #include "aligner_config.hpp" #include "aln_query.hpp" #include "aln_match.hpp" -namespace mtg::graph::align { +#include "common/logger.hpp" + +namespace mtg::graph::align_redone { template using ChainScores = std::vector>; @@ -50,37 +54,23 @@ void extend_chain(const AnchorChain &chain, const AlignmentCallback &callback, const std::function &terminate = []() { return false; }) { std::vector alns; + assert(chain.size()); alns.emplace_back(*chain.back().first); - for (auto it = chain.rbegin() + 1; it != chain.rend(); ++it) { + for (auto it = chain.rbegin(); it + 1 != chain.rend(); ++it) { std::vector next_alns; for (auto&& aln : alns) { - size_t dist = (it - 1)->second; - anchor_extender((it - 1)->first, it->first, std::move(aln), dist, 0, + std::cerr << aln << " begins with " << *it->first << std::endl; + assert(it->first->get_clipping() == aln.get_clipping()); + + size_t dist = it->second; + anchor_extender(it->first, (it + 1)->first, std::move(aln), dist, 0, [&](Alignment&& next_aln) { next_alns.emplace_back(std::move(next_aln)); } ); } + std::swap(next_alns, alns); } - // auto jt = score_traceback.begin(); - // std::vector alns; - // alns.emplace_back(*chain[0].first); - // ++jt; - - // for (auto it = chain.begin() + 1; it != chain.end(); ++it, ++jt) { - // assert(jt != score_traceback.end()); - // std::vector next_alns; - // for (auto&& aln : alns) { - // anchor_extender(it->first, (it - 1)->first, std::move(aln), it->second, *jt, - // [&](Alignment&& next_aln) { - // next_alns.emplace_back(std::move(next_aln)); - // } - // ); - // } - // std::swap(next_alns, alns); - // } - - assert(jt == score_traceback.end()); for (auto&& aln : alns) { if (terminate()) @@ -203,11 +193,15 @@ void chain_anchors(const Query &query, for (size_t i = 0; i < chain_scores.size(); ++i) { const auto &[score, last, dist] = chain_scores[i]; - if (score > 0) + if (score > 0) { + assert(last <= end); + assert(last >= begin); best_chains.emplace_back(-score, (begin + i)->get_orientation(), i); + } } std::sort(best_chains.begin(), best_chains.end()); + common::logger->info("Best chain score: {}", -std::get<0>(best_chains[0])); sdsl::bit_vector used(chain_scores.size()); for (auto [nscore, orientation, i] : best_chains) { @@ -227,6 +221,7 @@ void chain_anchors(const Query &query, while (last != end) { last_anchor = last; std::tie(score, last, dist) = chain_scores[last - begin]; + assert(last_anchor != end); chain.emplace_back(last_anchor, dist); scores.emplace_back(score); if (used[last_anchor - begin]) { @@ -244,6 +239,8 @@ void chain_anchors(const Query &query, used[a_ptr - begin] = true; } + common::logger->info("Found chain with score {}", scores.back()); + extend_chain(chain, scores, anchor_extender, callback, terminate); } } diff --git a/metagraph/src/graph/alignment_redone/aln_cigar.cpp b/metagraph/src/graph/alignment_redone/aln_cigar.cpp index a330ab41d8..042d706aaf 100644 --- a/metagraph/src/graph/alignment_redone/aln_cigar.cpp +++ b/metagraph/src/graph/alignment_redone/aln_cigar.cpp @@ -1,8 +1,9 @@ #include "aln_cigar.hpp" #include "kmer/alphabets.hpp" +#include "graph/representation/succinct/boss.hpp" -namespace mtg::graph::align { +namespace mtg::graph::align_redone { OperatorTable initialize_opt_table() { @@ -126,19 +127,16 @@ bool Cigar::is_valid(std::string_view reference, std::string_view query) const { switch (op.first) { case CLIPPED: { - if ((ref_it != reference.begin() || alt_it != query.begin()) - && (ref_it != reference.end() || alt_it != query.end())) { - if (alt_it > query.end() - op.second) { - std::cerr << "Query too short after " - << Cigar::opt_to_char(op.first) << std::endl - << to_string() << std::endl - << reference << std::endl - << query << std::endl; - return false; - } - - alt_it += op.second; + if (alt_it > query.end() - op.second) { + std::cerr << "Query too short after " + << Cigar::opt_to_char(op.first) << std::endl + << to_string() << std::endl + << reference << std::endl + << query << std::endl; + return false; } + + alt_it += op.second; } break; case MATCH: case MISMATCH: { @@ -193,7 +191,7 @@ bool Cigar::is_valid(std::string_view reference, std::string_view query) const { alt_it += op.second; } break; case DELETION: { - if (i && cigar_[i - 1].first == INSERTION) { + if (i && cigar_[i - 1].first == INSERTION && *ref_it != boss::BOSS::kSentinel) { std::cerr << "DELETION after INSERTION" << std::endl << to_string() << std::endl << reference << std::endl diff --git a/metagraph/src/graph/alignment_redone/aln_cigar.hpp b/metagraph/src/graph/alignment_redone/aln_cigar.hpp index 2c3a3941ae..6adfeb2ddb 100644 --- a/metagraph/src/graph/alignment_redone/aln_cigar.hpp +++ b/metagraph/src/graph/alignment_redone/aln_cigar.hpp @@ -8,7 +8,7 @@ #include #include -namespace mtg::graph::align { +namespace mtg::graph::align_redone { class Cigar { public: diff --git a/metagraph/src/graph/alignment_redone/aln_match.cpp b/metagraph/src/graph/alignment_redone/aln_match.cpp index d5ef048dfc..8fbdefcec9 100644 --- a/metagraph/src/graph/alignment_redone/aln_match.cpp +++ b/metagraph/src/graph/alignment_redone/aln_match.cpp @@ -2,7 +2,7 @@ #include "common/logger.hpp" -namespace mtg::graph::align { +namespace mtg::graph::align_redone { using mtg::common::logger; @@ -87,7 +87,7 @@ std::string spell_path(const DeBruijnGraph &graph, void Anchor::append(const Anchor &other, const DeBruijnGraph *graph) { assert(query_.data() == other.query_.data()); assert(query_.data() + query_.size() == other.query_.data() + other.query_.size()); - assert(label_id_ == other.label_id_); + assert(label_class_ == other.label_class_); if (empty()) { *this = other; @@ -117,7 +117,7 @@ void Anchor::append(const Anchor &other, const DeBruijnGraph *graph) { path_.insert(path_.end(), other.path_.begin(), other.path_.end()); suffix_ = other.suffix_; - assert(is_valid()); + assert(!graph || is_spelling_valid(*graph)); } } // namespace mtg::graph::align \ No newline at end of file diff --git a/metagraph/src/graph/alignment_redone/aln_match.hpp b/metagraph/src/graph/alignment_redone/aln_match.hpp index 4120ab51b1..e289ef5894 100644 --- a/metagraph/src/graph/alignment_redone/aln_match.hpp +++ b/metagraph/src/graph/alignment_redone/aln_match.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include #include "aln_cigar.hpp" #include "aligner_config.hpp" @@ -8,7 +9,7 @@ #include "graph/annotated_dbg.hpp" #include "annotation/int_matrix/base/int_matrix.hpp" -namespace mtg::graph::align { +namespace mtg::graph::align_redone { std::string spell_path(const DeBruijnGraph &graph, const std::vector &path); @@ -25,6 +26,8 @@ class Match { std::string_view get_query() const { return query_; } std::string_view get_seed() const { return seed_; } + size_t size() const { return path_.size(); } + const std::vector& get_path() const { return path_; } virtual std::string_view get_spelling() const = 0; @@ -55,9 +58,7 @@ class Match { std::vector&& path, score_t score) : query_(query), orientation_(orientation), seed_(query_.substr(begin, end - begin)), - path_(std::move(path)), score_(score) { - assert(is_spelling_valid()); - } + path_(std::move(path)), score_(score) {} std::string_view query_; bool orientation_; @@ -66,6 +67,16 @@ class Match { score_t score_; }; +inline std::ostream& operator<<(std::ostream &out, const Match &a) { + out << a.get_orientation() << "," + << a.get_clipping() << "," + << a.get_seed() << "," + << a.get_spelling() << "," + << a.get_end_clipping() << "," + << a.get_score(); + return out; +} + class Anchor : public Match { public: using label_class_t = uint64_t; @@ -87,13 +98,23 @@ class Anchor : public Match { label_class_t label_class = nlabel, int64_t coord = ncoord) : Match(query, begin, end, orientation, - std::move(path), config.match_score(query)), - label_class_(label_class), coord_(coord), suffix_(std::move(suffix)) {} + std::move(path), 0), + label_class_(label_class), coord_(coord), suffix_(std::move(suffix)) { + score_ = config.score_cigar(seed_, query_, generate_cigar()); + } std::string_view get_spelling() const override final { return seed_; } std::string get_path_spelling() const override final { return std::string(seed_) + suffix_; } std::string_view get_trim_spelling() const override final { return suffix_; } + Cigar generate_cigar() const { + Cigar cigar; + cigar.append(Cigar::CLIPPED, get_clipping()); + cigar.append(Cigar::MATCH, seed_.size()); + cigar.append(Cigar::CLIPPED, get_end_clipping()); + return cigar; + } + size_t get_end_trim() const override final { return suffix_.size(); } void append(const Anchor &other, const DeBruijnGraph *graph = nullptr); @@ -113,19 +134,38 @@ class Alignment : public Match { Alignment& operator=(const Alignment&) = default; Alignment& operator=(Alignment&&) = default; + Alignment(const DeBruijnGraph &graph, + std::string_view query, + bool orientation, + std::vector&& path, + const DBGAlignerConfig &config, + Cigar&& cigar, + size_t end_trim = 0, + Anchor::label_class_t label_class = Anchor::nlabel) + : Match(query, cigar.get_clipping(), query.size() - cigar.get_end_clipping(), + orientation, std::move(path), 0), + end_trim_(end_trim), + path_spelling_(spell_path(graph, path_)), + cigar_(std::move(cigar)), + label_classes_(path_.size(), label_class) { + assert(path_spelling_.size() - graph.get_k() + 1 == path_.size()); + assert(get_spelling().size() + get_end_trim() == path_spelling_.size()); + assert(cigar_.get_end_clipping() == get_end_clipping()); + score_ = config.score_cigar(get_spelling(), query_, cigar_); + } + Alignment(const Anchor &anchor) - : Match(anchor.get_query(), anchor.get_seed().data() - anchor.get_query().data(), - anchor.get_query().data() + anchor.get_query().size() - (anchor.get_seed().data() + anchor.get_seed().size()), + : Match(anchor.get_query(), + anchor.get_clipping(), + anchor.get_clipping() + anchor.get_seed().size(), anchor.get_orientation(), std::vector(anchor.get_path()), anchor.get_score()), end_trim_(anchor.get_end_trim()), path_spelling_(anchor.get_path_spelling()), - cigar_(std::string(get_clipping(), 'S') - + std::string(seed_.size(), '=') - + std::string(get_end_clipping(), 'S')), + cigar_(anchor.generate_cigar()), label_classes_(path_.size(), anchor.get_label_class()) { - assert(path_spelling_.size() == query_.size() + end_trim_); + assert(get_spelling().size() + get_end_trim() == path_spelling_.size()); } const Cigar& get_cigar() const { return cigar_; } @@ -146,7 +186,12 @@ class Alignment : public Match { size_t end_trim_; std::string path_spelling_; Cigar cigar_; - std::vector label_classes_; + std::vector label_classes_; }; +inline std::ostream& operator<<(std::ostream &out, const Alignment &a) { + out << static_cast(a) << "," << a.get_cigar().to_string(); + return out; +} + } // namespace mtg::graph::align \ No newline at end of file diff --git a/metagraph/src/graph/alignment_redone/aln_query.hpp b/metagraph/src/graph/alignment_redone/aln_query.hpp index d029099007..0964eb2154 100644 --- a/metagraph/src/graph/alignment_redone/aln_query.hpp +++ b/metagraph/src/graph/alignment_redone/aln_query.hpp @@ -5,7 +5,7 @@ #include "graph/representation/base/sequence_graph.hpp" #include "common/seq_tools/reverse_complement.hpp" -namespace mtg::graph::align { +namespace mtg::graph::align_redone { class Query { public: diff --git a/metagraph/src/graph/alignment_redone/aln_seeder.cpp b/metagraph/src/graph/alignment_redone/aln_seeder.cpp index 54d454d134..e32d4a54c9 100644 --- a/metagraph/src/graph/alignment_redone/aln_seeder.cpp +++ b/metagraph/src/graph/alignment_redone/aln_seeder.cpp @@ -4,7 +4,7 @@ #include "common/logger.hpp" #include "common/vector_map.hpp" -namespace mtg::graph::align { +namespace mtg::graph::align_redone { std::vector Seeder::get_inexact_anchors() const { auto anchors = get_anchors(); @@ -72,6 +72,21 @@ void call_distances(const DeBruijnGraph &graph, } } +template +struct template_parameter; + +template