// ==========================================================================
// SeqAn - The Library for Sequence Analysis
// ==========================================================================
// Copyright (c) 2006-2018, Knut Reinert, FU Berlin
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
// * Neither the name of Knut Reinert or the FU Berlin nor the names of
// its contributors may be used to endorse or promote products derived
// from this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL KNUT REINERT OR THE FU BERLIN BE LIABLE
// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
// OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
// DAMAGE.
//
// ==========================================================================
// Author: Andreas Gogol-Doering <andreas.doering@mdc-berlin.de>
// Author: Anne-Katrin Emde <anne-katrin.emde@fu-berlin.de>
// ==========================================================================
// Implementation of the Waterman-Eggert algorithm, sometimes also called
// Smith-Waterman algorithm with declumping.
// ==========================================================================
#ifndef SEQAN_INCLUDE_SEQAN_ALIGN_LOCAL_ALIGNMENT_WATERMAN_EGGERT_IMPL_H_
#define SEQAN_INCLUDE_SEQAN_ALIGN_LOCAL_ALIGNMENT_WATERMAN_EGGERT_IMPL_H_
namespace seqan {
// ============================================================================
// Forwards
// ============================================================================
// ============================================================================
// Tags, Classes, Enums
// ============================================================================
// ----------------------------------------------------------------------------
// Helper Class ScoreAndID
// ----------------------------------------------------------------------------
// TODO(holtgrew): Why is Pair<> not enough?
// TODO(holtgrew): Rename to ScoreAndID_ (with trailing underscore)
// Simple class that stores a value with an ID.
template <typename TValue, typename TID>
class ScoreAndID
{
public:
TValue value_;
TID id_;
ScoreAndID() : value_(std::numeric_limits<TValue>::min()), id_(std::numeric_limits<TValue>::max())
{}
ScoreAndID(TValue score, TID id_pos)
{
value_ = score;
id_ = id_pos;
}
inline bool operator>(ScoreAndID const & other) const
{
return value_ > other.value_;
}
inline bool operator<(ScoreAndID const & other) const
{
return value_ < other.value_;
}
};
// ----------------------------------------------------------------------------
// Class LocalAlignmentFinder
// ----------------------------------------------------------------------------
template <typename TScoreValue = int>
class LocalAlignmentFinder
{
public:
typedef Matrix<TScoreValue> TMatrix;
typedef typename Position<TMatrix>::Type TMatrixPosition;
typedef typename Size<TMatrix>::Type TSize;
typedef ScoreAndID<TScoreValue,TMatrixPosition> TPQEntry;
typedef Iter<TMatrix,PositionIterator> TMatrixIterator;
typedef PriorityType<TPQEntry> TPriorityQ;
typedef String<bool> TBoolMatrix;
//DP-matrix
TMatrix matrix;
//matrix that memorizes the cells from which not to go diagonal
TBoolMatrix forbidden;
//priority queue for quickly finding the maximum score in the DP-matrix
TPriorityQ pQ;
//position of maximum score (where traceback is started from)
TMatrixPosition bestEndPos;
//position where traceback ended and where declumping begins
TMatrixPosition bestBeginPos;
//traceback path that is set to forbidden while declumping
AlignTraceback<TSize> trace;
bool needReinit; //true: call "smithWaterman", false: call "smithWatermanGetNext"
LocalAlignmentFinder() : bestEndPos(0), bestBeginPos(0), needReinit(true)
{}
// TODO(holtgrew): Remove and replace all occurrences with default constructor.
template<typename TAlign>
LocalAlignmentFinder(TAlign const &) : bestEndPos(0), bestBeginPos(0), needReinit(true)
{}
};
// ============================================================================
// Metafunctions
// ============================================================================
// ============================================================================
// Functions
// ============================================================================
// ----------------------------------------------------------------------------
// Function _initLocalAlignmentFinder()
// ----------------------------------------------------------------------------
template<typename TSequenceH, typename TSequenceV, typename TScoreValue, typename TTag>
void
_initLocalAlignmentFinder(TSequenceH const & seqH,
TSequenceV const & seqV,
LocalAlignmentFinder<TScoreValue> & finder,
TTag) {
typedef LocalAlignmentFinder<TScoreValue> TFinder;
typedef typename TFinder::TMatrix TMatrix;
typedef typename Size<TMatrix>::Type TSize;
TSize len0 = length(seqH);
TSize len1 = length(seqV);
setDimension(finder.matrix, 2);
setLength(finder.matrix, 0, len0 + 1);
setLength(finder.matrix, 1, len1 + 1);
resize(finder.matrix);
resize(finder.forbidden, (len0 + 1) * (len1 + 1), false);
finder.bestEndPos = std::numeric_limits<typename TFinder::TMatrixPosition>::max();
finder.bestBeginPos = std::numeric_limits<typename TFinder::TMatrixPosition>::max();
}
// ----------------------------------------------------------------------------
// Function clear()
// ----------------------------------------------------------------------------
template <typename TScoreValue>
void clear(LocalAlignmentFinder<TScoreValue> & sw_finder)
{
sw_finder.needReinit = true;
}
// ----------------------------------------------------------------------------
// Function getScore()
// ----------------------------------------------------------------------------
template <typename TScoreValue>
TScoreValue getScore(LocalAlignmentFinder<TScoreValue> const & sw)
{
typedef LocalAlignmentFinder<TScoreValue> TFinder;
if(sw.bestEndPos != std::numeric_limits<typename TFinder::TMatrixPosition>::max())
return getValue(const_cast<typename TFinder::TMatrix &>(sw.matrix), sw.bestEndPos);
return 0;
}
// ----------------------------------------------------------------------------
// Function _smithWatermanGetMatrix()
// ----------------------------------------------------------------------------
template <typename TScoreValue, typename TScoreSpec, typename TStringH, typename TStringV>
TScoreValue
_smithWatermanGetMatrix(LocalAlignmentFinder<TScoreValue> & sw,
TStringH const & strH,
TStringV const & strV,
Score<TScoreValue, TScoreSpec> const & score_,
TScoreValue cutoff)
{
// typedefs
typedef Matrix<TScoreValue> TMatrix;
typedef typename Position<TMatrix>::Type TMatrixPosition;
typedef Iter<TMatrix,PositionIterator> TMatrixIterator;
typedef typename Iterator<TStringH const, Rooted>::Type TStringIteratorH;
//typedef typename Value<TStringH const>::Type TValueH;
typedef typename Iterator<TStringV const, Rooted>::Type TStringIteratorV;
typedef typename Value<TStringV const>::Type TValueV;
//-------------------------------------------------------------------------
//define some variables
// TSize str1_length = length(strH);
// TSize str2_length = length(strV);
TStringIteratorH x_begin = begin(strH) - 1;
TStringIteratorH x_end = end(strH) - 1;
TStringIteratorV y_begin = begin(strV) - 1;
TStringIteratorV y_end = end(strV) - 1;
TStringIteratorH x = x_end;
TStringIteratorV y;
TScoreValue score_gap = scoreGapExtend(score_);
TScoreValue h = 0;
TScoreValue v = 0;
TMatrixIterator col_ = end(sw.matrix) - 1;
TMatrixIterator finger1;
TMatrixIterator finger2;
//-------------------------------------------------------------------------
// init
finger1 = col_;
*finger1 = 0;
//std::cout <<" ";
for (x = x_end; x != x_begin; --x)
{
goPrevious(finger1, 0);
*finger1 = 0;
}
//-------------------------------------------------------------------------
//fill matrix
for (y = y_end; y != y_begin; --y)
{
TValueV cy = *y;
h = 0;
v = 0;
finger2 = col_; //points to last column
goPrevious(col_, 1); //points to this column
finger1 = col_;
*finger1 = v;
for (x = x_end; x != x_begin; --x)
{
goPrevious(finger1, 0);
goPrevious(finger2, 0);
TScoreValue currScore = score(score_, *x, cy);
if (*x == cy)
{
v = h + currScore;
h = *finger2;
}
else
{
TScoreValue s1 = h + currScore;
h = *finger2;
TScoreValue s2 = score_gap + ((h > v) ? h : v);
v = (s1 > s2) ? s1 : s2;
if (v < 0) v = 0;
}
*finger1 = v;
if (v >= cutoff)
{
push(sw.pQ,ScoreAndID<TScoreValue,TMatrixPosition>(v,position(finger1)));
}
}
}
// check if any scores >= cutoff were found
if(!empty(sw.pQ))
{
ScoreAndID<TScoreValue,TMatrixPosition> best = top(sw.pQ);
v = getValue(sw.matrix,best.id_);
sw.bestEndPos = best.id_;
}
else
v=0;
return v;
}
// ----------------------------------------------------------------------------
// Function _smithWatermanDeclump()
// ----------------------------------------------------------------------------
// declumping
template <typename TScoreValue,
typename TSequenceH, typename TGapsSpecH,
typename TSequenceV, typename TGapsSpecV,
typename TScoreSpec>
void
_smithWatermanDeclump(LocalAlignmentFinder<TScoreValue> & sw ,
Gaps<TSequenceH, TGapsSpecH> & gapsH,
Gaps<TSequenceV, TGapsSpecV> & gapsV,
Score<TScoreValue, TScoreSpec> const & score_)
{
//-------------------------------------------------------------------------
//typedefs
//typedef typename LocalAlignmentFinder<TScoreValue>::TMatrixPosition TMatrixPosition;
typedef typename LocalAlignmentFinder<TScoreValue>::TMatrix TMatrix;
typedef Iter<TMatrix, PositionIterator> TMatrixIterator;
typedef Gaps<TSequenceH, TGapsSpecH> TGapsH;
typedef typename Iterator<TGapsH, Rooted>::Type TGapsHIter;
typedef typename Iterator<TSequenceH, Rooted>::Type TSequenceHIter;
//typedef typename Value<TSequenceH>::Type TValueH;
typedef Gaps<TSequenceV, TGapsSpecV> TGapsV;
typedef typename Iterator<TGapsV, Rooted>::Type TGapsVIter;
typedef typename Iterator<TSequenceV, Rooted>::Type TSequenceVIter;
typedef typename Value<TSequenceV>::Type TValueV;
//-------------------------------------------------------------------------
//variables
// TRow row0 = row(align_,0);
// TRow row1 = row(align_,1);
// beginPosition == # leading gaps
// endPosition == length of clipped region without trailing gaps
// clippedEndPosition == source position of clipping end.
// TAlignIterator ali_it0_stop = iter(row0,beginPosition(row0));
// TAlignIterator ali_it1_stop = iter(row1,beginPosition(row1));
TGapsHIter ali_it0_stop = begin(gapsH);
TGapsVIter ali_it1_stop = begin(gapsV);
// SEQAN_ASSERT( endPosition(row0)- beginPosition(row0) == endPosition(row1)- beginPosition(row1) );
// TAlignIterator ali_it0 = iter(row0,endPosition(row0));
// TAlignIterator ali_it1 = iter(row1,endPosition(row1));
TGapsHIter ali_it0 = end(gapsH);
TGapsVIter ali_it1 = end(gapsV);
// TStringIterator x_begin = begin(source(row0))-1;
// TStringIterator y_begin = begin(source(row1))-1;
// TStringIterator x_end = iter(source(row0),clippedEndPosition(row0))-1;
// TStringIterator y_end = iter(source(row1),clippedEndPosition(row1))-1;
TSequenceHIter x_begin = begin(source(gapsH))-1;
TSequenceVIter y_begin = begin(source(gapsV))-1;
TSequenceHIter x_end = iter(source(gapsH), endPosition(gapsH) - 1);
TSequenceVIter y_end = iter(source(gapsV), endPosition(gapsV) - 1);
// TStringIterator x = x_end;
// TStringIterator y = y_end;
// TStringIterator x_stop = x_end;
TSequenceHIter x = x_end;
TSequenceVIter y = y_end;
TSequenceHIter x_stop = x_end;
TScoreValue score_gap = scoreGapExtend(score_);
TScoreValue h,v;
TMatrixIterator finger0 = iter(sw.matrix,sw.bestBeginPos);
TMatrixIterator end_col = finger0;
TMatrixIterator finger1 = finger0;
TMatrixIterator forbidden = finger0;
bool different = true;
bool forbidden_reached = true;
bool end_reached = false;
bool skip_row = false;
/* int str0_length = length(source(row(align_,0)))+1;
int str1_length = length(source(row(align_,1)))+1;
for(int i = 0; i <str1_length; ++i){
for(int j=0;j<str0_length;++j){
std::cout << getValue(sw.matrix,(str0_length*i)+j);
if(sw.forbidden[(str0_length*i)+j]==true)
std::cout <<"(1) ";
else
std::cout <<"(0) ";
}
std::cout <<"\n";
}*/
clearClipping(gapsH);
clearClipping(gapsV);
// setClippedBeginPosition(row(align_, 0),0);
// setClippedBeginPosition(row(align_, 1),0);
for (y = y_end; (y != y_begin) ; --y)
{
different = true;
//compute next "forbidden" cell (where you are not allowed to go diagonal)
if(forbidden_reached && !end_reached)
{
if(ali_it0==ali_it0_stop && ali_it1==ali_it1_stop)
{
end_reached = true;
}
if(!end_reached)
{
--ali_it0;
goPrevious(forbidden,1);
while(isGap(ali_it0)&& ali_it0!=ali_it0_stop)
{
skip_row = true;
--ali_it0;
--ali_it1;
goPrevious(forbidden,1);
}
--ali_it1;
goPrevious(forbidden,0);
while(isGap(ali_it1)&& ali_it1!=ali_it1_stop)
{
--ali_it1;
--ali_it0;
goPrevious(forbidden,0);
}
// mark the forbidden cell
sw.forbidden[position(forbidden)]=true;
forbidden_reached = false;
}
}
TValueV cy = *y;
h = *end_col;
finger1 = end_col; //points to last column
goPrevious(end_col, 1); //points to this column
finger0 = end_col;
v = *finger0;
x = x_end;
// declump the current row until the cells in the declumped and
// the cells in the original matrix are the same (or the border is reached)
// indicated by bool different
while (x != x_begin && different)
{
goPrevious(finger0, 0);
goPrevious(finger1, 0);
TScoreValue currScore = score(score_, *x, cy);
if (*x == cy && !(sw.forbidden[position(finger0)]))
{
v = h + currScore;
h = *finger1;
}
else
{
TScoreValue s1;
if(finger0 == forbidden)
{
skip_row = false;
forbidden_reached = true;
s1 = 0;
}
else
{
if(sw.forbidden[position(finger0)]) s1 = 0;
else s1 = h + currScore;
}
h = *finger1;
TScoreValue s2 = score_gap + ((h > v) ? h : v);
v = (s1 > s2) ? s1 : s2;
if (v < 0) v = 0;
}
// value is the same as in the original matrix
if(*finger0==v)
{
//x_stop is as far as we have to go at least
if(x<x_stop)
{
different=false;
// x_stop=x;
}
else
{
// x_end indicates where to start in the next row
if(x==x_end && ((!forbidden_reached && !skip_row)||end_reached))
{
--x_end;
goPrevious(end_col, 0);
}
}
}
if(x<x_stop)// && different)
{
x_stop=x;
}
*finger0 = v;
--x;
}
if(x_end==x_begin)
break;
}
// cout <<"...declumped.\n";
}
// ----------------------------------------------------------------------------
// Function _smithWatermanTrace()
// ----------------------------------------------------------------------------
// Traceback.
template <typename TSourceH, typename TGapsSpecH,
typename TSourceV, typename TGapsSpecV,
typename TScoreValue, typename TScoreSpec,
unsigned DIMENSION>
typename Iterator<Matrix<TScoreValue, DIMENSION>, Standard >::Type
_smithWatermanTrace(Gaps<TSourceH, TGapsSpecH> & gapsH,
Gaps<TSourceV, TGapsSpecV> & gapsV,
typename LocalAlignmentFinder<TScoreValue>::TBoolMatrix & fb_matrix,
Iter< Matrix<TScoreValue, DIMENSION>, PositionIterator > source_,
Score<TScoreValue, TScoreSpec> const & scoring_) {
//typedefs
typedef Iter<Matrix<TScoreValue, DIMENSION>, PositionIterator > TMatrixIterator;
typedef typename Position<Matrix<TScoreValue, DIMENSION> >::Type TPosition;
// typedef Segment<TTargetSource, InfixSegment> TTargetSourceSegment;
typedef typename Iterator<TSourceH, Standard>::Type TSourceIteratorH;
typedef typename Iterator<TSourceV, Standard>::Type TSourceIteratorV;
typedef Gaps<TSourceH, TGapsSpecH> TGapsH;
typedef Gaps<TSourceV, TGapsSpecV> TGapsV;
typedef typename Iterator<TGapsH, Rooted>::Type TTargetIteratorH;
typedef typename Iterator<TGapsV, Rooted>::Type TTargetIteratorV;
//-------------------------------------------------------------------------
//variables
TPosition pos_0 = coordinate(source_, 0);
TPosition pos_1 = coordinate(source_, 1);
TSourceH strH = source(gapsH);
TSourceV strV = source(gapsV);
TTargetIteratorH target_0 = iter(gapsH, pos_0);
TTargetIteratorV target_1 = iter(gapsV, pos_1);
TSourceIteratorH it_0 = iter(strH, pos_0, Standard());
TSourceIteratorH it_0_end = end(strH);
TSourceIteratorV it_1 = iter(strV, pos_1, Standard());
TSourceIteratorV it_1_end = end(strV);
TScoreValue score_gap = scoreGapExtend(scoring_);
//-------------------------------------------------------------------------
//follow the trace until 0 is reached
while ((*source_!=0) && (it_0 != it_0_end) && (it_1 != it_1_end))
{
bool gv;
bool gh;
bool forbidden = fb_matrix[position(source_)];
if (*it_0 == *it_1 && !forbidden)
{
gv = gh = true;
}
else
{
TMatrixIterator it_ = source_;
goNext(it_, 0);
TScoreValue v = *it_ + score_gap;
TScoreValue d;
if(forbidden)
d = 0;
else{
goNext(it_, 1);
d = *it_ + score(scoring_, *it_0, *it_1);
}
it_ = source_;
goNext(it_, 1);
TScoreValue h = *it_ + score_gap;
gv = (v >= h) | (d >= h);
gh = (h > v) | (d >= v);
}
if (gv)
{
++it_0;
goNext(source_, 0);
}
else
{
insertGap(target_0);
}
if (gh)
{
++it_1;
goNext(source_, 1);
}
else
{
insertGap(target_1);
}
++target_0;
++target_1;
}
// We have removed all gaps and clippings from gapsH and gapsV in the calling functions, so the following works.
// Note that we have to set the end position first.
// TODO(holtgrew): Use setBegin/EndPosition().
setClippedEndPosition(gapsH, toViewPosition(gapsH, position(it_0, strH)));
setClippedEndPosition(gapsV, toViewPosition(gapsV, position(it_1, strV)));
setClippedBeginPosition(gapsH, toViewPosition(gapsH, pos_0));
setClippedBeginPosition(gapsV, toViewPosition(gapsV, pos_1));
return source_;
}
// ----------------------------------------------------------------------------
// Function _getNextBestEndPosition()
// ----------------------------------------------------------------------------
// Adjust the priority queue of scores until the true maximum is found.
template <typename TScoreValue>
typename LocalAlignmentFinder<TScoreValue>::TMatrixPosition
_getNextBestEndPosition(LocalAlignmentFinder<TScoreValue> & sw ,
TScoreValue cutoff)
{
// get maximal score from priority queue
TScoreValue topScore = 0;
if (!empty(sw.pQ))
topScore = getValue(sw.matrix, top(sw.pQ).id_);
// check if matrix entry of topScore did not change while declumping
if (!empty(sw.pQ)) {
while (top(sw.pQ).value_ != topScore) {
if (topScore >= cutoff) {
((sw.pQ).heap[0]).value_ = topScore;
adjustTop(sw.pQ);
} else {
pop(sw.pQ);
}
if (!empty(sw.pQ)) topScore = getValue(sw.matrix, top(sw.pQ).id_);
else break;
}
}
// priority queue with top scores is empty
if(empty(sw.pQ)) {//||top(sw.pQ).value_<cutoff) {
sw.needReinit = true;
return 0;
}
typename LocalAlignmentFinder<TScoreValue>::TMatrixPosition ret_pos = top(sw.pQ).id_;
sw.bestEndPos = ret_pos;
pop(sw.pQ);
return ret_pos;
}
// ----------------------------------------------------------------------------
// Function _smithWaterman()
// ----------------------------------------------------------------------------
// Wrapper that computes the matrix and does the backtracking for the best alignment
template <typename TSourceH, typename TGapsSpecH,
typename TSourceV, typename TGapsSpecV,
typename TScoreValue, typename TScoreSpec>
TScoreValue
_smithWaterman(Gaps<TSourceH, TGapsSpecH> & gapsH,
Gaps<TSourceV, TGapsSpecV> & gapsV,
LocalAlignmentFinder<TScoreValue> & sw_finder,
Score<TScoreValue, TScoreSpec> const & score_,
TScoreValue cutoff)
{
// TODO(holtgrew): This sourceSegment() stuff is confusing... Do we *really* need this?
// Clear gaps and clipping from result Gaps structures.
clearGaps(gapsH);
clearGaps(gapsV);
clearClipping(gapsH);
clearClipping(gapsV);
_initLocalAlignmentFinder(sourceSegment(gapsH), sourceSegment(gapsV), sw_finder, WatermanEggert());
TScoreValue ret = _smithWatermanGetMatrix(sw_finder, sourceSegment(gapsH), sourceSegment(gapsV), score_,cutoff);
if(ret==0)
return ret;
sw_finder.needReinit = false;
typedef Iter<typename LocalAlignmentFinder<TScoreValue>::TMatrix,PositionIterator > TMatrixIterator;
TMatrixIterator best_begin;
// TODO(holtgrew): What does the following comment mean?
// TODO: sw_finder statt kram
best_begin = _smithWatermanTrace(gapsH, gapsV, sw_finder.forbidden,iter(sw_finder.matrix,(top(sw_finder.pQ)).id_), score_);
sw_finder.bestBeginPos = position(best_begin);
pop(sw_finder.pQ);
return ret;
}
// ----------------------------------------------------------------------------
// Function _smithWatermanGetNext()
// ----------------------------------------------------------------------------
// Wrapper that declumps the matrix and traces back the next best alignment
template <typename TSequenceH, typename TGapsSpecH, typename TSequenceV, typename TGapsSpecV, typename TScoreValue, typename TScoreSpec>
TScoreValue
_smithWatermanGetNext(Gaps<TSequenceH, TGapsSpecH> & gapsH,
Gaps<TSequenceV, TGapsSpecV> & gapsV,
LocalAlignmentFinder<TScoreValue> & sw_finder ,
Score<TScoreValue, TScoreSpec> const & score_,
TScoreValue cutoff)
{
_smithWatermanDeclump(sw_finder, gapsH, gapsV, score_);
clearGaps(gapsH);
clearGaps(gapsV);
clearClipping(gapsH);
clearClipping(gapsV);
typename LocalAlignmentFinder<TScoreValue>::TMatrixPosition next_best_end;
next_best_end = _getNextBestEndPosition(sw_finder,cutoff);
if(next_best_end==0)
return 0;
typename LocalAlignmentFinder<TScoreValue>::TMatrixIterator next_best_begin;
next_best_begin= _smithWatermanTrace(gapsH, gapsV, sw_finder.forbidden,iter(sw_finder.matrix,next_best_end), score_);
sw_finder.bestBeginPos = position(next_best_begin);
return getValue(sw_finder.matrix,next_best_end);
}
} // namespace seqan
#endif // #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_LOCAL_ALIGNMENT_WATERMAN_EGGERT_IMPL_H_
