Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[FIX] Wrong alignment score #3043 #3098

Merged
merged 1 commit into from
Nov 14, 2022
+118 −26
Merged

[FIX] Wrong alignment score #3043 #3098

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
[FIX] Alignment Wrong alignment score #3043 
Added carry bit to up and left open trace directions in order to capture original open signal even if it is overwritten by the orhtogonal direction.
In the particular issue the left extension overrites the up open signal in the last column and fifth row.
The traceback now follows the vertical or horizontal direction unitl it encounters the carry bit signal.
@rrahn @eseiler
rrahn authored and eseiler committed Nov 14, 2022
commit 4fe548913e96d3f908dd524bd3dc13b48f87bfa4
10 changes: 7 additions & 3 deletions include/seqan3/alignment/matrix/detail/trace_directions.hpp
View file
Original file line number Diff line number Diff line change
@@ -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
9 changes: 4 additions & 5 deletions include/seqan3/alignment/matrix/detail/trace_iterator_base.hpp
View file
Original file line number Diff line number Diff line change
@@ -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<bool>(old_dir & trace_directions::up_open))
if (static_cast<bool>(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<bool>(old_dir & trace_directions::left_open))
if (static_cast<bool>(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<bool>(dir & trace_directions::up) || static_cast<bool>(dir & trace_directions::up_open))
else if (static_cast<bool>(dir & trace_directions::up))
{
current_direction = trace_directions::up;
}
else if (static_cast<bool>(dir & trace_directions::left)
|| static_cast<bool>(dir & trace_directions::left_open))
else if (static_cast<bool>(dir & trace_directions::left))
{
current_direction = trace_directions::left;
}
30 changes: 22 additions & 8 deletions include/seqan3/alignment/pairwise/alignment_algorithm.hpp
View file
Original file line number Diff line number Diff line change
@@ -743,14 +743,28 @@ 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);
}
}

8 changes: 5 additions & 3 deletions include/seqan3/alignment/pairwise/detail/policy_affine_gap_with_trace_recursion.hpp
View file
Original file line number Diff line number Diff line change
@@ -78,9 +78,11 @@ class policy_affine_gap_with_trace_recursion : protected policy_affine_gap_recur
diagonal_score = (diagonal_score < vertical_score)
? (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(), diagonal_score);
diagonal_score =
(diagonal_score < horizontal_score)
? (best_trace = previous_cell.horizontal_trace() | (best_trace & trace_directions::carry_up_open),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(best_trace & trace_directions::carry_left_open) is not needed in the above (L 79) case, because the horizontal trace is always set

horizontal_score)
: (best_trace |= previous_cell.horizontal_trace(), diagonal_score);

score_type tmp = diagonal_score + gap_open_score;
vertical_score += gap_extension_score;
6 changes: 4 additions & 2 deletions include/seqan3/alignment/pairwise/policy/affine_gap_policy.hpp
View file
Original file line number Diff line number Diff line change
@@ -97,8 +97,10 @@ 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)
: (trace_cell.current |= trace_cell.r_left, tmp);
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
{
4 changes: 2 additions & 2 deletions test/unit/alignment/matrix/detail/debug_stream_trace_directions_test.cpp
View file
Original file line number Diff line number Diff line change
@@ -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;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Changes as to not need to touch every test

static constexpr seqan3::detail::trace_directions L = seqan3::detail::trace_directions::carry_left_open;

TEST(debug_stream_test, ascii)
{
4 changes: 2 additions & 2 deletions test/unit/alignment/pairwise/fixture/alignment_fixture.hpp
View file
Original file line number Diff line number Diff line change
@@ -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;
70 changes: 70 additions & 0 deletions test/unit/alignment/pairwise/fixture/global_affine_unbanded.hpp
View file
Original file line number Diff line number Diff line change
@@ -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<int32_t> 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
// ----------------------------------------------------------------------------
3 changes: 2 additions & 1 deletion test/unit/alignment/pairwise/global_affine_unbanded_test.cpp
View file
Original file line number Diff line number Diff line change
@@ -28,7 +28,8 @@ 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,