diff --git a/include/seqan3/alignment/matrix/detail/trace_directions.hpp b/include/seqan3/alignment/matrix/detail/trace_directions.hpp index 21f0aaa6196..4cf13d931cc 100644 --- a/include/seqan3/alignment/matrix/detail/trace_directions.hpp +++ b/include/seqan3/alignment/matrix/detail/trace_directions.hpp @@ -32,13 +32,17 @@ enum struct trace_directions : uint8_t //!\brief Trace comes from the diagonal entry. diagonal = 0b00001, //!\brief Trace comes from the above entry, while opening the gap. - up_open = 0b00010, + up_open = 0b00110, //!\brief Trace comes from the above entry. up = 0b00100, //!\brief Trace comes from the left entry, while opening the gap. - left_open = 0b01000, + left_open = 0b11000, //!\brief Trace comes from the left entry. - left = 0b10000 + left = 0b10000, + //!\brief Carry bit for the last up open even if it is not the maximum value. + carry_up_open = 0b00010, + //!\brief Carry bit for the last left open even if it is not the maximum value. + carry_left_open = 0b01000 }; } // namespace seqan3::detail diff --git a/include/seqan3/alignment/matrix/detail/trace_iterator_base.hpp b/include/seqan3/alignment/matrix/detail/trace_iterator_base.hpp index 6d9307a99aa..6342e812208 100644 --- a/include/seqan3/alignment/matrix/detail/trace_iterator_base.hpp +++ b/include/seqan3/alignment/matrix/detail/trace_iterator_base.hpp @@ -159,14 +159,14 @@ class trace_iterator_base { derived().go_up(matrix_iter); // Set new trace direction if last position was up_open. - if (static_cast(old_dir & trace_directions::up_open)) + if (static_cast(old_dir & trace_directions::carry_up_open)) set_trace_direction(*matrix_iter); } else if (current_direction == trace_directions::left) { derived().go_left(matrix_iter); // Set new trace direction if last position was left_open. - if (static_cast(old_dir & trace_directions::left_open)) + if (static_cast(old_dir & trace_directions::carry_left_open)) set_trace_direction(*matrix_iter); } else @@ -259,12 +259,11 @@ class trace_iterator_base { current_direction = trace_directions::diagonal; } - else if (static_cast(dir & trace_directions::up) || static_cast(dir & trace_directions::up_open)) + else if (static_cast(dir & trace_directions::up)) { current_direction = trace_directions::up; } - else if (static_cast(dir & trace_directions::left) - || static_cast(dir & trace_directions::left_open)) + else if (static_cast(dir & trace_directions::left)) { current_direction = trace_directions::left; } diff --git a/include/seqan3/alignment/pairwise/alignment_algorithm.hpp b/include/seqan3/alignment/pairwise/alignment_algorithm.hpp index 6ca4f3d2dd0..45f1dd81244 100644 --- a/include/seqan3/alignment/pairwise/alignment_algorithm.hpp +++ b/include/seqan3/alignment/pairwise/alignment_algorithm.hpp @@ -743,14 +743,23 @@ class alignment_algorithm : if constexpr (traits_t::compute_sequence_alignment) { auto trace_matrix_it = trace_debug_matrix.begin() + offset; - std::ranges::copy(column - | std::views::transform( - [](auto const & tpl) - { - using std::get; - return get<1>(tpl).current; - }), - trace_debug_matrix.begin() + offset); + std::ranges::copy(column | std::views::transform([](auto const & tpl) + { + using std::get; + auto trace = get<1>(tpl).current; + + if (auto _up = (trace & trace_directions::up_open); _up == trace_directions::carry_up_open) + trace = trace ^ trace_directions::carry_up_open; // remove silent up open signal + else if (_up == trace_directions::up_open) + trace = trace ^ trace_directions::up; // display up open only with single bit. + + if (auto _left = (trace & trace_directions::left_open); _left == trace_directions::carry_left_open) + trace = trace ^ trace_directions::carry_left_open; // remove silent left open signal + else if (_left == trace_directions::left_open) + trace = trace ^ trace_directions::left; // display left open only with single bit. + + return trace; + }), trace_debug_matrix.begin() + offset); } } diff --git a/include/seqan3/alignment/pairwise/detail/policy_affine_gap_with_trace_recursion.hpp b/include/seqan3/alignment/pairwise/detail/policy_affine_gap_with_trace_recursion.hpp index 3bb4ffe344a..f88d02962d8 100644 --- a/include/seqan3/alignment/pairwise/detail/policy_affine_gap_with_trace_recursion.hpp +++ b/include/seqan3/alignment/pairwise/detail/policy_affine_gap_with_trace_recursion.hpp @@ -79,7 +79,9 @@ class policy_affine_gap_with_trace_recursion : protected policy_affine_gap_recur ? (best_trace = previous_cell.vertical_trace(), vertical_score) : (best_trace |= previous_cell.vertical_trace(), diagonal_score); diagonal_score = (diagonal_score < horizontal_score) - ? (best_trace = previous_cell.horizontal_trace(), horizontal_score) + ? (best_trace = previous_cell.horizontal_trace() | + (best_trace & trace_directions::carry_up_open), + horizontal_score) : (best_trace |= previous_cell.horizontal_trace(), diagonal_score); score_type tmp = diagonal_score + gap_open_score; diff --git a/include/seqan3/alignment/pairwise/policy/affine_gap_policy.hpp b/include/seqan3/alignment/pairwise/policy/affine_gap_policy.hpp index 605af482121..6a970739248 100644 --- a/include/seqan3/alignment/pairwise/policy/affine_gap_policy.hpp +++ b/include/seqan3/alignment/pairwise/policy/affine_gap_policy.hpp @@ -97,7 +97,9 @@ class affine_gap_policy { tmp = (tmp < score_cell.up) ? (trace_cell.current = trace_cell.up, score_cell.up) : (trace_cell.current = trace_directions::diagonal | trace_cell.up, tmp); - tmp = (tmp < score_cell.r_left) ? (trace_cell.current = trace_cell.r_left, score_cell.r_left) + tmp = (tmp < score_cell.r_left) ? (trace_cell.current = trace_cell.r_left | (trace_cell.current & + trace_directions::carry_up_open), + score_cell.r_left) : (trace_cell.current |= trace_cell.r_left, tmp); } else diff --git a/test/unit/alignment/matrix/detail/debug_stream_trace_directions_test.cpp b/test/unit/alignment/matrix/detail/debug_stream_trace_directions_test.cpp index a238c792da8..98c8cf4a8e6 100644 --- a/test/unit/alignment/matrix/detail/debug_stream_trace_directions_test.cpp +++ b/test/unit/alignment/matrix/detail/debug_stream_trace_directions_test.cpp @@ -18,8 +18,8 @@ static constexpr seqan3::detail::trace_directions N = seqan3::detail::trace_dire static constexpr seqan3::detail::trace_directions D = seqan3::detail::trace_directions::diagonal; static constexpr seqan3::detail::trace_directions u = seqan3::detail::trace_directions::up; static constexpr seqan3::detail::trace_directions l = seqan3::detail::trace_directions::left; -static constexpr seqan3::detail::trace_directions U = seqan3::detail::trace_directions::up_open; -static constexpr seqan3::detail::trace_directions L = seqan3::detail::trace_directions::left_open; +static constexpr seqan3::detail::trace_directions U = seqan3::detail::trace_directions::carry_up_open; +static constexpr seqan3::detail::trace_directions L = seqan3::detail::trace_directions::carry_left_open; TEST(debug_stream_test, ascii) { diff --git a/test/unit/alignment/pairwise/fixture/alignment_fixture.hpp b/test/unit/alignment/pairwise/fixture/alignment_fixture.hpp index 83057218416..e5d2bd79911 100644 --- a/test/unit/alignment/pairwise/fixture/alignment_fixture.hpp +++ b/test/unit/alignment/pairwise/fixture/alignment_fixture.hpp @@ -19,8 +19,8 @@ static constexpr auto N = seqan3::detail::trace_directions::none; static constexpr auto D = seqan3::detail::trace_directions::diagonal; static constexpr auto u = seqan3::detail::trace_directions::up; static constexpr auto l = seqan3::detail::trace_directions::left; -static constexpr auto U = seqan3::detail::trace_directions::up_open; -static constexpr auto L = seqan3::detail::trace_directions::left_open; +static constexpr auto U = seqan3::detail::trace_directions::carry_up_open; +static constexpr auto L = seqan3::detail::trace_directions::carry_left_open; static constexpr auto DU = D | U; static constexpr auto UL = U | L; static constexpr auto DL = D | L; diff --git a/test/unit/alignment/pairwise/fixture/global_affine_unbanded.hpp b/test/unit/alignment/pairwise/fixture/global_affine_unbanded.hpp index 9f972410a35..b86b1321e41 100644 --- a/test/unit/alignment/pairwise/fixture/global_affine_unbanded.hpp +++ b/test/unit/alignment/pairwise/fixture/global_affine_unbanded.hpp @@ -417,6 +417,76 @@ static auto dna4_match_4_mismatch_5_gap_1_open_10_both_empty = []() }; }(); +// [ISSUE #3043] Wrong alignment score https://github.com/seqan/seqan3/issues/3043 +static auto issue_3043 = []() +{ + using seqan3::operator""_dna4; + + auto make_scheme = [] () + { + seqan3::nucleotide_scoring_scheme scheme{}; + scheme.score('A'_dna4, 'A'_dna4) = 145; + scheme.score('A'_dna4, 'T'_dna4) = -144; + scheme.score('A'_dna4, 'G'_dna4) = -153; + scheme.score('A'_dna4, 'C'_dna4) = -152; + + scheme.score('T'_dna4, 'A'_dna4) = -144; + scheme.score('T'_dna4, 'T'_dna4) = 144; + scheme.score('T'_dna4, 'G'_dna4) = -143; + scheme.score('T'_dna4, 'C'_dna4) = -140; + + scheme.score('G'_dna4, 'A'_dna4) = -153; + scheme.score('G'_dna4, 'T'_dna4) = -143; + scheme.score('G'_dna4, 'G'_dna4) = 144; + scheme.score('G'_dna4, 'C'_dna4) = -145; + + scheme.score('C'_dna4, 'A'_dna4) = -152; + scheme.score('C'_dna4, 'T'_dna4) = -140; + scheme.score('C'_dna4, 'G'_dna4) = -145; + scheme.score('C'_dna4, 'C'_dna4) = 144; + + return scheme; + }; + + return alignment_fixture + { + "GGCAAGAA"_dna4, + "CGAAGC"_dna4, + seqan3::align_cfg::method_global{} | seqan3::align_cfg::gap_cost_affine{seqan3::align_cfg::open_score{-52}, + seqan3::align_cfg::extension_score{-58}} + | seqan3::align_cfg::scoring_scheme{make_scheme()}, + 74, + "GGCAAGAA--", + "--C--GAAGC", + /*.sequence1_begin_position = */ 0u, + /*.sequence2_begin_position = */ 0u, + /*.sequence1_end_position = */ 8u, + /*.sequence2_end_position = */ 6u, + std::vector + { + // e + /*e*/ 0,-110,-168,-226,-284,-342,-400,-458,-516, + -110,-145,-255, -24,-134,-192,-250,-308,-366, + -168, 34, -1,-111,-169,-227, -48,-158,-216, + -226, -76,-111,-153, 34, -24,-134, 97, -13, + -284,-134,-169,-250, -8, 179, 69, 11, 242, + -342,-140, 10,-100,-118, 69, 323, 213, 155, + -400,-250,-100, 154, 44, 11, 213, 171, 74 + }, + std::vector + { + // e,G ,G ,C ,A ,A ,G ,A ,A , + /*e*/ N,L ,l ,l ,l ,l ,l ,l ,l , + /*C*/ U,DUL,DUL,DUl,L ,l ,l ,l ,l , + /*G*/ u,DUL,DuL,L ,l ,l ,DUl,L ,l , + /*A*/ u,UL ,UL ,DuL,DUL,DUL,l ,DUl,DUL, + /*A*/ u,uL ,uL ,uL ,DUl,DUL,L ,DUl,DUl, + /*G*/ u,DuL,DuL,L ,Ul ,Ul ,DUL,L ,l , + /*C*/ u,uL ,UL ,DUL,L ,ul ,Ul ,DUL,uL + } + }; +}(); + // ---------------------------------------------------------------------------- // alignment fixtures using amino acid alphabet // ---------------------------------------------------------------------------- diff --git a/test/unit/alignment/pairwise/global_affine_unbanded_test.cpp b/test/unit/alignment/pairwise/global_affine_unbanded_test.cpp index 732015b3940..0bd930b5e4a 100644 --- a/test/unit/alignment/pairwise/global_affine_unbanded_test.cpp +++ b/test/unit/alignment/pairwise/global_affine_unbanded_test.cpp @@ -28,7 +28,10 @@ using pairwise_global_affine_unbanded_testing_types = ::testing::Types< pairwise_alignment_fixture< &seqan3::test::alignment::fixture::global::affine::unbanded::dna4_match_4_mismatch_5_gap_1_open_10_seq2_empty>, pairwise_alignment_fixture< - &seqan3::test::alignment::fixture::global::affine::unbanded::dna4_match_4_mismatch_5_gap_1_open_10_both_empty>>; + &seqan3::test::alignment::fixture::global::affine::unbanded::dna4_match_4_mismatch_5_gap_1_open_10_both_empty>, + pairwise_alignment_fixture< + &seqan3::test::alignment::fixture::global::affine::unbanded::issue_3043> + >; INSTANTIATE_TYPED_TEST_SUITE_P(pairwise_global_affine_unbanded, pairwise_alignment_test,