SeqAn3
The Modern C++ library for sequence analysis.
alignment_trace_algorithms.hpp
Go to the documentation of this file.
1 // -----------------------------------------------------------------------------------------------------
2 // Copyright (c) 2006-2019, Knut Reinert & Freie Universität Berlin
3 // Copyright (c) 2016-2019, Knut Reinert & MPI für molekulare Genetik
4 // This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License
5 // shipped with this file and also available at: https://github.com/seqan/seqan3/blob/master/LICENSE.md
6 // -----------------------------------------------------------------------------------------------------
7 
13 #pragma once
14 
15 #include <deque>
16 #include <vector>
17 
25 
26 namespace seqan3::detail
27 {
28 
36  template <typename trace_matrix_t>
38  requires Matrix<remove_cvref_t<trace_matrix_t>> &&
41 inline alignment_coordinate alignment_front_coordinate(trace_matrix_t && matrix,
42  alignment_coordinate const back_coordinate)
43 {
44  constexpr auto D = trace_directions::diagonal;
45  constexpr auto L = trace_directions::left;
46  constexpr auto U = trace_directions::up;
47  size_t row = back_coordinate.second + 1;
48  size_t col = back_coordinate.first + 1;
49 
50  assert(row < matrix.rows());
51  assert(col < matrix.cols());
52 
53  while (true)
54  {
55  trace_directions dir = matrix.at(row, col);
56  if ((dir & L) == L)
57  {
58  col = std::max<size_t>(col, 1) - 1;
59  }
60  else if ((dir & U) == U)
61  {
62  row = std::max<size_t>(row, 1) - 1;
63  }
64  else if ((dir & D) == D)
65  {
66  row = std::max<size_t>(row, 1) - 1;
67  col = std::max<size_t>(col, 1) - 1;
68  }
69  else
70  {
71 #ifndef NDEBUG
72  if (!(row == 0 || col == 0))
73  throw std::logic_error{"Unknown seqan3::trace_direction in an inner cell of the trace matrix."};
74 #endif
75  break;
76  }
77  }
78 
79  return {column_index_type{col}, row_index_type{row}};
80 }
81 
95 template <
96  TupleLike alignment_t,
97  typename database_t,
98  typename query_t,
99  typename trace_matrix_t>
101  requires Matrix<remove_cvref_t<trace_matrix_t>> &&
102  std::Same<typename remove_cvref_t<trace_matrix_t>::entry_type, trace_directions> &&
103  detail::all_satisfy_aligned_seq<detail::tuple_type_list_t<alignment_t>>
105 inline alignment_t alignment_trace(database_t && database,
106  query_t && query,
107  trace_matrix_t && matrix,
108  alignment_coordinate const back_coordinate,
109  alignment_coordinate const front_coordinate)
110 {
111  constexpr auto N = trace_directions::none;
112  constexpr auto D = trace_directions::diagonal;
113  constexpr auto L = trace_directions::left;
114  constexpr auto U = trace_directions::up;
115  size_t col = back_coordinate.first + 1;
116  size_t row = back_coordinate.second + 1;
117 
118  assert(row <= query.size());
119  assert(col <= database.size());
120  assert(row < matrix.rows());
121  assert(col < matrix.cols());
122 
123  alignment_t aligned_seq{};
124  assign_unaligned(std::get<0>(aligned_seq), view::slice(database, front_coordinate.first, col));
125  assign_unaligned(std::get<1>(aligned_seq), view::slice(query, front_coordinate.second, row));
126  auto end_aligned_db = std::ranges::cend(std::get<0>(aligned_seq));
127  auto end_aligned_qy = std::ranges::cend(std::get<1>(aligned_seq));
128 
129  if (matrix.at(0, 0) != N)
130  throw std::logic_error{"End trace must be NONE"};
131 
132  while (true)
133  {
134  trace_directions dir = matrix.at(row, col);
135  if ((dir & L) == L)
136  {
137  col = std::max<size_t>(col, 1) - 1;
138  --end_aligned_db;
139  insert_gap(std::get<1>(aligned_seq), end_aligned_qy);
140  }
141  else if ((dir & U) == U)
142  {
143  row = std::max<size_t>(row, 1) - 1;
144  insert_gap(std::get<0>(aligned_seq), end_aligned_db);
145  --end_aligned_qy;
146  }
147  else if ((dir & D) == D)
148  {
149  row = std::max<size_t>(row, 1) - 1;
150  col = std::max<size_t>(col, 1) - 1;
151  --end_aligned_db;
152  --end_aligned_qy;
153  }
154  else
155  {
156 #ifndef NDEBUG
157  if (!(row == 0 || col == 0))
158  throw std::logic_error{"Unknown seqan3::trace_direction in an inner cell of the trace matrix."};
159 #endif
160  break;
161  }
162  }
163 
164  return aligned_seq;
165 }
166 
167 } // namespace seqan3::detail
void assign_unaligned(aligned_seq_t &aligned_seq, unaligned_sequence_type &&unaligned_seq)
An implementation of seqan3::AlignedSequence::assign_unaligned_sequence for sequence containers...
Definition: aligned_sequence_concept.hpp:345
aligned_seq_t::iterator insert_gap(aligned_seq_t &aligned_seq, typename aligned_seq_t::const_iterator pos_it)
An implementation of seqan3::AlignedSequence::insert_gap for sequence containers. ...
Definition: aligned_sequence_concept.hpp:229
Includes the AlignedSequence and the related insert_gap and erase_gap functions to enable stl contain...
constexpr auto slice
A view adaptor that returns a half-open interval on the underlying range.
Definition: slice.hpp:144
Provides seqan3::gapped.
Provides seqan3::view::all.
Provides the declaration of seqan3::detail::trace_directions.
The concept std::Same<T, U> is satisfied if and only if T and U denote the same type.
Definition: aligned_sequence_concept.hpp:35
Provides seqan3::detail::alignment_coordinate.
Provides various transformation traits used by the range module.
::ranges::cend cend
Alias for ranges::cend. Returns an iterator to the end of a range.
Definition: ranges:214
Provides seqan3::detail::Matrix.