Changeset 3332 for trunk/yat/omic/Pileup.h

Timestamp:

Oct 23, 2014, 1:18:15 PM (8 years ago)

Author:

Peter

Message:

In Pileup store reads in a std::list (rather than std::deque) to allow
easy erase as soon as we are one passed last base. Erasing reads
immediately implies we don't need a filter_iterator, but can use
list::const_iterator directly.

Fix a bug in Pileup::Entry, which held a BamRead? and CigarIterator? and
the latter is essentially a uint32_t* that is returned from
BamRead::cigar(). This caused problems in copy/assignment since the
BamRead? was copied by value and the CigarIterator? didn't point to
BamRead? in this class but to BamRead? in rhs. The problem was solved by
replacing the BamRead? by a shared_ptr<BamRead?>.

refs #806

File:

• trunk/yat/omic/Pileup.h (modified) (16 diffs)

trunk/yat/omic/Pileup.h

@@ -24 +24 @@
 
 #include "BamRead.h"
+#include "CigarIterator.h"
 
 #include "yat/utility/yat_assert.h"
 
-#include <boost/iterator/filter_iterator.hpp>
+#include <boost/shared_ptr.hpp>
 
 #include <algorithm>
-#include <deque>
+#include <list>
 
 namespace theplu {
@@ -51 +52 @@
     public:
       /**
+         \brief default constructor
+       */
+      Entry(void) {}
+
+      /**
          Create an Entry with read \a b
        */
-      Entry(const yat::omic::BamRead& b);
+      Entry(const BamRead& b);
 
       /**
          \return const reference to BamRead
        */
-      const yat::omic::BamRead& bam(void) const { return bam_; }
-
-      /**
-         \return cigar element at current position
-       */
-      uint32_t cigar(void) const
+      const BamRead& bam(void) const
       {
-        YAT_ASSERT(c_index_ < bam_.core().n_cigar);
-        return bam_.cigar(c_index_);
-      }
+        YAT_ASSERT(bam_.get());
+        return *bam_;
+      }
+
+      /**
+         \return cigar operation at current position
+       */
+      uint8_t cigar_op(void) const { return *cigar_; }
+
+
+      /**
+         \return length of current cigar element
+       */
+      uint32_t cigar_oplen(void) const
+      {
+        YAT_ASSERT(bam_.get());
+        YAT_ASSERT(cigar_.base() >= bam_->cigar());
+        return bam_cigar_oplen(*cigar_.base());
+      }
+
 
       /**
@@ -74 +92 @@
       unsigned short qual(void) const
       {
-        YAT_ASSERT(qpos_ < bam_.sequence_length());
-        return bam_.qual(qpos_);
+        YAT_ASSERT(bam_.get());
+        YAT_ASSERT(qpos_ < bam_->sequence_length());
+        return bam_->qual(qpos_);
       }
 
@@ -94 +113 @@
 
       /**
-         \return true if beyond last (aligned) base, i.e., we don't
-         count clipped bases
-       */
-      bool end_of_sequence(void) const
-      {
-        return qpos_ >= bam_.sequence_length();
-      }
-
-      /**
          \return true if this position was deleted in this query
        */
       bool is_del(void) const
       {
-        YAT_ASSERT(c_index_ < bam_.core().n_cigar);
-        return bam_cigar_type(cigar())==2;
+        return cigar_type()==2;
       }
 
@@ -116 +125 @@
       size_t qpos(void) const { return qpos_; }
     private:
-      yat::omic::BamRead bam_;
+      boost::shared_ptr<BamRead> bam_;
       // index of base pointed to
       size_t qpos_;
-      size_t c_index_;
-      size_t c_index2_;
+      CigarIterator cigar_;
       size_t skip_;
-      void step_cigar(void);
-      uint32_t cigar_op(void) const
-      {
-        return bam_cigar_op(cigar());
-      }
+      // return 0, 1, 2, or 3 see bam_cigar_type in yat/utility/Cigar.h
+      uint8_t cigar_type(void) const { return bam_cigar_type(cigar_op()); }
       friend class Pileup;
     };
 
-  private:
-
-    /**
-       Private class used in const_iterator.
-     */
-    class Filter
-    {
-    public:
-      Filter(void) {}
-      Filter(int32_t tid, int32_t pos) : tid_(tid), pos_(pos) {}
-      // return true if Entry read overlaps with tid:pos
-      bool operator()(const Entry& read) const
-      {
-        return read.bam().pos() <= pos_ && pos_ < read.bam().end()
-          && read.bam().tid() == tid_;
-      }
-    private:
-      int32_t tid_;
-      int32_t pos_;
-    };
-
   public:
     /**
@@ -156 +140 @@
        overlap with current position (as deined by tid() and pos()).
      */
-    typedef boost::filter_iterator<
-    Filter
-    , typename std::deque<Entry>::const_iterator> const_iterator;
+    typedef typename std::list<Entry>::const_iterator const_iterator;
 
     /**
@@ -211 +193 @@
     bool skip_ref(void) const;
   private:
-    void step_cigar(void);
+    typedef typename std::list<Entry>::iterator iterator;
 
     Iterator bam_;
@@ -217 +199 @@
     int32_t pos_;
     int32_t tid_;
-    std::deque<Entry> data_;
+    std::list<Entry> data_;
     uint32_t skip_ref_;
   };
@@ -246 +228 @@
   typename Pileup<Iterator>::const_iterator Pileup<Iterator>::begin(void) const
   {
-    return boost::make_filter_iterator(Filter(tid(), pos()),
-                                       data_.begin(), data_.end());
+    return data_.begin();
   }
 
@@ -254 +235 @@
   typename Pileup<Iterator>::const_iterator Pileup<Iterator>::end(void) const
   {
-    return boost::make_filter_iterator(Filter(),
-                                       data_.end(), data_.end());
+    return data_.end();
   }
 
@@ -272 +252 @@
       // in insertion
       if (skip_ref()) {
-        for (size_t i=0; i<data_.size(); ++i) {
-          // ignore reads left of current position
-          if (data_[i].end_of_sequence())
-            continue;
-          if (bam_cigar_type(data_[i].cigar())==1) {
-            data_[i].increment();
+        for (iterator it = data_.begin(); it!=data_.end(); ++it) {
+          if (it->cigar_type() == 1) {
+            it->increment();
           }
-          --data_[i].skip_;
+          --it->skip_;
         }
         --skip_ref_;
@@ -286 +263 @@
       else {
         ++pos_;
-        for (size_t i=0; i<data_.size(); ++i) {
-          if (data_[i].end_of_sequence())
+        typename std::list<Entry>::iterator read = data_.begin();
+        while (read != data_.end()) {
+          read->increment();
+          // remove entry if read is entirely left of pos
+          if (read->bam().end() <= pos_) { // FIXME speedup possible
+            data_.erase(read++);
             continue;
-          data_[i].increment();
+          }
           // check if we're stepping into insertion
-          if (!data_[i].end_of_sequence() && (data_[i].cigar_op() == BAM_CINS))
-            skip_ref_ = std::max(skip_ref_, bam_cigar_oplen(data_[i].cigar()));
+          if (read->cigar_op() == BAM_CINS) {
+            skip_ref_ = std::max(skip_ref_, read->cigar_oplen());
+          }
+          ++read;
         }
+
         // stepped into an insertion
-        if (skip_ref_) {
-          for (size_t i=0; i<data_.size(); ++i)
-            data_[i].skip_ = skip_ref_;
-        }
+        if (skip_ref_)
+          for (iterator iter=data_.begin(); iter!=data_.end(); ++iter)
+            iter->skip_ = skip_ref_;
       }
     }
@@ -310 +293 @@
       ++bam_;
     }
-
-    // pop reads that are end of sequence
-    while (data_.size() && data_[0].end_of_sequence()) {
-      data_.pop_front();
-    }
   }
 
@@ -326 +304 @@
 
   template<typename Iterator>
-  Pileup<Iterator>::Entry::Entry(const yat::omic::BamRead& b)
-    : bam_(b), qpos_(0), c_index_(0), c_index2_(0), skip_(0)
-  {
-    if (bam_.core().n_cigar==0)
+  Pileup<Iterator>::Entry::Entry(const BamRead& b)
+    : bam_(new BamRead(b)), qpos_(0), cigar_(bam_->cigar()), skip_(0)
+  {
+    YAT_ASSERT(cigar_.base() == bam_->cigar());
+    if (bam_->core().n_cigar==0)
       return;
     // iterate cigar to first match
-    while (bam_cigar_type(bam_.cigar(c_index_)) != 3) {
-      YAT_ASSERT(bam_cigar_type(bam_.cigar(c_index_))!=2);
-      if (bam_cigar_type(bam_.cigar(c_index_)) == 1) {
-        qpos_ += b.cigar_oplen(c_index_);
-      }
-      ++c_index_;
-      if (c_index_ >= bam_.core().n_cigar) {
-        qpos_ = bam_.sequence_length();
-        return;
-      }
-      YAT_ASSERT(c_index_ < bam_.core().n_cigar);
+    while (qpos_ < bam_->sequence_length() && cigar_type() != 3) {
+      // type is either 0 (hardclipped) or 1 (softclipped)
+      YAT_ASSERT(cigar_type()!=2);
+      if (cigar_type() == 1)
+        ++qpos_;
+      ++cigar_;
     }
   }
@@ -350 +324 @@
   void Pileup<Iterator>::Entry::increment(void)
   {
-    if (bam_cigar_type(cigar()) & 1)
+    YAT_ASSERT(cigar_.base() >= bam_->cigar());
+    if (cigar_type() & 1)
       ++qpos_;
-    step_cigar();
+    ++cigar_;
+    YAT_ASSERT(cigar_.base() >= bam_->cigar());
   }
 
@@ -360 +336 @@
   {
     if (skip_==0) {
-      if (bam_cigar_type(cigar()) & 1)
-        return bam_.sequence(qpos_);
-      else
-        return 0;
+      if (cigar_type() & 1)
+        return bam_->sequence(qpos_);
+      return 0;
     }
     if (cigar_op()!=BAM_CINS)
       return 0;
-    return bam_.sequence(qpos_);
-  }
-
-
-  template<typename Iterator>
-  void Pileup<Iterator>::Entry::step_cigar(void)
-  {
-    ++c_index2_;
-    YAT_ASSERT(c_index_<bam_.core().n_cigar);
-    if (c_index2_ >= bam_.cigar_oplen(c_index_)) {
-      c_index2_ = 0;
-      ++c_index_;
-    }
+    return bam_->sequence(qpos_);
   }