Ignore:
Timestamp:
Oct 10, 2014, 10:30:20 AM (7 years ago)
Author:
Peter
Message:

redefine Pileup::const_iterator so it avoids reads that do not overlap current position. refs #806

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/yat/omic/Pileup.h

    r3317 r3327  
    2525#include "BamRead.h"
    2626
     27#include "yat/utility/yat_assert.h"
     28
     29#include <boost/iterator/filter_iterator.hpp>
     30
    2731#include <algorithm>
    28 #include <cassert>
    2932#include <deque>
    3033
     
    4043  {
    4144  public:
    42 
    4345    /**
    4446       Essentially a BamRead but with additional information regarding
     
    6365      uint32_t cigar(void) const
    6466      {
    65         assert(c_index_ < bam_.core().n_cigar);
     67        YAT_ASSERT(c_index_ < bam_.core().n_cigar);
    6668        return bam_.cigar(c_index_);
    6769      }
     
    7274      unsigned short qual(void) const
    7375      {
    74         assert(qpos_ < bam_.sequence_length());
     76        YAT_ASSERT(qpos_ < bam_.sequence_length());
    7577        return bam_.qual(qpos_);
    7678      }
     
    8890         reference that has coverage.
    8991       */
     92      // FIXME should this be public?
    9093      void increment(void);
    9194
     
    104107      bool is_del(void) const
    105108      {
    106         assert(c_index_ < bam_.core().n_cigar);
     109        YAT_ASSERT(c_index_ < bam_.core().n_cigar);
    107110        return bam_cigar_type(cigar())==2;
    108111      }
     
    127130    };
    128131
    129     /**
    130        \brief Const iterator that can be used to iterate over reads at
    131        current position
    132      */
    133     typedef typename std::deque<Entry>::const_iterator const_iterator;
     132  private:
     133    class Filter
     134    {
     135    public:
     136      Filter(void) {}
     137      Filter(int32_t tid, int32_t pos) : tid_(tid), pos_(pos) {}
     138      // return true if Entry read overlaps with tid:pos
     139      bool operator()(const Entry& read) const
     140      {
     141        return read.bam().pos() <= pos_ && pos_ < read.bam().end()
     142          && read.bam().tid() == tid_;
     143      }
     144    private:
     145      int32_t tid_;
     146      int32_t pos_;
     147    };
     148
     149  public:
     150    /**
     151       Const iterator that can be used to iterate over reads that
     152       overlap with current position (as deined by tid() and pos()).
     153     */
     154    typedef boost::filter_iterator<
     155    Filter
     156    , typename std::deque<Entry>::const_iterator> const_iterator;
    134157
    135158    /**
     
    144167
    145168    /**
    146        first Entry in pileup
    147      */
    148     const_iterator begin(void) const { return data_.begin(); }
    149 
    150     /**
    151        last Entry in Pileup
     169       \return an iterator that points to the first Entry in Pileup.
     170
     171       Iteration is done in same order as range provided in
     172       constructor. Reads that do not overlap with current position
     173       (tid:pos) are ignored.
     174     */
     175    const_iterator begin(void) const;
     176
     177    /**
     178       \return an iterator that points one past the last element in
     179       the Pileup.
    152180     */
    153181    const_iterator end(void) const;
     
    214242
    215243  template<typename Iterator>
     244  typename Pileup<Iterator>::const_iterator Pileup<Iterator>::begin(void) const
     245  {
     246    Filter filter(tid(), pos());
     247    typename std::deque<Entry>::const_iterator iter = data_.begin();
     248    typename std::deque<Entry>::const_iterator end = data_.end();
     249    return boost::make_filter_iterator(filter, iter, end);
     250  }
     251
     252
     253  template<typename Iterator>
    216254  typename Pileup<Iterator>::const_iterator Pileup<Iterator>::end(void) const
    217255  {
    218     if (data_.empty())
    219       return data_.end();
    220     // if last entry is passed pos, then don't include it in [begin, end) range
    221     if ((data_.back().bam().tid() > tid() ||
    222          data_.back().bam().pos() > pos())) {
    223       return data_.begin() + (data_.size()-1);
    224     }
    225 
    226     return data_.end();
     256    Filter filter(tid(), pos());
     257    typename std::deque<Entry>::const_iterator end = data_.end();
     258    return boost::make_filter_iterator(filter, end, end);
    227259  }
    228260
     
    302334    // iterate cigar to first match
    303335    while (bam_cigar_type(bam_.cigar(c_index_)) != 3) {
    304       assert(bam_cigar_type(bam_.cigar(c_index_))!=2);
     336      YAT_ASSERT(bam_cigar_type(bam_.cigar(c_index_))!=2);
    305337      if (bam_cigar_type(bam_.cigar(c_index_)) == 1) {
    306338        qpos_ += b.cigar_oplen(c_index_);
     
    311343        return;
    312344      }
    313       assert(c_index_ < bam_.core().n_cigar);
     345      YAT_ASSERT(c_index_ < bam_.core().n_cigar);
    314346    }
    315347  }
     
    344376  {
    345377    ++c_index2_;
    346     assert(c_index_<bam_.core().n_cigar);
     378    YAT_ASSERT(c_index_<bam_.core().n_cigar);
    347379    if (c_index2_ >= bam_.cigar_oplen(c_index_)) {
    348380      c_index2_ = 0;
Note: See TracChangeset for help on using the changeset viewer.