diff --git a/include/seqan3/alignment/pairwise/alignment_algorithm.hpp b/include/seqan3/alignment/pairwise/alignment_algorithm.hpp index c5c5718ddf..898be56550 100644 --- a/include/seqan3/alignment/pairwise/alignment_algorithm.hpp +++ b/include/seqan3/alignment/pairwise/alignment_algorithm.hpp @@ -619,7 +619,11 @@ class alignment_algorithm : row_index_type{this->alignment_state.optimum.row_index}}; // At some point this needs to be refactored so that it is not necessary to adapt the coordinate. if constexpr (traits_t::is_banded) + { res.end_positions.second += res.end_positions.first - this->trace_matrix.band_col_index; + res.end_positions.first = this->to_original_sequence1_position(res.end_positions.first); + res.end_positions.second = this->to_original_sequence2_position(res.end_positions.second); + } } if constexpr (traits_t::compute_begin_positions) @@ -631,8 +635,10 @@ class alignment_algorithm : detail::row_index_type{this->alignment_state.optimum.row_index}, detail::column_index_type{this->alignment_state.optimum.column_index}}; auto trace_res = builder(this->trace_matrix.trace_path(optimum_coordinate)); - res.begin_positions.first = trace_res.first_sequence_slice_positions.first; - res.begin_positions.second = trace_res.second_sequence_slice_positions.first; + res.begin_positions.first = + this->to_original_sequence1_position(trace_res.first_sequence_slice_positions.first); + res.begin_positions.second = + this->to_original_sequence2_position(trace_res.second_sequence_slice_positions.first); if constexpr (traits_t::compute_sequence_alignment) res.alignment = std::move(trace_res.alignment); @@ -697,8 +703,10 @@ class alignment_algorithm : if constexpr (traits_t::compute_end_positions) { - res.end_positions.first = this->alignment_state.optimum.column_index[simd_index]; - res.end_positions.second = this->alignment_state.optimum.row_index[simd_index]; + res.end_positions.first = + this->to_original_sequence1_position(this->alignment_state.optimum.column_index[simd_index]); + res.end_positions.second = + this->to_original_sequence2_position(this->alignment_state.optimum.row_index[simd_index]); } callback(std::move(res)); diff --git a/include/seqan3/alignment/pairwise/policy/alignment_matrix_policy.hpp b/include/seqan3/alignment/pairwise/policy/alignment_matrix_policy.hpp index 29c85083b5..d7d50924de 100644 --- a/include/seqan3/alignment/pairwise/policy/alignment_matrix_policy.hpp +++ b/include/seqan3/alignment/pairwise/policy/alignment_matrix_policy.hpp @@ -140,25 +140,24 @@ class alignment_matrix_policy * band starts in the origin and ends in the sink. */ template - constexpr auto slice_sequences(sequence1_t & sequence1, - sequence2_t & sequence2, - align_cfg::band_fixed_size const & band) const noexcept + constexpr auto + slice_sequences(sequence1_t & sequence1, sequence2_t & sequence2, align_cfg::band_fixed_size const & band) noexcept { size_t seq1_size = std::ranges::distance(sequence1); size_t seq2_size = std::ranges::distance(sequence2); auto trim_sequence1 = [&]() constexpr { - size_t begin_pos = std::max(band.lower_diagonal - 1, 0); + seq1_slice_offset = std::max(band.lower_diagonal - 1, 0); size_t end_pos = std::min(band.upper_diagonal + seq2_size, seq1_size); - return sequence1 | views::slice(begin_pos, end_pos); + return sequence1 | views::slice(seq1_slice_offset, end_pos); }; auto trim_sequence2 = [&]() constexpr { - size_t begin_pos = std::abs(std::min(band.upper_diagonal + 1, 0)); + seq2_slice_offset = std::abs(std::min(band.upper_diagonal + 1, 0)); size_t end_pos = std::min(seq1_size - band.lower_diagonal, seq2_size); - return sequence2 | views::slice(begin_pos, end_pos); + return sequence2 | views::slice(seq2_slice_offset, end_pos); }; return std::tuple{trim_sequence1(), trim_sequence2()}; @@ -192,10 +191,37 @@ class alignment_matrix_policy ++trace_matrix_iter; } + /*!\brief Converts a sliced position in the alignment matrix to the corresponding position in the original sequence 1. + * + * \param position The position in the sliced alignment matrix. + * \return The corresponding position in the original sequence 1. + * + * Only changes the position if the sequence was sliced due to band configuration. + */ + constexpr size_t to_original_sequence1_position(size_t position) const noexcept + { + return position + seq1_slice_offset; + } + + /*!\brief Converts a sliced position in the alignment matrix to the corresponding position in the original sequence 2. + * + * \param position The position in the sliced alignment matrix. + * \return The corresponding position in the original sequence 2. + * + * Only changes the position if the sequence was sliced due to band configuration. + */ + constexpr size_t to_original_sequence2_position(size_t position) const noexcept + { + return position + seq2_slice_offset; + } + score_matrix_t score_matrix{}; //!< The scoring matrix. trace_matrix_t trace_matrix{}; //!< The trace matrix if needed. typename score_matrix_t::iterator score_matrix_iter{}; //!< The matrix iterator over the score matrix. typename trace_matrix_t::iterator trace_matrix_iter{}; //!< The matrix iterator over the trace matrix. + + size_t seq1_slice_offset{}; //!< The offset of the first sequence slice. + size_t seq2_slice_offset{}; //!< The offset of the second sequence slice; }; } // namespace seqan3::detail diff --git a/test/unit/alignment/pairwise/global_affine_banded_test.cpp b/test/unit/alignment/pairwise/global_affine_banded_test.cpp index 0b396f1a26..fd5494bc2e 100644 --- a/test/unit/alignment/pairwise/global_affine_banded_test.cpp +++ b/test/unit/alignment/pairwise/global_affine_banded_test.cpp @@ -68,3 +68,35 @@ TEST_F(pairwise_global_affine_banded, invalid_band_last_cell_not_covered) fixture.config | seqan3::align_cfg::output_score{}); EXPECT_THROW(result_range.begin(), seqan3::invalid_alignment_configuration); } + +TEST(banded_alignment_issue3266_test, wrong_begin_and_end_position) +{ + using namespace seqan3; + using namespace seqan3::literals; + + auto const configGeneral = align_cfg::scoring_scheme{nucleotide_scoring_scheme{match_score{1}, mismatch_score{-1}}} + | align_cfg::gap_cost_affine{align_cfg::open_score{0}, align_cfg::extension_score{-1}} + | align_cfg::method_global{align_cfg::free_end_gaps_sequence1_leading{true}, + align_cfg::free_end_gaps_sequence2_leading{true}, + align_cfg::free_end_gaps_sequence1_trailing{true}, + align_cfg::free_end_gaps_sequence2_trailing{true}}; + + auto const configBanded = + configGeneral | align_cfg::band_fixed_size{align_cfg::lower_diagonal{-40}, align_cfg::upper_diagonal{-20}}; + + //0 1 2 3 4 + //01234567890123456789012345678901234567890 + std::pair p{"CGTCTA"_dna4, "AAACCCGGGTTTAAACCCGGGTTTCGTGTACCCCCCCCCCC"_dna4}; + // CGTCTA + + auto general_results = align_pairwise(p, configGeneral); + auto general_res = *std::ranges::begin(general_results); + auto banded_results = align_pairwise(p, configBanded); + auto banded_res = *std::ranges::begin(banded_results); + + EXPECT_EQ(general_res.score(), banded_res.score()); + EXPECT_EQ(general_res.sequence1_begin_position(), banded_res.sequence1_begin_position()); + EXPECT_EQ(general_res.sequence2_begin_position(), banded_res.sequence2_begin_position()); + EXPECT_EQ(general_res.sequence1_end_position(), banded_res.sequence1_end_position()); + EXPECT_EQ(general_res.sequence2_end_position(), banded_res.sequence2_end_position()); +}