Changeset 3032

Timestamp:

Apr 27, 2013, 5:06:52 PM (10 years ago)

Author:

Peter

Message:

refs #746. Simplify cigar function using logic similar to sequence(2). Avoid bug as pointers to data get invalidated when calling reserve. Forgot to set data_len and avoid compiler warning using char as index to array

File:

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

trunk/yat/omic/BamRead.cc

@@ -151 +151 @@
   void BamRead::cigar(const std::vector<uint32_t>& c)
   {
-    // special case when cigar size is not changed
-    if (c.size() == core().n_cigar) {
-      if (c.size())
-        memcpy(bam1_cigar(bam_), &c[0], 4*c.size());
-      return;
-    }
-
-    int m_data = bam_->m_data;
-    int data_len = bam_->data_len + 4*c.size() - 4*core().n_cigar;
-    // double capacity if needed
-    if (m_data < data_len) {
-      m_data = data_len;
-      kroundup32(m_data);
-    }
-    uint8_t* data = (uint8_t*)calloc(m_data, 1);
+    int offset = 4*c.size() - 4*core().n_cigar;
+    reserve(bam_->data_len + offset);
     /*
       data comes as qname, cigar, rest
      */
-    // first copy qname
-    memcpy(data, bam_->data, core().l_qname);
+    // first move remaining data to new location
+    if (offset)
+      memmove(bam1_seq(bam_)+offset, bam1_seq(bam_),
+              bam_->data_len - 4*core().n_cigar - core().l_qname);
+
     // copy new cigar
     if (c.size())
       // each cigar element is 4 bytes
-      memcpy(data + core().l_qname, &c[0], 4*c.size());
-    // copy remaining data
-    memcpy(data + core().l_qname + 4*c.size(), bam1_seq(bam_),
-           data_len - 4*c.size() - core().l_qname);
-    bam_->m_data = m_data;
-    bam_->data_len = data_len;
+      memcpy(bam1_cigar(bam_), &c[0], 4*c.size());
+    bam_->data_len += offset;
     core().n_cigar = c.size();
-    free(bam_->data);
-    bam_->data = data;
+    assert(bam_->data_len <= bam_->m_data);
   }
 
@@ -242 +227 @@
     // |-- 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;
+    int aux_pos = bam1_seq(bam_) - bam_->data + (seq.size()+1)/2 + qual.size();
+    int len = aux_pos + bam_->l_aux;
     reserve(len);
     // move aux to its new location
-    memmove(static_cast<void*>(new_aux),
+    memmove(static_cast<void*>(bam_->data+aux_pos),
             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]]);
+      sequence(i, bam_nt16_table[static_cast<size_t>(seq[i])]);
     core().l_qseq = seq.size();
     // copy quality
@@ -257 +241 @@
       memcpy(static_cast<void*>(bam1_qual(bam_)),
              static_cast<const void*>(&qual[0]), qual.size());
+    bam_->data_len = len;
+    assert(bam_->data_len <= bam_->m_data);
   }