Changeset 3200 for trunk/yat/utility/Aligner.cc

Timestamp:

May 3, 2014, 1:36:31 PM (9 years ago)

Author:

Peter

Message:

implement first version of CIGAR (refs #785). Note, this revision does not compile without samtools (should be fixed shortly).

File:

• trunk/yat/utility/Aligner.cc (modified) (3 diffs)

trunk/yat/utility/Aligner.cc

@@ -26 +26 @@
 #include "Matrix.h"
 
+
+// FIXME fix this so it works also --without-samtools
+#include <yat/omic/config_bam.h>
+#include YAT_BAM_HEADER
+
+// FIXME we should include yat/omic/BamRead.h so we get #defines for
+// e.g. bam_cigar_op also if not included in bam.h
+
 #include <cassert>
 #include <limits>
+
+#include <iostream> // debug
 
 namespace theplu {
@@ -118 +128 @@
 
 
+  const Aligner::Cigar Aligner::cigar(size_t row, size_t column) const
+  {
+    Cigar cig;
+    while (row || column) {
+      switch (alignment(row, column)) {
+      case(diagonal):
+        assert(row);
+        --row;
+        assert(column);
+        --column;
+        cig.push_front(BAM_CMATCH);
+        break;
+      case (right):
+        assert(column);
+        --column;
+        cig.push_front(BAM_CINS);
+        break;
+      case (down):
+        assert(row);
+        --row;
+        cig.push_front(BAM_CDEL);
+        break;
+      default:
+        return cig;
+      }
+    }
+    return cig;
+  }
+
+
   Aligner::direction& Aligner::directions(size_t i, size_t j)
   {
@@ -123 +163 @@
   }
 
+
+  // ================ Cigar =======================
+
+  void Aligner::Cigar::clear(void)
+  {
+    cigar_.clear();
+  }
+
+
+  uint8_t Aligner::Cigar::op(size_t i) const
+  {
+    assert(i<cigar_.size());
+    return bam_cigar_op(cigar_[i]);
+  }
+
+
+  char Aligner::Cigar::opchr(size_t i) const
+  {
+    assert(i<cigar_.size());
+    return bam_cigar_opchr(cigar_[i]);
+  }
+
+
+  uint32_t Aligner::Cigar::oplen(size_t i) const
+  {
+    assert(i<cigar_.size());
+    return bam_cigar_oplen(cigar_[i]);
+  }
+
+
+  void Aligner::Cigar::pop_back(uint32_t len)
+  {
+    while (len) {
+      assert(size());
+      size_t i = size()-1;
+      if (len < oplen(i)) {
+        cigar_[i] -= (len<<BAM_CIGAR_SHIFT);
+        return;
+      }
+      else {
+        len -= oplen(i);
+        cigar_.pop_back();
+      }
+    }
+  }
+
+
+  void Aligner::Cigar::pop_front(uint32_t len)
+  {
+    while (len) {
+      assert(size());
+      if (len < oplen(0)) {
+        cigar_[0] -= (len<<BAM_CIGAR_SHIFT);
+        return;
+      }
+      else {
+        len -= oplen(0);
+        cigar_.pop_front();
+      }
+    }
+  }
+
+
+  void Aligner::Cigar::push_back(uint8_t op, uint32_t len)
+  {
+    if (cigar_.empty() || this->op(size()-1)!=op)
+      cigar_.push_front(bam_cigar_gen(len, op));
+    else
+      cigar_[size()-1] += (len<<BAM_CIGAR_SHIFT);
+  }
+
+
+  void Aligner::Cigar::push_front(uint8_t op, uint32_t len)
+  {
+    if (cigar_.empty() || this->op(0)!=op)
+      cigar_.push_front(bam_cigar_gen(len, op));
+    else
+      cigar_[0] += (len<<BAM_CIGAR_SHIFT);
+  }
+
+
+  size_t Aligner::Cigar::size(void) const
+  {
+    return cigar_.size();
+  }
+
+
+  uint32_t Aligner::Cigar::operator[](size_t i) const
+  {
+    assert(i<cigar_.size());
+    return cigar_[i];
+  }
+
+
+  std::ostream& operator<<(std::ostream& os, const Aligner::Cigar& cigar)
+  {
+    for (size_t i=0; i<cigar.size(); ++i)
+      os << cigar.oplen(i) << cigar.opchr(i);
+    return os;
+  }
+
 }}} // of namespace utility, yat, and theplu