1 // ========================================================================== 2 // SeqAn - The Library for Sequence Analysis 3 // ========================================================================== 4 // Copyright (c) 2006-2018, 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: David Weese <david.weese@fu-berlin.de> 33 // ========================================================================== 34 35 #ifndef INCLUDE_SEQAN_BAM_IO_CIGAR_H_ 36 #define INCLUDE_SEQAN_BAM_IO_CIGAR_H_ 37 38 namespace seqan { 39 40 // ============================================================================ 41 // Forwards 42 // ============================================================================ 43 44 // ============================================================================ 45 // Tags, Classes, Enums 46 // ============================================================================ 47 48 // ---------------------------------------------------------------------------- 49 // struct CigarElement 50 // ---------------------------------------------------------------------------- 51 52 /*! 53 * @class CigarElement 54 * @headerfile <seqan/bam_io.h> 55 * @brief One entry of a CIGAR string. 56 * 57 * @signature template <[typename TOperation[, typename TCount]]> 58 * class CigarElement; 59 * 60 * @tparam TOperation Type to use for storing operations, defaults to <tt>char</tt>. 61 * @tparam TCount Type to use for storing counts, defaults to <tt>unsigned</tt>. 62 */ 63 64 /*! 65 * @fn CigarElement::CigarElement 66 * @brief Constructor 67 * 68 * @signature CigarElement::CigarElement(); 69 * @signature CigarElement::CigarElement(operation, count); 70 * 71 * @param[in] operation The operation to use, of type <tt>TOperation</tt>. 72 * @param[in] count The operation count, of type <tt>TCount</tt>. 73 * 74 * @section Remarks 75 * 76 * The default constructor initialized both @link CigarElement::operation @endlink and @link CigarElement::count 77 * @endlink with <tt>0</tt>. 78 */ 79 80 /*! 81 * @var TCount CigarElement::count; 82 * 83 * @brief The number of operations. 84 */ 85 86 /*! 87 * @var TOperation CigarElement::operation; 88 * 89 * @brief The described operation. 90 */ 91 92 template <typename TOperation_ = char, typename TCount_ = unsigned> 93 struct CigarElement 94 { 95 typedef TOperation_ TOperation; 96 typedef TCount_ TCount; 97 98 TOperation operation; 99 TCount count; 100 CigarElementCigarElement101 CigarElement() : operation(0), count(0) {} 102 CigarElementCigarElement103 CigarElement(TOperation o, TCount c): 104 operation(o), 105 count(c) {} 106 }; 107 108 // ============================================================================ 109 // Metafunctions 110 // ============================================================================ 111 112 template <typename TOperation, typename TCount> 113 struct Size<CigarElement<TOperation, TCount> > 114 { 115 typedef TCount Type; 116 }; 117 118 // ============================================================================ 119 // Functions 120 // ============================================================================ 121 122 template <typename TOperation, typename TCount> 123 inline bool operator>(CigarElement<TOperation, TCount> const & lhs, 124 CigarElement<TOperation, TCount> const & rhs) 125 { 126 return lhs.operation > rhs.operation || (lhs.operation == rhs.operation && (lhs.count) > (rhs.count)); 127 } 128 129 template <typename TOperation, typename TCount> 130 inline bool operator<(CigarElement<TOperation, TCount> const & lhs, 131 CigarElement<TOperation, TCount> const & rhs) 132 { 133 return lhs.operation < rhs.operation || (lhs.operation == rhs.operation && (lhs.count) < (rhs.count)); 134 } 135 136 template <typename TOperation, typename TCount> 137 inline bool operator==(CigarElement<TOperation, TCount> const & lhs, 138 CigarElement<TOperation, TCount> const & rhs) 139 { 140 return lhs.operation == rhs.operation && (lhs.count) == (rhs.count); 141 } 142 143 // ---------------------------------------------------------------------------- 144 // toBamCigarElement() 145 // ---------------------------------------------------------------------------- 146 147 template <typename TOperation, typename TCount> 148 [[deprecated("The behavior is not clear. This function should not be used anymore.")]] 149 uint32_t toBamCigarElement(CigarElement<TOperation, TCount> const & cigarElement) 150 { 151 char operation = 0; 152 switch (cigarElement.operation) { 153 case 'X': operation += 1; SEQAN_FALLTHROUGH 154 case '=': operation += 1; SEQAN_FALLTHROUGH 155 case 'P': operation += 1; SEQAN_FALLTHROUGH 156 case 'H': operation += 1; SEQAN_FALLTHROUGH 157 case 'S': operation += 1; SEQAN_FALLTHROUGH 158 case 'N': operation += 1; SEQAN_FALLTHROUGH 159 case 'D': operation += 1; SEQAN_FALLTHROUGH 160 case 'I': operation += 1; SEQAN_FALLTHROUGH 161 case 'M': break; 162 } 163 return (cigarElement.count << 4) | operation; 164 } 165 166 // ---------------------------------------------------------------------------- 167 // getMDString() 168 // ---------------------------------------------------------------------------- 169 170 template < 171 typename TMDString, 172 typename TGaps1, 173 typename TGaps2> 174 inline unsigned 175 getMDString( 176 TMDString &md, 177 TGaps1 &gaps1, // typically reference 178 TGaps2 &gaps2) // typically read 179 { 180 typedef typename Value<TMDString>::Type TMDChar; 181 typedef typename Value<typename Host<TGaps1>::Type>::Type TVal1; 182 typedef typename Value<typename Host<TGaps2>::Type>::Type TVal2; 183 184 typename Iterator<TGaps1>::Type it1 = begin(gaps1); 185 typename Iterator<TGaps2>::Type it2 = begin(gaps2); 186 char op, lastOp = ' '; 187 unsigned numOps = 0; 188 unsigned errors = 0; 189 190 clear(md); 191 for (; !atEnd(it1) && !atEnd(it2); goNext(it1), goNext(it2)) 192 { 193 if (isGap(it1)) 194 { 195 if (!isGap(it2)) 196 ++errors; 197 continue; // insertion to the reference (gaps1) 198 // op = 'I'; // ignore insertions completely 199 } 200 if (isGap(it2)) 201 { 202 ++errors; 203 // if (op == 'I') // ignore paddings 204 // continue; 205 op = 'D'; // deletion from the reference (gaps1) 206 } 207 else 208 { 209 if ((TVal1)*it1 == (TVal2)*it2) 210 { 211 op = 'M'; 212 } 213 else 214 { 215 op = 'R'; 216 ++errors; 217 } 218 } 219 220 // append match run 221 if (lastOp != op) 222 { 223 if (lastOp == 'M') 224 { 225 std::stringstream num; 226 num << numOps; 227 append(md, num.str()); 228 } 229 numOps = 0; 230 } 231 232 // append deleted/replaced reference character 233 if (op != 'M') 234 { 235 // add ^ for deleted reference bases (from non-deletion to deletion) 236 if (op == 'D' && lastOp != 'D') 237 appendValue(md, '^'); 238 // add 0 for each replaced base that doesn't follow a match (for samtools/BWA compatibility) 239 else if (op == 'R' && lastOp != 'M') 240 appendValue(md, '0'); 241 appendValue(md, convert<TMDChar>(*it1)); 242 } 243 244 lastOp = op; 245 ++numOps; 246 } 247 SEQAN_ASSERT_EQ(atEnd(it1), atEnd(it2)); 248 if (lastOp == 'M') 249 { 250 std::stringstream num; 251 num << numOps; 252 append(md, num.str()); 253 } 254 return errors; 255 } 256 257 // ---------------------------------------------------------------------------- 258 // getCigarString() 259 // ---------------------------------------------------------------------------- 260 261 template < 262 typename TCigar, 263 typename TGaps1, 264 typename TGaps2, 265 typename TThresh> 266 inline void 267 getCigarString( 268 TCigar &cigar, 269 TGaps1 &gaps1, // typically reference 270 TGaps2 &gaps2, // typically read 271 TThresh splicedGapThresh) 272 { 273 typename Iterator<TGaps1>::Type it1 = begin(gaps1); 274 typename Iterator<TGaps2>::Type it2 = begin(gaps2); 275 // typedef typename Value<typename Host<TGaps1>::Type>::Type TVal1; 276 // typedef typename Value<typename Host<TGaps2>::Type>::Type TVal2; 277 278 clear(cigar); 279 char op, lastOp = ' '; 280 unsigned numOps = 0; 281 282 // std::cout << "gaps1\t" << gaps1 << std::endl; 283 // std::cout << "gaps2\t" << gaps2 << "\t" << clippedBeginPosition(gaps2) << std::endl; 284 for (; !atEnd(it1) && !atEnd(it2); goNext(it1), goNext(it2)) 285 { 286 if (isGap(it1)) 287 { 288 if (isGap(it2)) 289 op = 'P'; 290 else if (isClipped(it2)) 291 op = '?'; 292 else 293 op = 'I'; 294 } 295 else if (isClipped(it1)) 296 { 297 op = '?'; 298 } 299 else 300 { 301 if (isGap(it2)) 302 op = 'D'; 303 else if (isClipped(it2)) 304 op = 'S'; 305 else 306 op = 'M'; 307 // op = ((TVal1)*it1 == (TVal2)*it2)? '=': 'X'; 308 } 309 310 // append CIGAR operation 311 if (lastOp != op) 312 { 313 if (lastOp == 'D' && numOps >= (unsigned)splicedGapThresh) 314 lastOp = 'N'; 315 if (numOps > 0) 316 { 317 std::stringstream num; 318 num << numOps; 319 append(cigar, num.str()); 320 appendValue(cigar, lastOp); 321 } 322 numOps = 0; 323 lastOp = op; 324 } 325 ++numOps; 326 } 327 // if (atEnd(it1) != atEnd(it2)) 328 // std::cerr << "Invalid pairwise alignment:" << std::endl << gaps1 << std::endl << gaps2 << std::endl; 329 SEQAN_CHECK(atEnd(it1) == atEnd(it2), "Cannot get CIGAR from invalid pairwise alignment!"); 330 if (lastOp == 'D' && numOps >= (unsigned)splicedGapThresh) 331 lastOp = 'N'; 332 if (numOps > 0) 333 { 334 std::stringstream num; 335 num << numOps; 336 append(cigar, num.str()); 337 appendValue(cigar, lastOp); 338 } 339 } 340 341 template < 342 typename TCigar, 343 typename TGaps1, 344 typename TGaps2> 345 inline void 346 getCigarString( 347 TCigar &cigar, 348 TGaps1 &gaps1, // typically reference 349 TGaps2 &gaps2) // typically read 350 { 351 return getCigarString(cigar, gaps1, gaps2, 20); 352 } 353 354 template < 355 typename TOperation, 356 typename TCount, 357 typename TSpec, 358 typename TGaps1, 359 typename TGaps2, 360 typename TThresh> 361 inline void 362 getCigarString( 363 String<CigarElement<TOperation, TCount>, TSpec> &cigar, 364 TGaps1 &gaps1, 365 TGaps2 &gaps2, 366 TThresh splicedGapThresh) 367 { 368 typename Iterator<TGaps1>::Type it1 = begin(gaps1); 369 typename Iterator<TGaps2>::Type it2 = begin(gaps2); 370 // typedef typename Value<typename Host<TGaps1>::Type>::Type TVal1; 371 // typedef typename Value<typename Host<TGaps2>::Type>::Type TVal2; 372 373 clear(cigar); 374 char op = '?', lastOp = ' '; 375 unsigned numOps = 0; 376 377 // std::cout << gaps1 << std::endl; 378 // std::cout << gaps2 << std::endl; 379 for (; !atEnd(it1) && !atEnd(it2); goNext(it1), goNext(it2)) 380 { 381 if (isGap(it1)) 382 { 383 if (isGap(it2)) 384 op = 'P'; 385 else if (isClipped(it2)) 386 op = '?'; 387 else 388 op = 'I'; 389 } 390 else if (isClipped(it1)) 391 { 392 op = '?'; 393 } 394 else 395 { 396 if (isGap(it2)) 397 op = 'D'; 398 else if (isClipped(it2)) 399 op = 'S'; 400 else 401 // op = ((TVal1)*it1 == (TVal2)*it2)? '=': 'X'; 402 op = 'M'; 403 } 404 if (lastOp != op) 405 { 406 if (lastOp == 'D' && numOps >= (unsigned)splicedGapThresh) 407 lastOp = 'N'; 408 if (numOps > 0) 409 appendValue(cigar, CigarElement<>(lastOp, numOps)); 410 numOps = 0; 411 lastOp = op; 412 } 413 ++numOps; 414 } 415 SEQAN_ASSERT_EQ(atEnd(it1), atEnd(it2)); 416 if (lastOp == 'D' && numOps >= (unsigned)splicedGapThresh) 417 lastOp = 'N'; 418 if (numOps > 0) 419 appendValue(cigar, CigarElement<>(op, numOps)); 420 } 421 422 // ---------------------------------------------------------------------------- 423 // alignAndGetCigarString() 424 // ---------------------------------------------------------------------------- 425 426 template < 427 typename TCigar, typename TMDString, typename TContig, typename TReadSeq, 428 typename TAlignedRead, typename TErrors > 429 inline void 430 alignAndGetCigarString( 431 TCigar &cigar, TMDString &md, TContig &, TReadSeq &, 432 TAlignedRead &, TErrors &, Nothing const &) 433 { 434 cigar = "*"; 435 clear(md); 436 } 437 438 struct BamAlignFunctorEditDistance 439 { 440 typedef String<GapAnchor<int> > TGapAnchors; 441 442 TGapAnchors contigAnchors, readAnchors; 443 444 template <typename TGaps1, typename TGaps2, typename TErrors> 445 inline int 446 align(TGaps1 &gaps1, TGaps2 &gaps2, TErrors maxErrors) 447 { 448 return -globalAlignment( 449 gaps1, gaps2, 450 Score<short, EditDistance>(), 451 -(int)maxErrors, (int)maxErrors 452 ); 453 } 454 }; 455 456 struct BamAlignFunctorSemiGlobalGotoh 457 { 458 typedef String<GapAnchor<int> > TGapAnchors; 459 460 Score<int> score; 461 TGapAnchors contigAnchors, readAnchors; 462 463 BamAlignFunctorSemiGlobalGotoh(Score<int> score_) : 464 score(score_) 465 {} 466 467 template <typename TGaps1, typename TGaps2, typename TErrors> 468 inline int 469 align(TGaps1 &gaps1, TGaps2 &gaps2, TErrors maxErrors) 470 { 471 return globalAlignment( 472 gaps1, gaps2, score, 473 AlignConfig<true, false, false, true>(), 474 -(int)maxErrors, (int)maxErrors, 475 Gotoh() 476 ) / scoreMismatch(score); 477 } 478 }; 479 480 struct BamAlignFunctorDefault 481 { 482 }; 483 484 template < 485 typename TCigar, typename TMDString, typename TContigInfix, typename TReadSeq, 486 typename TAlignedRead, typename TErrors, typename TAlignFunctor> 487 inline void 488 _alignAndGetCigarString( 489 TCigar &cigar, TMDString &md, TContigInfix const &contigInfix, TReadSeq const &fwdReadSeq, 490 TAlignedRead &, TErrors &errors, TAlignFunctor &functor) 491 { 492 typedef Gaps<TContigInfix, AnchorGaps<typename TAlignFunctor::TGapAnchors> > TContigGaps; 493 typedef Gaps<TReadSeq, AnchorGaps<typename TAlignFunctor::TGapAnchors> > TReadGaps; 494 495 clear(functor.contigAnchors); 496 clear(functor.readAnchors); 497 498 TContigGaps contigGaps(contigInfix, functor.contigAnchors); 499 TReadGaps readGaps(fwdReadSeq, functor.readAnchors); 500 501 // if there is already an alignment between contigInfix and fwdReadSeq with 0 or 1 error then 502 // we don't to realign as it contains no gaps 503 if (!(errors == 0 || (errors == 1 && length(contigInfix) == length(fwdReadSeq)))) 504 errors = functor.align(contigGaps, readGaps, errors); 505 506 getCigarString(cigar, contigGaps, readGaps); 507 TErrors mdErrors = getMDString(md, contigGaps, readGaps); 508 509 ignoreUnusedVariableWarning(mdErrors); 510 SEQAN_ASSERT_EQ(errors, mdErrors); 511 } 512 513 template < 514 typename TCigar, typename TMDString, typename TContig, typename TReadSeq, 515 typename TAlignedRead, typename TErrors, typename TAlignFunctor> 516 inline void 517 alignAndGetCigarString( 518 TCigar &cigar, TMDString &md, TContig const &contig, TReadSeq const &readSeq, 519 TAlignedRead &alignedRead, TErrors &errors, TAlignFunctor &functor) 520 { 521 typedef typename TContig::TContigSeq TContigSeq; 522 typedef typename Infix<TContigSeq const>::Type TContigInfix; 523 524 TContigInfix contigInfix; 525 526 if (alignedRead.beginPos <= alignedRead.endPos) 527 { 528 contigInfix = infix(contig.seq, alignedRead.beginPos, alignedRead.endPos); 529 _alignAndGetCigarString(cigar, md, contigInfix, readSeq, alignedRead, errors, functor); 530 } 531 else 532 { 533 contigInfix = infix(contig.seq, alignedRead.endPos, alignedRead.beginPos); 534 _alignAndGetCigarString(cigar, md, contigInfix, reverseComplementString(readSeq), alignedRead, errors, functor); 535 } 536 } 537 538 template < 539 typename TCigar, typename TMDString, typename TContig, typename TReadSeq, 540 typename TAlignedRead, typename TErrors> 541 inline void 542 alignAndGetCigarString( 543 TCigar &cigar, TMDString &md, TContig const &contig, TReadSeq const &readSeq, 544 TAlignedRead &alignedRead, TErrors &errors, BamAlignFunctorDefault &) 545 { 546 typedef typename TContig::TContigSeq TContigSeq; 547 typedef Gaps<TContigSeq, AnchorGaps<typename TContig::TGapAnchors> > TContigGaps; 548 typedef typename ReverseComplementString<TReadSeq const>::Type TRefCompReadSeq; 549 typedef Gaps<TReadSeq, AnchorGaps<typename TAlignedRead::TGapAnchors> > TReadGaps; 550 typedef Gaps<TRefCompReadSeq, AnchorGaps<typename TAlignedRead::TGapAnchors> > TRCReadGaps; 551 552 TContigGaps contigGaps(contig.seq, contig.gaps); 553 554 if (alignedRead.beginPos <= alignedRead.endPos) 555 { 556 setClippedBeginPosition(contigGaps, alignedRead.beginPos); 557 setClippedEndPosition(contigGaps, alignedRead.endPos); 558 559 TReadGaps readGaps(readSeq, alignedRead.gaps); 560 561 getCigarString(cigar, contigGaps, readGaps); 562 errors = getMDString(md, contigGaps, readGaps); 563 } 564 else 565 { 566 setClippedBeginPosition(contigGaps, alignedRead.endPos); 567 setClippedEndPosition(contigGaps, alignedRead.beginPos); 568 569 TRCReadGaps readGaps(reverseComplementString(readSeq), alignedRead.gaps); 570 571 getCigarString(cigar, contigGaps, readGaps); 572 errors = getMDString(md, contigGaps, readGaps); 573 } 574 } 575 576 // ---------------------------------------------------------------------------- 577 // _getClippedLength() 578 // ---------------------------------------------------------------------------- 579 580 template <typename TCigarString, typename TNum> 581 inline void _getClippedLength(TNum & sum, TCigarString const & cigar) 582 { 583 typedef typename Iterator<TCigarString const, Standard>::Type TCigarIter; 584 585 TCigarIter it = begin(cigar, Standard()); 586 TCigarIter itEnd = end(cigar, Standard()); 587 588 sum = 0; 589 for (; it != itEnd; ++it) 590 if (getValue(it).operation != 'S' && getValue(it).operation != 'H') 591 sum += getValue(it).count; 592 } 593 594 // ---------------------------------------------------------------------------- 595 // _getLengthInRef() 596 // ---------------------------------------------------------------------------- 597 598 template <typename TCigarString, typename TNum> 599 inline void _getLengthInRef(TNum & sum, TCigarString const & cigar) 600 { 601 typedef typename Iterator<TCigarString const, Standard>::Type TCigarIter; 602 603 TCigarIter it = begin(cigar, Standard()); 604 TCigarIter itEnd = end(cigar, Standard()); 605 606 sum = 0; 607 for (; it != itEnd; ++it) 608 if (getValue(it).operation != 'S' && getValue(it).operation != 'H' && getValue(it).operation != 'I') 609 sum += getValue(it).count; 610 } 611 612 // ---------------------------------------------------------------------------- 613 // _getQueryLength() 614 // ---------------------------------------------------------------------------- 615 616 template <typename TCigarString> 617 inline typename Size<typename Value<TCigarString>::Type>::Type 618 _getQueryLength(TCigarString const & cigar) 619 { 620 typedef typename Iterator<TCigarString const, Standard>::Type TCigarIter; 621 typedef typename Size<typename Value<TCigarString>::Type>::Type TSize; 622 TCigarIter it = begin(cigar, Standard()); 623 TCigarIter itEnd = end(cigar, Standard()); 624 625 TSize len = 0; 626 for (; it != itEnd; ++it) 627 if (getValue(it).operation != 'D' && getValue(it).operation != 'H' && getValue(it).operation != 'N' && getValue(it).operation != 'P') 628 len += getValue(it).count; 629 return len; 630 } 631 632 // ---------------------------------------------------------------------------- 633 // cigarToGapAnchorRead() 634 // ---------------------------------------------------------------------------- 635 636 template <typename TGaps, typename TCigarString> 637 unsigned cigarToGapAnchorRead(TGaps & gaps, TCigarString const & cigar) 638 { 639 typename Iterator<TGaps>::Type it = begin(gaps); 640 bool atBegin = true; 641 unsigned beginGaps = 0; 642 for (unsigned i = 0; i < length(cigar); ++i) 643 { 644 switch (cigar[i].operation) 645 { 646 case 'D': 647 case 'N': 648 case 'P': 649 if (atBegin) 650 beginGaps += cigar[i].count; 651 insertGaps(it, cigar[i].count); 652 SEQAN_FALLTHROUGH 653 case 'I': 654 case 'M': 655 case 'S': 656 it += cigar[i].count; 657 atBegin = false; 658 } 659 } 660 return beginGaps; 661 } 662 663 // ---------------------------------------------------------------------------- 664 // cigarToGapAnchorContig() 665 // ---------------------------------------------------------------------------- 666 667 template<typename TCigarString, typename TGaps> 668 unsigned cigarToGapAnchorContig(TGaps & gaps, TCigarString const & cigar) 669 { 670 typename Iterator<TGaps>::Type it = begin(gaps); 671 bool atBegin = true; 672 unsigned beginGaps = 0; 673 for (unsigned i = 0; i < length(cigar); ++i) 674 { 675 switch (cigar[i].operation) 676 { 677 case 'I': 678 case 'P': 679 if (atBegin) 680 beginGaps += cigar[i].count; 681 insertGaps(it, cigar[i].count); 682 SEQAN_FALLTHROUGH 683 case 'D': 684 case 'M': 685 case 'N': 686 case 'S': 687 it += cigar[i].count; 688 atBegin = false; 689 } 690 } 691 return beginGaps; 692 } 693 694 } // namespace seqan 695 696 #endif // #ifndef INCLUDE_SEQAN_BAM_IO_CIGAR_H_ 697