Changeset 3031 for trunk/yat/omic/BamRead.cc

Timestamp:

Apr 27, 2013, 1:24:13 PM (10 years ago)

Author:

Peter

Message:

new function BamRead::sequence(2). refs #746

File:

• trunk/yat/omic/BamRead.cc (modified) (2 diffs)

trunk/yat/omic/BamRead.cc

@@ -206 +206 @@
 
 
+  void BamRead::reserve(int n)
+  {
+    assert(bam_);
+    if (bam_->m_data >= n)
+      return;
+    // at least double capacity
+    bam_->m_data = std::max(bam_->m_data<<1, n);
+    bam_->data = (uint8_t*)realloc(bam_->data, bam_->m_data);
+    if (!bam_->data) {
+      throw utility::runtime_error("BamRead::reserve failed");
+    }
+  }
+
+
   std::string BamRead::sequence(void) const
   {
@@ -221 +235 @@
 
 
+  void BamRead::sequence(const std::string& seq,
+                         const std::vector<uint8_t>& qual)
+  {
+    assert(seq.size()==qual.size());
+    // data consist of
+    // |-- qname -->-- cigar -->-- seq -->-- qual -->-- aux -->
+
+    uint8_t* new_aux = bam1_seq(bam_) + (seq.size()+1)/2 + qual.size();
+
+    int len = new_aux + bam_->l_aux - bam_->data;
+    reserve(len);
+    // move aux to its new location
+    memmove(static_cast<void*>(new_aux),
+            static_cast<const void*>(bam1_aux(bam_)), bam_->l_aux);
+    // copy sequence
+    for (size_t i=0; i<seq.size(); ++i)
+      sequence(i, bam_nt16_table[seq[i]]);
+    core().l_qseq = seq.size();
+    // copy quality
+    if (qual.size()) // FIXME is this needed or does memcpy work with n=0?
+      memcpy(static_cast<void*>(bam1_qual(bam_)),
+             static_cast<const void*>(&qual[0]), qual.size());
+  }
+
+
   void BamRead::sequence(size_t i, uint8_t seq)
   {