Changeset 3031

Timestamp:

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

Author:

Peter

Message:

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

Location:

trunk

Files:

• test/bam.cc (modified) (2 diffs)
• yat/omic/BamRead.cc (modified) (2 diffs)
• yat/omic/BamRead.h (modified) (2 diffs)

trunk/test/bam.cc

@@ -120 +120 @@
 }
 
+
+void test_sequence(test::Suite& suite, const BamRead& b)
+{
+  suite.out() << "test sequence\n";
+  BamRead bam(b);
+  std::string seq = b.sequence();
+  std::vector<uint8_t> qual;
+  qual.reserve(seq.size());
+  for (size_t i=0; i<seq.size(); ++i)
+    qual.push_back(bam.qual(i));
+  assert(seq.size()==100);
+  bam.sequence(seq, qual);
+
+  if (bam.sequence()!=seq) {
+    suite.add(false);
+    suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
+    suite.err() << "expected: " << seq << "\n";
+  }
+
+  // just for consistency
+  std::vector<uint32_t> cig;
+  bam.cigar(cig);
+  bam.core().flag &= ~BAM_FUNMAP;
+
+  // trim off one base
+  seq.resize(seq.size()-1);
+  qual.resize(qual.size()-1);
+  bam.sequence(seq, qual);
+
+  if (bam.sequence()!=seq) {
+    suite.add(false);
+    suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
+    suite.err() << "expected: " << seq << "\n";
+  }
+
+  // extend sequence
+  seq.resize(120, 'A');
+  qual.resize(120, 0);
+  bam.sequence(seq, qual);
+  if (bam.sequence()!=seq) {
+    suite.add(false);
+    suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
+    suite.err() << "expected: " << seq << "\n";
+  }
+}
+
+
 void test1(test::Suite& suite)
 {
@@ -165 +212 @@
   test_cigar(suite, bam, in.header());
   test_aux(suite, bam);
-}
-#endif
+  test_sequence(suite, bam);
+}
+#endif

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)
   {

trunk/yat/omic/BamRead.h

@@ -289 +289 @@
 
     /**
+       \brief set sequence and quality
+
+       \since New in yat 0.11
+     */
+    void sequence(const std::string& seq, const std::vector<uint8_t>& qual);
+
+    /**
        \see core().l_qseq
      */
@@ -306 +313 @@
 
   private:
+    // ensure capacity of data pointer is (at least) n
+    void reserve(int n);
+
     bam1_t* bam_;