yat/omic/Pileup.h

Code
Comments
Other
Rev Date Author Line
3310 25 Aug 14 peter 1 #ifndef theplu_yat_omic_pileup
3310 25 Aug 14 peter 2 #define theplu_yat_omic_pileup
3310 25 Aug 14 peter 3
3310 25 Aug 14 peter 4 // $Id$
3310 25 Aug 14 peter 5
3310 25 Aug 14 peter 6 /*
3999 08 Oct 20 peter 7   Copyright (C) 2014, 2016, 2017, 2018, 2019, 2020 Peter Johansson
3310 25 Aug 14 peter 8
3310 25 Aug 14 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
3310 25 Aug 14 peter 10
3310 25 Aug 14 peter 11   The yat library is free software; you can redistribute it and/or
3310 25 Aug 14 peter 12   modify it under the terms of the GNU General Public License as
3310 25 Aug 14 peter 13   published by the Free Software Foundation; either version 3 of the
3310 25 Aug 14 peter 14   License, or (at your option) any later version.
3310 25 Aug 14 peter 15
3310 25 Aug 14 peter 16   The yat library is distributed in the hope that it will be useful,
3310 25 Aug 14 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
3310 25 Aug 14 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3310 25 Aug 14 peter 19   General Public License for more details.
3310 25 Aug 14 peter 20
3310 25 Aug 14 peter 21   You should have received a copy of the GNU General Public License
3758 17 Oct 18 peter 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
3310 25 Aug 14 peter 23 */
3310 25 Aug 14 peter 24
3310 25 Aug 14 peter 25 #include "BamRead.h"
3310 25 Aug 14 peter 26
3336 24 Oct 14 peter 27 #include "yat/utility/CigarIterator.h"
3327 10 Oct 14 peter 28 #include "yat/utility/yat_assert.h"
3327 10 Oct 14 peter 29
3335 23 Oct 14 peter 30 #include <boost/iterator/iterator_traits.hpp>
3517 16 Aug 16 peter 31 #include <boost/iterator/iterator_concepts.hpp>
3335 23 Oct 14 peter 32 #include <boost/concept/assert.hpp>
3327 10 Oct 14 peter 33
3310 25 Aug 14 peter 34 #include <algorithm>
3332 23 Oct 14 peter 35 #include <list>
4019 06 Nov 20 peter 36 #include <memory>
3310 25 Aug 14 peter 37
3310 25 Aug 14 peter 38 namespace theplu {
3310 25 Aug 14 peter 39 namespace yat {
3310 25 Aug 14 peter 40 namespace omic {
3310 25 Aug 14 peter 41
3317 19 Sep 14 peter 42   /**
3335 23 Oct 14 peter 43      Pileup takes a range of BamReads and pile them up on the current
3335 23 Oct 14 peter 44      position. User can investigate the pile of reads at current
3335 23 Oct 14 peter 45      (locus) position via begin() and end(), which provide access to
3335 23 Oct 14 peter 46      reads at current position. The position can be incremented and in
3335 23 Oct 14 peter 47      that way the user can investigate each position of the genome.
3335 23 Oct 14 peter 48
3317 19 Sep 14 peter 49      \since new in yat 0.13
3335 23 Oct 14 peter 50
3335 23 Oct 14 peter 51      Type Requirments:
3517 16 Aug 16 peter 52      - \c Iterator is a \readable_iterator
3517 16 Aug 16 peter 53      - \c Iterator is a \single_pass_iterator
3335 23 Oct 14 peter 54      - \c value_type is convertible to BamRead
3317 19 Sep 14 peter 55    */
3310 25 Aug 14 peter 56   template<typename Iterator>
3310 25 Aug 14 peter 57   class Pileup
3310 25 Aug 14 peter 58   {
3310 25 Aug 14 peter 59   public:
3310 25 Aug 14 peter 60     /**
3310 25 Aug 14 peter 61        Essentially a BamRead but with additional information regarding
3310 25 Aug 14 peter 62        base we currently pointing to
3310 25 Aug 14 peter 63      */
3310 25 Aug 14 peter 64     class Entry
3310 25 Aug 14 peter 65     {
3310 25 Aug 14 peter 66     public:
3310 25 Aug 14 peter 67       /**
3332 23 Oct 14 peter 68          \brief default constructor
3332 23 Oct 14 peter 69        */
3332 23 Oct 14 peter 70       Entry(void) {}
3332 23 Oct 14 peter 71
3332 23 Oct 14 peter 72       /**
3334 23 Oct 14 peter 73          Create an Entry with read \a b. Entry is initialized to point
3334 23 Oct 14 peter 74          to first aligned base (b.pos()).
3310 25 Aug 14 peter 75        */
3332 23 Oct 14 peter 76       Entry(const BamRead& b);
3310 25 Aug 14 peter 77
3310 25 Aug 14 peter 78       /**
3328 10 Oct 14 peter 79          \return const reference to BamRead
3310 25 Aug 14 peter 80        */
3332 23 Oct 14 peter 81       const BamRead& bam(void) const
3332 23 Oct 14 peter 82       {
3332 23 Oct 14 peter 83         YAT_ASSERT(bam_.get());
3332 23 Oct 14 peter 84         return *bam_;
3332 23 Oct 14 peter 85       }
3310 25 Aug 14 peter 86
3317 19 Sep 14 peter 87       /**
3332 23 Oct 14 peter 88          \return cigar operation at current position
3317 19 Sep 14 peter 89        */
3332 23 Oct 14 peter 90       uint8_t cigar_op(void) const { return *cigar_; }
3332 23 Oct 14 peter 91
3332 23 Oct 14 peter 92
3332 23 Oct 14 peter 93       /**
3332 23 Oct 14 peter 94          \return length of current cigar element
3332 23 Oct 14 peter 95        */
3332 23 Oct 14 peter 96       uint32_t cigar_oplen(void) const
3310 25 Aug 14 peter 97       {
3332 23 Oct 14 peter 98         YAT_ASSERT(bam_.get());
3332 23 Oct 14 peter 99         YAT_ASSERT(cigar_.base() >= bam_->cigar());
3332 23 Oct 14 peter 100         return bam_cigar_oplen(*cigar_.base());
3310 25 Aug 14 peter 101       }
3310 25 Aug 14 peter 102
3332 23 Oct 14 peter 103
3310 25 Aug 14 peter 104       /**
3310 25 Aug 14 peter 105          \return base quality of current position
3310 25 Aug 14 peter 106        */
3310 25 Aug 14 peter 107       unsigned short qual(void) const
3310 25 Aug 14 peter 108       {
3332 23 Oct 14 peter 109         YAT_ASSERT(bam_.get());
3332 23 Oct 14 peter 110         YAT_ASSERT(qpos_ < bam_->sequence_length());
3332 23 Oct 14 peter 111         return bam_->qual(qpos_);
3310 25 Aug 14 peter 112       }
3310 25 Aug 14 peter 113
3310 25 Aug 14 peter 114       /**
3646 18 May 17 peter 115          Returned base is encoded in 4 bits and can be converted to
3646 18 May 17 peter 116          char using nt16_str(uint8_t).
3646 18 May 17 peter 117
3310 25 Aug 14 peter 118          \return base of current position
3310 25 Aug 14 peter 119        */
3310 25 Aug 14 peter 120       uint8_t sequence(void) const;
3310 25 Aug 14 peter 121
3310 25 Aug 14 peter 122       /**
3333 23 Oct 14 peter 123          Increment CIGAR. If current CIGAR element is consuming query
3333 23 Oct 14 peter 124          (not deletion), increment to next base on query.
3310 25 Aug 14 peter 125        */
3310 25 Aug 14 peter 126       void increment(void);
3310 25 Aug 14 peter 127
3310 25 Aug 14 peter 128       /**
3310 25 Aug 14 peter 129          \return true if this position was deleted in this query
3310 25 Aug 14 peter 130        */
3310 25 Aug 14 peter 131       bool is_del(void) const
3310 25 Aug 14 peter 132       {
3332 23 Oct 14 peter 133         return cigar_type()==2;
3310 25 Aug 14 peter 134       }
3310 25 Aug 14 peter 135
3317 19 Sep 14 peter 136       /**
3317 19 Sep 14 peter 137          \return index of base at this position
3317 19 Sep 14 peter 138        */
3310 25 Aug 14 peter 139       size_t qpos(void) const { return qpos_; }
3310 25 Aug 14 peter 140     private:
4019 06 Nov 20 peter 141       std::shared_ptr<BamRead> bam_;
3310 25 Aug 14 peter 142       // index of base pointed to
3310 25 Aug 14 peter 143       size_t qpos_;
3336 24 Oct 14 peter 144       yat::utility::CigarIterator<const uint32_t*> cigar_;
3310 25 Aug 14 peter 145       size_t skip_;
3332 23 Oct 14 peter 146       // return 0, 1, 2, or 3 see bam_cigar_type in yat/utility/Cigar.h
3332 23 Oct 14 peter 147       uint8_t cigar_type(void) const { return bam_cigar_type(cigar_op()); }
3310 25 Aug 14 peter 148       friend class Pileup;
3310 25 Aug 14 peter 149     };
3310 25 Aug 14 peter 150
3327 10 Oct 14 peter 151   public:
3317 19 Sep 14 peter 152     /**
3327 10 Oct 14 peter 153        Const iterator that can be used to iterate over reads that
3647 22 May 17 peter 154        overlap with current position (as defined by tid() and pos()).
3317 19 Sep 14 peter 155      */
3332 23 Oct 14 peter 156     typedef typename std::list<Entry>::const_iterator const_iterator;
3310 25 Aug 14 peter 157
3310 25 Aug 14 peter 158     /**
3801 04 May 19 peter 159        Same as const_iterator but iterates reversed.
3801 04 May 19 peter 160
3801 04 May 19 peter 161        \see rbegin(void) rend(void)
3801 04 May 19 peter 162
3801 04 May 19 peter 163        \since new in yat 0.17
3801 04 May 19 peter 164      */
3801 04 May 19 peter 165     typedef typename std::list<Entry>::const_reverse_iterator
3801 04 May 19 peter 166     const_reverse_iterator;
3801 04 May 19 peter 167
3801 04 May 19 peter 168     /**
3310 25 Aug 14 peter 169        Create an empty Pileup
3310 25 Aug 14 peter 170      */
3310 25 Aug 14 peter 171     Pileup(void) {}
3310 25 Aug 14 peter 172
3310 25 Aug 14 peter 173     /**
3334 23 Oct 14 peter 174        Create a Pileup that will use reads in range [first, last). The
3334 23 Oct 14 peter 175        range must be sorted as defined by BamLessPos.
3334 23 Oct 14 peter 176
3334 23 Oct 14 peter 177        Iterator is a \readable_iterator.
3310 25 Aug 14 peter 178      */
3310 25 Aug 14 peter 179     Pileup(Iterator first, Iterator last);
3310 25 Aug 14 peter 180
3310 25 Aug 14 peter 181     /**
3327 10 Oct 14 peter 182        \return an iterator that points to the first Entry in Pileup.
3327 10 Oct 14 peter 183
3327 10 Oct 14 peter 184        Iteration is done in same order as range provided in
3327 10 Oct 14 peter 185        constructor. Reads that do not overlap with current position
3327 10 Oct 14 peter 186        (tid:pos) are ignored.
3310 25 Aug 14 peter 187      */
3327 10 Oct 14 peter 188     const_iterator begin(void) const;
3310 25 Aug 14 peter 189
3310 25 Aug 14 peter 190     /**
3327 10 Oct 14 peter 191        \return an iterator that points one past the last element in
3327 10 Oct 14 peter 192        the Pileup.
3310 25 Aug 14 peter 193      */
3310 25 Aug 14 peter 194     const_iterator end(void) const;
3310 25 Aug 14 peter 195
3310 25 Aug 14 peter 196     /**
3310 25 Aug 14 peter 197        \return true if there are more positions to increment into
3310 25 Aug 14 peter 198      */
3310 25 Aug 14 peter 199     bool good(void) const { return !(bam_ == bam_end_ && data_.empty()); }
3310 25 Aug 14 peter 200
3310 25 Aug 14 peter 201     /**
3334 23 Oct 14 peter 202        \brief step to next position.
3334 23 Oct 14 peter 203
3334 23 Oct 14 peter 204        If any Entry has a insertion at current position, increment
3334 23 Oct 14 peter 205        each Entry that has an insertion. Otherwise increment to next
3334 23 Oct 14 peter 206        position and update each Entry accordingly.
3335 23 Oct 14 peter 207
3335 23 Oct 14 peter 208        If the position is not covered by any read, the Pileup
3335 23 Oct 14 peter 209        fastforward to next covered locus.
3310 25 Aug 14 peter 210      */
3310 25 Aug 14 peter 211     void increment(void);
3310 25 Aug 14 peter 212
3310 25 Aug 14 peter 213     /**
3802 12 May 19 peter 214        \return number of reads in pileup. Reads with a deletion of the
3802 12 May 19 peter 215        current position is also counted.
3802 12 May 19 peter 216
3802 12 May 19 peter 217        \since new in yat 0.17
3802 12 May 19 peter 218      */
3802 12 May 19 peter 219     size_t size(void) const;
3802 12 May 19 peter 220
3802 12 May 19 peter 221     /**
3801 04 May 19 peter 222        \return an iterator that points to the last Entry in Pileup.
3801 04 May 19 peter 223
3801 04 May 19 peter 224        Iteration is done in reversed order; othwerwise behaves as begin(void)
3801 04 May 19 peter 225
3801 04 May 19 peter 226        \since new in yat 0.17
3801 04 May 19 peter 227      */
3801 04 May 19 peter 228     const_reverse_iterator rbegin(void) const;
3801 04 May 19 peter 229
3801 04 May 19 peter 230     /**
3801 04 May 19 peter 231        \return const reverse iterator to reverse end
3801 04 May 19 peter 232
3801 04 May 19 peter 233        \see end(void)
3801 04 May 19 peter 234        \see rbegin(void)
3801 04 May 19 peter 235
3801 04 May 19 peter 236        \since new in yat 0.17
3801 04 May 19 peter 237      */
3801 04 May 19 peter 238     const_reverse_iterator rend(void) const;
3801 04 May 19 peter 239
3801 04 May 19 peter 240     /**
3328 10 Oct 14 peter 241        \return template id as defined in bam file header
3310 25 Aug 14 peter 242      */
3310 25 Aug 14 peter 243     int32_t tid(void) const { return tid_; }
3310 25 Aug 14 peter 244
3310 25 Aug 14 peter 245     /**
3310 25 Aug 14 peter 246        \return position
3970 13 Aug 20 peter 247
3970 13 Aug 20 peter 248        When in a skip_ref region, the position of the next reference
3970 13 Aug 20 peter 249        is returned.
3310 25 Aug 14 peter 250      */
3310 25 Aug 14 peter 251     int32_t pos(void) const { return pos_; }
3310 25 Aug 14 peter 252
3310 25 Aug 14 peter 253     /**
3310 25 Aug 14 peter 254        \return true if position does not exist in reference, i.e.,
3310 25 Aug 14 peter 255        inserted region.
3310 25 Aug 14 peter 256      */
3310 25 Aug 14 peter 257     bool skip_ref(void) const;
3310 25 Aug 14 peter 258   private:
3332 23 Oct 14 peter 259     typedef typename std::list<Entry>::iterator iterator;
3310 25 Aug 14 peter 260
3310 25 Aug 14 peter 261     Iterator bam_;
3310 25 Aug 14 peter 262     Iterator bam_end_;
3310 25 Aug 14 peter 263     int32_t pos_;
3310 25 Aug 14 peter 264     int32_t tid_;
3332 23 Oct 14 peter 265     std::list<Entry> data_;
3310 25 Aug 14 peter 266     uint32_t skip_ref_;
3310 25 Aug 14 peter 267   };
3310 25 Aug 14 peter 268
3310 25 Aug 14 peter 269
3317 19 Sep 14 peter 270   /**
3317 19 Sep 14 peter 271      \relates Pileup
3317 19 Sep 14 peter 272
3317 19 Sep 14 peter 273      \since new in yat 0.13
3317 19 Sep 14 peter 274    */
3310 25 Aug 14 peter 275   template<typename Iterator>
3310 25 Aug 14 peter 276   Pileup<Iterator> make_pileup(Iterator first, Iterator last)
3310 25 Aug 14 peter 277   {
3310 25 Aug 14 peter 278     return Pileup<Iterator>(first, last);
3310 25 Aug 14 peter 279   }
3310 25 Aug 14 peter 280
3310 25 Aug 14 peter 281   // implementations
3310 25 Aug 14 peter 282
3310 25 Aug 14 peter 283   template<typename Iterator>
3310 25 Aug 14 peter 284   Pileup<Iterator>::Pileup(Iterator first, Iterator last)
3328 10 Oct 14 peter 285     : bam_(first), bam_end_(last), pos_(0), tid_(0), skip_ref_(0)
3310 25 Aug 14 peter 286   {
3517 16 Aug 16 peter 287     BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<Iterator>));
3517 16 Aug 16 peter 288     BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<Iterator>));
3335 23 Oct 14 peter 289     using boost::Convertible;
3335 23 Oct 14 peter 290     typedef typename boost::iterator_value<Iterator>::type value_type;
3335 23 Oct 14 peter 291     BOOST_CONCEPT_ASSERT((Convertible<value_type, BamRead>));
3310 25 Aug 14 peter 292     increment();
3310 25 Aug 14 peter 293   }
3310 25 Aug 14 peter 294
3310 25 Aug 14 peter 295
3310 25 Aug 14 peter 296   template<typename Iterator>
3327 10 Oct 14 peter 297   typename Pileup<Iterator>::const_iterator Pileup<Iterator>::begin(void) const
3327 10 Oct 14 peter 298   {
3332 23 Oct 14 peter 299     return data_.begin();
3327 10 Oct 14 peter 300   }
3327 10 Oct 14 peter 301
3327 10 Oct 14 peter 302
3327 10 Oct 14 peter 303   template<typename Iterator>
3801 04 May 19 peter 304   typename Pileup<Iterator>::const_reverse_iterator
3801 04 May 19 peter 305   Pileup<Iterator>::rbegin(void) const
3801 04 May 19 peter 306   {
3801 04 May 19 peter 307     return data_.rbegin();
3801 04 May 19 peter 308   }
3801 04 May 19 peter 309
3801 04 May 19 peter 310
3801 04 May 19 peter 311   template<typename Iterator>
3310 25 Aug 14 peter 312   typename Pileup<Iterator>::const_iterator Pileup<Iterator>::end(void) const
3310 25 Aug 14 peter 313   {
3332 23 Oct 14 peter 314     return data_.end();
3310 25 Aug 14 peter 315   }
3310 25 Aug 14 peter 316
3310 25 Aug 14 peter 317
3310 25 Aug 14 peter 318   template<typename Iterator>
3801 04 May 19 peter 319   typename Pileup<Iterator>::const_reverse_iterator
3801 04 May 19 peter 320   Pileup<Iterator>::rend(void) const
3801 04 May 19 peter 321   {
3801 04 May 19 peter 322     return data_.rend();
3801 04 May 19 peter 323   }
3801 04 May 19 peter 324
3801 04 May 19 peter 325
3801 04 May 19 peter 326   template<typename Iterator>
3310 25 Aug 14 peter 327   void Pileup<Iterator>::increment(void)
3310 25 Aug 14 peter 328   {
3334 23 Oct 14 peter 329     if (!data_.empty()) {
3310 25 Aug 14 peter 330       // in insertion
3310 25 Aug 14 peter 331       if (skip_ref()) {
3332 23 Oct 14 peter 332         for (iterator it = data_.begin(); it!=data_.end(); ++it) {
3332 23 Oct 14 peter 333           if (it->cigar_type() == 1) {
3332 23 Oct 14 peter 334             it->increment();
3310 25 Aug 14 peter 335           }
3332 23 Oct 14 peter 336           --it->skip_;
3310 25 Aug 14 peter 337         }
3310 25 Aug 14 peter 338         --skip_ref_;
3310 25 Aug 14 peter 339       }
3310 25 Aug 14 peter 340       // walk one base forward on reference
3310 25 Aug 14 peter 341       else {
3310 25 Aug 14 peter 342         ++pos_;
3332 23 Oct 14 peter 343         typename std::list<Entry>::iterator read = data_.begin();
3332 23 Oct 14 peter 344         while (read != data_.end()) {
3332 23 Oct 14 peter 345           // remove entry if read is entirely left of pos
3916 31 May 20 peter 346           if (read->bam().end() <= pos_) {
3332 23 Oct 14 peter 347             data_.erase(read++);
3310 25 Aug 14 peter 348             continue;
3332 23 Oct 14 peter 349           }
3916 31 May 20 peter 350           read->increment();
3310 25 Aug 14 peter 351           // check if we're stepping into insertion
3332 23 Oct 14 peter 352           if (read->cigar_op() == BAM_CINS) {
3332 23 Oct 14 peter 353             skip_ref_ =  std::max(skip_ref_, read->cigar_oplen());
3332 23 Oct 14 peter 354           }
3332 23 Oct 14 peter 355           ++read;
3310 25 Aug 14 peter 356         }
3332 23 Oct 14 peter 357
3310 25 Aug 14 peter 358         // stepped into an insertion
3332 23 Oct 14 peter 359         if (skip_ref_)
3332 23 Oct 14 peter 360           for (iterator iter=data_.begin(); iter!=data_.end(); ++iter)
3332 23 Oct 14 peter 361             iter->skip_ = skip_ref_;
3310 25 Aug 14 peter 362       }
3310 25 Aug 14 peter 363     }
3310 25 Aug 14 peter 364
3334 23 Oct 14 peter 365
3334 23 Oct 14 peter 366     // fast forward to next covered locus
3334 23 Oct 14 peter 367     if (data_.empty()) {
3334 23 Oct 14 peter 368       if (bam_==bam_end_)
3334 23 Oct 14 peter 369         return;
3334 23 Oct 14 peter 370       tid_ = bam_->core().tid;
3334 23 Oct 14 peter 371       pos_ = bam_->core().pos;
3334 23 Oct 14 peter 372     }
3334 23 Oct 14 peter 373
3310 25 Aug 14 peter 374     // read data
3310 25 Aug 14 peter 375     while (bam_!=bam_end_ && bam_->core().tid==tid_ && bam_->pos()<=pos_) {
3310 25 Aug 14 peter 376       // skip unmapped reads
3310 25 Aug 14 peter 377       if (!(bam_->core().flag & BAM_FUNMAP)) {
3310 25 Aug 14 peter 378         data_.push_back(Entry(*bam_));
3310 25 Aug 14 peter 379       }
3310 25 Aug 14 peter 380       ++bam_;
3310 25 Aug 14 peter 381     }
3310 25 Aug 14 peter 382   }
3310 25 Aug 14 peter 383
3310 25 Aug 14 peter 384
3310 25 Aug 14 peter 385   template<typename Iterator>
3310 25 Aug 14 peter 386   bool Pileup<Iterator>::skip_ref(void) const
3310 25 Aug 14 peter 387   {
3310 25 Aug 14 peter 388     return skip_ref_;
3310 25 Aug 14 peter 389   }
3310 25 Aug 14 peter 390
3310 25 Aug 14 peter 391
3310 25 Aug 14 peter 392   template<typename Iterator>
3802 12 May 19 peter 393   size_t Pileup<Iterator>::size(void) const
3802 12 May 19 peter 394   {
3802 12 May 19 peter 395     return data_.size();
3802 12 May 19 peter 396   }
3802 12 May 19 peter 397
3802 12 May 19 peter 398
3802 12 May 19 peter 399   template<typename Iterator>
3332 23 Oct 14 peter 400   Pileup<Iterator>::Entry::Entry(const BamRead& b)
3332 23 Oct 14 peter 401     : bam_(new BamRead(b)), qpos_(0), cigar_(bam_->cigar()), skip_(0)
3310 25 Aug 14 peter 402   {
3332 23 Oct 14 peter 403     YAT_ASSERT(cigar_.base() == bam_->cigar());
3332 23 Oct 14 peter 404     if (bam_->core().n_cigar==0)
3310 25 Aug 14 peter 405       return;
3310 25 Aug 14 peter 406     // iterate cigar to first match
3332 23 Oct 14 peter 407     while (qpos_ < bam_->sequence_length() && cigar_type() != 3) {
3332 23 Oct 14 peter 408       // type is either 0 (hardclipped) or 1 (softclipped)
3332 23 Oct 14 peter 409       YAT_ASSERT(cigar_type()!=2);
3332 23 Oct 14 peter 410       if (cigar_type() == 1)
3332 23 Oct 14 peter 411         ++qpos_;
3332 23 Oct 14 peter 412       ++cigar_;
3310 25 Aug 14 peter 413     }
3310 25 Aug 14 peter 414   }
3310 25 Aug 14 peter 415
3310 25 Aug 14 peter 416
3310 25 Aug 14 peter 417   template<typename Iterator>
3310 25 Aug 14 peter 418   void Pileup<Iterator>::Entry::increment(void)
3310 25 Aug 14 peter 419   {
3332 23 Oct 14 peter 420     YAT_ASSERT(cigar_.base() >= bam_->cigar());
3332 23 Oct 14 peter 421     if (cigar_type() & 1)
3310 25 Aug 14 peter 422       ++qpos_;
3332 23 Oct 14 peter 423     ++cigar_;
3332 23 Oct 14 peter 424     YAT_ASSERT(cigar_.base() >= bam_->cigar());
3310 25 Aug 14 peter 425   }
3310 25 Aug 14 peter 426
3310 25 Aug 14 peter 427
3310 25 Aug 14 peter 428   template<typename Iterator>
3310 25 Aug 14 peter 429   uint8_t Pileup<Iterator>::Entry::sequence(void) const
3310 25 Aug 14 peter 430   {
3310 25 Aug 14 peter 431     if (skip_==0) {
3332 23 Oct 14 peter 432       if (cigar_type() & 1)
3332 23 Oct 14 peter 433         return bam_->sequence(qpos_);
3332 23 Oct 14 peter 434       return 0;
3310 25 Aug 14 peter 435     }
3310 25 Aug 14 peter 436     if (cigar_op()!=BAM_CINS)
3310 25 Aug 14 peter 437       return 0;
3332 23 Oct 14 peter 438     return bam_->sequence(qpos_);
3310 25 Aug 14 peter 439   }
3310 25 Aug 14 peter 440
3310 25 Aug 14 peter 441 }}}
3310 25 Aug 14 peter 442 #endif

Generated on Sat Oct 14 02:22:53 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)