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