1 // ========================================================================== 2 // SeqAn - The Library for Sequence Analysis 3 // ========================================================================== 4 // Copyright (c) 2006-2015, Knut Reinert, FU Berlin 5 // All rights reserved. 6 // 7 // Redistribution and use in source and binary forms, with or without 8 // modification, are permitted provided that the following conditions are met: 9 // 10 // * Redistributions of source code must retain the above copyright 11 // notice, this list of conditions and the following disclaimer. 12 // * Redistributions in binary form must reproduce the above copyright 13 // notice, this list of conditions and the following disclaimer in the 14 // documentation and/or other materials provided with the distribution. 15 // * Neither the name of Knut Reinert or the FU Berlin nor the names of 16 // its contributors may be used to endorse or promote products derived 17 // from this software without specific prior written permission. 18 // 19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 20 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 21 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 22 // ARE DISCLAIMED. IN NO EVENT SHALL KNUT REINERT OR THE FU BERLIN BE LIABLE 23 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 24 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 25 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 26 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 27 // LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 28 // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH 29 // DAMAGE. 30 // 31 // ========================================================================== 32 // Author: Rene Rahn <rene.rahn@fu-berlin.de> 33 // ========================================================================== 34 // Implements the traceback algorithm. 35 // ========================================================================== 36 37 #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_DP_TRACEBACK_IMPL_H_ 38 #define SEQAN_INCLUDE_SEQAN_ALIGN_DP_TRACEBACK_IMPL_H_ 39 40 // TODO(holtgrew): GapsRight traceback is currently untested. 41 // TODO(rmaerker): Change Tracback to TraceConfig<TGapsPlacement, TPathSpec> | TraceBackOff 42 43 namespace seqan { 44 45 // ============================================================================ 46 // Forwards 47 // ============================================================================ 48 49 // ============================================================================ 50 // Tags, Classes, Enums 51 // ============================================================================ 52 53 // ---------------------------------------------------------------------------- 54 // Class TracebackCoordinator_ 55 // ---------------------------------------------------------------------------- 56 57 template <typename TPosition> 58 class TracebackCoordinator_ 59 { 60 public: 61 TPosition _currColumn; 62 TPosition _currRow; 63 TPosition _endColumn; 64 TPosition _endRow; 65 TPosition _breakpoint1; // First breakpoint where banded trace switches to unbanded trace. 66 TPosition _breakpoint2; // Second breakpoint where unbanded trace switches back to banded trace. Only if begin of upper diagonal is bigger than end of lower diagonal. 67 bool _isInBand; 68 69 template <typename TBandFlag, typename TSizeH, typename TSizeV> TracebackCoordinator_(TPosition currColumn,TPosition currRow,DPBandConfig<TBandFlag> const & band,TSizeH seqHSize,TSizeV seqVSize)70 TracebackCoordinator_(TPosition currColumn, 71 TPosition currRow, 72 DPBandConfig<TBandFlag> const & band, 73 TSizeH seqHSize, 74 TSizeV seqVSize) 75 : _currColumn(currColumn), 76 _currRow(currRow), 77 _endColumn(0u), 78 _endRow(0u), 79 _breakpoint1(0u), 80 _breakpoint2(0u), 81 _isInBand(false) 82 { 83 _initTracebackCoordinator(*this, band, seqHSize, seqVSize); 84 } 85 86 template <typename TBandFlag, typename TSizeH, typename TSizeV> TracebackCoordinator_(TPosition currColumn,TPosition currRow,TPosition endColumn,TPosition endRow,DPBandConfig<TBandFlag> const & band,TSizeH seqHSize,TSizeV seqVSize)87 TracebackCoordinator_(TPosition currColumn, 88 TPosition currRow, 89 TPosition endColumn, 90 TPosition endRow, 91 DPBandConfig<TBandFlag> const & band, 92 TSizeH seqHSize, 93 TSizeV seqVSize) 94 : _currColumn(currColumn), 95 _currRow(currRow), 96 _endColumn(endColumn), 97 _endRow(endRow), 98 _breakpoint1(0u), 99 _breakpoint2(0u), 100 _isInBand(false) 101 { 102 _initTracebackCoordinator(*this, band, seqHSize, seqVSize); 103 } 104 }; 105 106 // ============================================================================ 107 // Metafunctions 108 // ============================================================================ 109 110 // ---------------------------------------------------------------------------- 111 // Metafunction PreferGapsAtEnd_ 112 // ---------------------------------------------------------------------------- 113 114 // Checks whether the gaps at the end should be preferred over a matching area. 115 template <typename TDPProfile> 116 struct PreferGapsAtEnd_ : False{}; 117 118 template <typename TAlgorithm, typename TTracebackSpec> 119 struct PreferGapsAtEnd_<DPProfile_<TAlgorithm, AffineGaps, TTracebackSpec > > : True{}; 120 121 template <typename TAlgorithm, typename TTraceSpec> 122 struct PreferGapsAtEnd_<DPProfile_<TAlgorithm, LinearGaps, TracebackOn<TracebackConfig_<TTraceSpec, GapsRight> > > > : True{}; 123 124 125 // ============================================================================ 126 // Functions 127 // ============================================================================ 128 129 // ---------------------------------------------------------------------------- 130 // Function _hasReachedEnd() 131 // ---------------------------------------------------------------------------- 132 133 template <typename TPosition> 134 inline bool 135 _hasReachedEnd(TracebackCoordinator_<TPosition> const & coordinator) 136 { 137 return coordinator._currColumn <= coordinator._endColumn || coordinator._currRow <= coordinator._endRow; 138 } 139 140 // ---------------------------------------------------------------------------- 141 // Function _initTracebackCoordinator() 142 // ---------------------------------------------------------------------------- 143 144 template <typename TPosition, typename TBandFlag, typename TSizeH, typename TSizeV> 145 inline void 146 _initTracebackCoordinator(TracebackCoordinator_<TPosition> & coordinator, 147 DPBandConfig<TBandFlag> const & band, 148 TSizeH seqHSize, 149 TSizeV seqVSize) 150 { 151 typedef typename Position<DPBandConfig<TBandFlag> >::Type TBandPosition; 152 if (IsSameType<TBandFlag, BandOn>::VALUE) 153 { 154 // Adapt the current column value when the lower diagonal is positive (shift right in horizontal direction). 155 if (lowerDiagonal(band) >= 0) 156 coordinator._currColumn += static_cast<TPosition>(lowerDiagonal(band)); 157 // Adapt the current row value when the current column comes after the upper diagonal (shift down in vertical direction). 158 if (static_cast<TBandPosition>(coordinator._currColumn) > upperDiagonal(band)) 159 coordinator._currRow += coordinator._currColumn - upperDiagonal(band); 160 // Adapt the end row value when the end column comes after the upper diagonal (shift down in vertical direction). 161 if (static_cast<TBandPosition>(coordinator._endColumn) > upperDiagonal(band)) 162 coordinator._endRow += coordinator._endColumn - upperDiagonal(band); 163 164 coordinator._breakpoint1 = _min(seqHSize, static_cast<TSizeH>(_max(0, upperDiagonal(band)))); 165 coordinator._breakpoint2 = _min(seqHSize, static_cast<TSizeH>(_max(0, static_cast<TBandPosition>(seqVSize) + 166 lowerDiagonal(band)))); 167 // Update the current row if the current column is before the upper diagoal or the first column where the maximal band size is reached. 168 if (coordinator._currColumn < _min(coordinator._breakpoint1, coordinator._breakpoint2)) 169 coordinator._currRow -= _min(coordinator._breakpoint1, coordinator._breakpoint2) - coordinator._currColumn; 170 coordinator._isInBand = true; 171 } 172 } 173 174 // ---------------------------------------------------------------------------- 175 // Function _isInBand() 176 // ---------------------------------------------------------------------------- 177 178 template <typename TPosition> 179 inline bool 180 _isInBand(TracebackCoordinator_<TPosition> const & coordinator) 181 { 182 if (!coordinator._isInBand) 183 return coordinator._isInBand; 184 return (coordinator._currColumn > coordinator._breakpoint1 || coordinator._currColumn <= coordinator._breakpoint2); 185 } 186 187 188 // ---------------------------------------------------------------------------- 189 // Function _doTracebackGoDiagonal() 190 // ---------------------------------------------------------------------------- 191 192 template <typename TTarget, typename TDPTraceMatrixNavigator, typename TTraceValue, typename TSize, typename TPosition, 193 typename TGapCosts> 194 inline void 195 _doTracebackGoDiagonal(TTarget & target, 196 TDPTraceMatrixNavigator & matrixNavigator, 197 TTraceValue & traceValue, 198 TTraceValue & lastTraceValue, 199 TSize & fragmentLength, 200 TracebackCoordinator_<TPosition> & tracebackCoordinator, 201 TGapCosts const &) 202 { 203 if (!(lastTraceValue & TraceBitMap_::DIAGONAL)) // the old trace value was not diagonal 204 { 205 _recordSegment(target, tracebackCoordinator._currColumn, tracebackCoordinator._currRow, fragmentLength, 206 lastTraceValue); 207 208 lastTraceValue = TraceBitMap_::DIAGONAL; 209 fragmentLength = 0; 210 } 211 _traceDiagonal(matrixNavigator, _isInBand(tracebackCoordinator)); 212 traceValue = value(matrixNavigator); 213 --tracebackCoordinator._currColumn; 214 --tracebackCoordinator._currRow; 215 ++fragmentLength; 216 } 217 218 // ---------------------------------------------------------------------------- 219 // Function _doTracebackGoVertical() 220 // ---------------------------------------------------------------------------- 221 222 template <typename TTarget, typename TDPTraceMatrixNavigator, typename TTraceValue, typename TSize, typename TPosition, 223 typename TGapCosts> 224 inline void 225 _doTracebackGoVertical(TTarget & target, 226 TDPTraceMatrixNavigator & matrixNavigator, 227 TTraceValue & traceValue, 228 TTraceValue & lastTraceValue, 229 TSize & fragmentLength, 230 TracebackCoordinator_<TPosition> & tracebackCoordinator, 231 TGapCosts const &) 232 { 233 if (!(lastTraceValue & TraceBitMap_::VERTICAL)) // the old trace value was not diagonal 234 { 235 _recordSegment(target, tracebackCoordinator._currColumn, tracebackCoordinator._currRow, fragmentLength, 236 lastTraceValue); 237 238 lastTraceValue = TraceBitMap_::VERTICAL; 239 fragmentLength = 0; 240 } 241 // We are in a vertical gap. So continue after we reach the end of the vertical gap. 242 if (IsSameType<TGapCosts, AffineGaps>::VALUE) 243 { 244 while ((!(traceValue & TraceBitMap_::VERTICAL_OPEN) || (traceValue & TraceBitMap_::VERTICAL)) && (tracebackCoordinator._currRow != 1)) 245 { 246 _traceVertical(matrixNavigator, _isInBand(tracebackCoordinator)); 247 traceValue = value(matrixNavigator); 248 --tracebackCoordinator._currRow; 249 ++fragmentLength; 250 } 251 // We have to ensure, that we do not continue in vertical direction if we reached a vertical_open sign. 252 _traceVertical(matrixNavigator, _isInBand(tracebackCoordinator)); 253 // Forbid continuing in vertical direction. 254 traceValue = value(matrixNavigator); // & (TraceBitMap_::NO_VERTICAL_TRACEBACK | TraceBitMap_::MAX_FROM_HORIZONTAL_MATRIX); 255 --tracebackCoordinator._currRow; 256 ++fragmentLength; 257 } 258 else 259 { 260 _traceVertical(matrixNavigator, _isInBand(tracebackCoordinator)); 261 traceValue = value(matrixNavigator); 262 --tracebackCoordinator._currRow; 263 ++fragmentLength; 264 } 265 } 266 267 // ---------------------------------------------------------------------------- 268 // Function _doTracebackMaxFromVertical() 269 // ---------------------------------------------------------------------------- 270 271 template <typename TTarget, typename TDPTraceMatrixNavigator, typename TTraceValue, typename TSize, typename TPosition, 272 typename TGapCosts> 273 inline void 274 _doTracebackMaxFromVertical(TTarget & target, 275 TDPTraceMatrixNavigator & matrixNavigator, 276 TTraceValue & traceValue, 277 TTraceValue & lastTraceValue, 278 TSize & fragmentLength, 279 TracebackCoordinator_<TPosition> & tracebackCoordinator, 280 TGapCosts const &) 281 { 282 if (!(lastTraceValue & TraceBitMap_::VERTICAL)) // the old trace value was not diagonal 283 { 284 _recordSegment(target, tracebackCoordinator._currColumn, tracebackCoordinator._currRow, fragmentLength, 285 lastTraceValue); 286 lastTraceValue = TraceBitMap_::VERTICAL; 287 fragmentLength = 0; 288 } 289 _traceVertical(matrixNavigator, _isInBand(tracebackCoordinator)); 290 // Forbid continuing in vertical direction. 291 traceValue = value(matrixNavigator); // & (TraceBitMap_::NO_VERTICAL_TRACEBACK | TraceBitMap_::MAX_FROM_HORIZONTAL_MATRIX); 292 --tracebackCoordinator._currRow; 293 ++fragmentLength; 294 } 295 296 // ---------------------------------------------------------------------------- 297 // Function _doTracebackGoHorizontal() 298 // ---------------------------------------------------------------------------- 299 300 template <typename TTarget, typename TDPTraceMatrixNavigator, typename TTraceValue, typename TSize, typename TPosition, 301 typename TGapCosts> 302 inline void 303 _doTracebackGoHorizontal(TTarget & target, 304 TDPTraceMatrixNavigator & matrixNavigator, 305 TTraceValue & traceValue, 306 TTraceValue & lastTraceValue, 307 TSize & fragmentLength, 308 TracebackCoordinator_<TPosition> & tracebackCoordinator, 309 TGapCosts const &) 310 { 311 if (!(lastTraceValue & TraceBitMap_::HORIZONTAL)) // the old trace value was not diagonal 312 { 313 _recordSegment(target, tracebackCoordinator._currColumn, tracebackCoordinator._currRow, fragmentLength, 314 lastTraceValue); 315 316 lastTraceValue = TraceBitMap_::HORIZONTAL; 317 fragmentLength = 0; 318 } 319 if (IsSameType<TGapCosts, AffineGaps>::VALUE) 320 { 321 while ((!(traceValue & TraceBitMap_::HORIZONTAL_OPEN) || (traceValue & TraceBitMap_::HORIZONTAL)) && (tracebackCoordinator._currColumn != 1)) 322 { 323 _traceHorizontal(matrixNavigator, _isInBand(tracebackCoordinator)); 324 traceValue = value(matrixNavigator); 325 --tracebackCoordinator._currColumn; 326 ++fragmentLength; 327 } 328 _traceHorizontal(matrixNavigator, _isInBand(tracebackCoordinator)); 329 // Forbid continuing in horizontal direction. 330 traceValue = value(matrixNavigator); // & (TraceBitMap_::NO_HORIZONTAL_TRACEBACK | TraceBitMap_::MAX_FROM_VERTICAL_MATRIX); 331 --tracebackCoordinator._currColumn; 332 ++fragmentLength; 333 } 334 else 335 { 336 _traceHorizontal(matrixNavigator, _isInBand(tracebackCoordinator)); 337 traceValue = value(matrixNavigator); 338 --tracebackCoordinator._currColumn; 339 ++fragmentLength; 340 } 341 } 342 343 // ---------------------------------------------------------------------------- 344 // Function _doTracebackMaxFromHorizontal() 345 // ---------------------------------------------------------------------------- 346 347 template <typename TTarget, typename TDPTraceMatrixNavigator, typename TTraceValue, typename TSize, typename TPosition, 348 typename TGapCosts> 349 inline void 350 _doTracebackMaxFromHorizontal(TTarget & target, 351 TDPTraceMatrixNavigator & matrixNavigator, 352 TTraceValue & traceValue, 353 TTraceValue & lastTraceValue, 354 TSize & fragmentLength, 355 TracebackCoordinator_<TPosition> & tracebackCoordinator, 356 TGapCosts const &) 357 { 358 if (!(lastTraceValue & TraceBitMap_::HORIZONTAL)) // the old trace value was not diagonal 359 { 360 _recordSegment(target, tracebackCoordinator._currColumn, tracebackCoordinator._currRow, fragmentLength, 361 lastTraceValue); 362 lastTraceValue = TraceBitMap_::HORIZONTAL; 363 fragmentLength = 0; 364 } 365 _traceHorizontal(matrixNavigator, _isInBand(tracebackCoordinator)); 366 // Forbid continuing in horizontal direction. 367 traceValue = value(matrixNavigator); // & (TraceBitMap_::NO_HORIZONTAL_TRACEBACK | TraceBitMap_::MAX_FROM_VERTICAL_MATRIX); 368 --tracebackCoordinator._currColumn; 369 ++fragmentLength; 370 } 371 372 // ---------------------------------------------------------------------------- 373 // Function _doTraceback() 374 // ---------------------------------------------------------------------------- 375 376 template <typename TTarget, typename TDPTraceMatrixNavigator, typename TTraceValue, typename TSize, typename TPosition, 377 typename TGapCosts, typename TIsGapsLeft> 378 inline void 379 _doTraceback(TTarget & target, 380 TDPTraceMatrixNavigator & matrixNavigator, 381 TTraceValue & traceValue, 382 TTraceValue & lastTraceValue, 383 TSize & fragmentLength, 384 TracebackCoordinator_<TPosition> & tracebackCoordinator, 385 TGapCosts const & gapsCost, 386 TIsGapsLeft const & /*isGapsLeft*/) 387 { 388 if (TIsGapsLeft::VALUE) // Gaps should be placed on the left. 389 { 390 if (traceValue & TraceBitMap_::DIAGONAL) 391 { 392 _doTracebackGoDiagonal(target, matrixNavigator, traceValue, lastTraceValue, fragmentLength, tracebackCoordinator, gapsCost); 393 } // In case of Gotoh we prefer the longest possible way in this direction. 394 else if (traceValue & TraceBitMap_::MAX_FROM_VERTICAL_MATRIX && traceValue & TraceBitMap_::VERTICAL) 395 { 396 _doTracebackGoVertical(target, matrixNavigator, traceValue, lastTraceValue, fragmentLength, tracebackCoordinator, gapsCost); 397 } 398 else if (traceValue & TraceBitMap_::MAX_FROM_VERTICAL_MATRIX && traceValue & TraceBitMap_::VERTICAL_OPEN) 399 { 400 _doTracebackMaxFromVertical(target, matrixNavigator, traceValue, lastTraceValue, fragmentLength, tracebackCoordinator, gapsCost); 401 } 402 else if (traceValue & TraceBitMap_::MAX_FROM_HORIZONTAL_MATRIX && traceValue & TraceBitMap_::HORIZONTAL) 403 { 404 _doTracebackGoHorizontal(target, matrixNavigator, traceValue, lastTraceValue, fragmentLength, tracebackCoordinator, gapsCost); 405 } 406 else if (traceValue & TraceBitMap_::MAX_FROM_HORIZONTAL_MATRIX && traceValue & TraceBitMap_::HORIZONTAL_OPEN) 407 { 408 _doTracebackMaxFromHorizontal(target, matrixNavigator, traceValue, lastTraceValue, fragmentLength, tracebackCoordinator, gapsCost); 409 } // In case of Gotoh we prefer the longest possible way in this direction. 410 else // the trace back is either NONE or something else 411 { 412 if (traceValue == TraceBitMap_::NONE) 413 { 414 return; 415 } 416 SEQAN_ASSERT_FAIL("Reached undefined traceback value!"); 417 } 418 } 419 else // Gaps should be placed on the right. 420 { 421 if (traceValue & TraceBitMap_::MAX_FROM_VERTICAL_MATRIX && traceValue & TraceBitMap_::VERTICAL) 422 { 423 _doTracebackGoVertical(target, matrixNavigator, traceValue, lastTraceValue, fragmentLength, tracebackCoordinator, gapsCost); 424 } 425 else if (traceValue & TraceBitMap_::MAX_FROM_VERTICAL_MATRIX && traceValue & TraceBitMap_::VERTICAL_OPEN) 426 { 427 _doTracebackMaxFromVertical(target, matrixNavigator, traceValue, lastTraceValue, fragmentLength, tracebackCoordinator, gapsCost); 428 } 429 else if (traceValue & TraceBitMap_::MAX_FROM_HORIZONTAL_MATRIX && traceValue & TraceBitMap_::HORIZONTAL) 430 { 431 _doTracebackGoHorizontal(target, matrixNavigator, traceValue, lastTraceValue, fragmentLength, tracebackCoordinator, gapsCost); 432 } 433 else if (traceValue & TraceBitMap_::MAX_FROM_HORIZONTAL_MATRIX && traceValue & TraceBitMap_::HORIZONTAL_OPEN) 434 { 435 _doTracebackMaxFromHorizontal(target, matrixNavigator, traceValue, lastTraceValue, fragmentLength, tracebackCoordinator, gapsCost); 436 } // In case of Gotoh we prefer the longest possible way in this direction. 437 else if (traceValue & TraceBitMap_::DIAGONAL) 438 { 439 _doTracebackGoDiagonal(target, matrixNavigator, traceValue, lastTraceValue, fragmentLength, tracebackCoordinator, gapsCost); 440 } // In case of Gotoh we prefer the longest possible way in this direction. 441 else // the trace back is either NONE or something else 442 { 443 if (traceValue == TraceBitMap_::NONE) 444 { 445 return; 446 } 447 SEQAN_ASSERT_FAIL("Reached undefined traceback value!"); 448 } 449 } 450 } 451 452 // ---------------------------------------------------------------------------- 453 // Function _retrieveInitialTraceDirection() 454 // ---------------------------------------------------------------------------- 455 456 template <typename TTraceValue, typename TDPProfile> 457 inline TTraceValue 458 _retrieveInitialTraceDirection(TTraceValue & traceValue, TDPProfile const & /*dpProfile*/) 459 { 460 if (PreferGapsAtEnd_<TDPProfile>::VALUE) 461 { 462 if (traceValue & TraceBitMap_::MAX_FROM_VERTICAL_MATRIX) 463 { 464 traceValue &= (TraceBitMap_::VERTICAL | TraceBitMap_::VERTICAL_OPEN | TraceBitMap_::MAX_FROM_VERTICAL_MATRIX); 465 return TraceBitMap_::VERTICAL; 466 } 467 else if (traceValue & TraceBitMap_::MAX_FROM_HORIZONTAL_MATRIX) 468 { 469 traceValue &= (TraceBitMap_::HORIZONTAL | TraceBitMap_::HORIZONTAL_OPEN | TraceBitMap_::MAX_FROM_HORIZONTAL_MATRIX); 470 return TraceBitMap_::HORIZONTAL; 471 } 472 return TraceBitMap_::DIAGONAL; // We set the last value to the 473 } 474 475 if (traceValue & TraceBitMap_::DIAGONAL) 476 return TraceBitMap_::DIAGONAL; 477 if (traceValue & (TraceBitMap_::VERTICAL | TraceBitMap_::MAX_FROM_VERTICAL_MATRIX)) 478 return TraceBitMap_::VERTICAL; 479 if (traceValue & (TraceBitMap_::HORIZONTAL | TraceBitMap_::MAX_FROM_HORIZONTAL_MATRIX)) 480 return TraceBitMap_::HORIZONTAL; 481 482 return TraceBitMap_::NONE; 483 } 484 485 // ---------------------------------------------------------------------------- 486 // Function _computeTraceback() 487 // ---------------------------------------------------------------------------- 488 489 template <typename TTarget, typename TDPTraceMatrixNavigator, typename TSequenceH, typename TSequenceV, 490 typename TBandFlag, typename TAlgorithm, typename TGapCosts, typename TTracebackSpec> 491 void _computeTraceback(TTarget & target, 492 TDPTraceMatrixNavigator & matrixNavigator, 493 unsigned maxHostPosition, 494 TSequenceH const & seqH, 495 TSequenceV const & seqV, 496 DPBandConfig<TBandFlag> const & band, 497 DPProfile_<TAlgorithm, TGapCosts, TTracebackSpec> const & dpProfile) 498 { 499 typedef typename Container<TDPTraceMatrixNavigator>::Type TContainer; 500 typedef typename Size<TContainer>::Type TSize; 501 typedef typename Position<TContainer>::Type TPosition; 502 typedef typename TraceBitMap_::TTraceValue TTraceValue; 503 typedef typename Size<TSequenceH>::Type TSizeH; 504 typedef typename Size<TSequenceV>::Type TSizeV; 505 506 if (IsSameType<TTracebackSpec, TracebackOff>::VALUE) 507 return; 508 509 // Determine whether or not we place gaps to the left. 510 typedef typename IsGapsLeft_<TTracebackSpec>::Type TIsGapsLeft; 511 512 TSizeH seqHSize = length(seqH); 513 TSizeV seqVSize = length(seqV); 514 515 // Set the navigator to the position where the maximum was found. 516 _setToPosition(matrixNavigator, maxHostPosition); 517 518 SEQAN_ASSERT_LEQ(coordinate(matrixNavigator, +DPMatrixDimension_::HORIZONTAL), seqHSize); 519 SEQAN_ASSERT_LEQ(coordinate(matrixNavigator, +DPMatrixDimension_::VERTICAL), seqVSize); 520 521 TTraceValue traceValue = value(matrixNavigator); 522 TTraceValue lastTraceValue = _retrieveInitialTraceDirection(traceValue, dpProfile); 523 524 TracebackCoordinator_<TPosition> tracebackCoordinator(coordinate(matrixNavigator, +DPMatrixDimension_::HORIZONTAL), 525 coordinate(matrixNavigator, +DPMatrixDimension_::VERTICAL), 526 band, seqHSize, seqVSize); 527 528 if (TraceTail_<TAlgorithm>::VALUE) 529 { 530 if (tracebackCoordinator._currRow != seqVSize) 531 _recordSegment(target, seqHSize, tracebackCoordinator._currRow, seqVSize - tracebackCoordinator._currRow, 532 +TraceBitMap_::VERTICAL); 533 if (tracebackCoordinator._currColumn != seqHSize) 534 _recordSegment(target, tracebackCoordinator._currColumn, tracebackCoordinator._currRow, seqHSize - 535 tracebackCoordinator._currColumn, +TraceBitMap_::HORIZONTAL); 536 } 537 538 TSize fragmentLength = 0; 539 while (!_hasReachedEnd(tracebackCoordinator) && traceValue != TraceBitMap_::NONE) 540 _doTraceback(target, matrixNavigator, traceValue, lastTraceValue, fragmentLength, tracebackCoordinator, TGapCosts(), TIsGapsLeft()); 541 542 543 // Record last detected fragment. 544 _recordSegment(target, tracebackCoordinator._currColumn, tracebackCoordinator._currRow, fragmentLength, lastTraceValue); 545 if (TraceHead_<TAlgorithm>::VALUE) 546 { 547 // Record leading gaps if any. 548 if (tracebackCoordinator._currRow != 0u) 549 _recordSegment(target, 0, 0, tracebackCoordinator._currRow, +TraceBitMap_::VERTICAL); 550 if (tracebackCoordinator._currColumn != 0u) 551 _recordSegment(target, 0, 0, tracebackCoordinator._currColumn, +TraceBitMap_::HORIZONTAL); 552 } 553 } 554 555 // Needed as a delegation method to allow invocation of both methods with host position and dpScout. 556 template <typename TTarget, typename TDPTraceMatrixNavigator, typename TDPCell, typename TScoutSpec, 557 typename TSequenceH, typename TSequenceV, typename TBandFlag, typename TAlgorithm, typename TGapCosts, 558 typename TTracebackSpec> 559 void _computeTraceback(TTarget & target, 560 TDPTraceMatrixNavigator & matrixNavigator, 561 DPScout_<TDPCell, TScoutSpec> const & dpScout, 562 TSequenceH const & seqH, 563 TSequenceV const & seqV, 564 DPBandConfig<TBandFlag> const & band, 565 DPProfile_<TAlgorithm, TGapCosts, TTracebackSpec> const & dpProfile) 566 { 567 _computeTraceback(target, matrixNavigator, maxHostPosition(dpScout), seqH, seqV, band, dpProfile); 568 } 569 570 } // namespace seqan 571 572 #endif // #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_DP_TRACEBACK_IMPL_H_ 573