Changeset 3351

Timestamp:

Nov 21, 2014, 1:20:55 PM (8 years ago)

Author:

Peter

Message:

make BamRead? compile with htslib. refs #794

Location:

trunk/yat/omic

Files:

• BamRead.cc (modified) (16 diffs)
• BamRead.h (modified) (3 diffs)

trunk/yat/omic/BamRead.cc

@@ -25 +25 @@
 #include "config_bam.h"
 
-#include YAT_BAM_HEADER
+#include YAT_SAM_HEADER
 
 #include "yat/utility/Exception.h"
@@ -41 +41 @@
 #endif
 
+// these #defines are no longer available in htslib, define them using
+// htslib version
+#ifndef bam1_seq
+#define bam1_seq(b) (bam_get_seq(b))
+#endif
+
+#ifndef bam1_seqi
+#define bam1_seqi(b, i) (bam_seqi(b, i))
+#endif
+
+#ifndef bam1_cigar
+#define bam1_cigar(b) (bam_get_cigar(b))
+#endif
+
+#ifndef bam1_qname
+#define bam1_qname(b) (bam_get_qname(b))
+#endif
+
+#ifndef bam1_qual
+#define bam1_qual(b) (bam_get_qual(b))
+#endif
+
 namespace theplu {
 namespace yat {
@@ -73 +95 @@
   {
     assert(bam_);
+#if YAT_HAVE_HTSLIB
+    return bam_get_aux(bam_);
+#else
     return bam1_aux(bam_);
+#endif
   }
 
@@ -107 +133 @@
   {
     assert(bam_);
+#if YAT_HAVE_HTSLIB
+    return bam_get_l_aux(bam_);
+#else
     return bam_->l_aux;
+#endif
   }
 
@@ -143 +173 @@
   const uint32_t* BamRead::cigar(void) const
   {
+#if YAT_HAVE_HTSLIB
+    return bam_get_cigar(bam_);
+#else
     return bam1_cigar(bam_);
+#endif
   }
 
@@ -150 +184 @@
   {
     assert(i<core().n_cigar);
-    return bam1_cigar(bam_)[i];
+    return cigar()[i];
   }
 
@@ -172 +206 @@
     std::ostringstream os;
     for (uint32_t i=0; i<core().n_cigar; ++i) {
-      os << bam_cigar_oplen(bam1_cigar(bam_)[i])
-         << bam_cigar_opchr(bam1_cigar(bam_)[i]);
+      os << bam_cigar_oplen(cigar(i))
+         << bam_cigar_opchr(cigar(i));
     }
     return os.str();
@@ -192 +226 @@
   {
     int offset = 4*c.size() - 4*core().n_cigar;
-    reserve(bam_->data_len + offset);
+    
+    reserve(data_size() + offset);
     /*
       data come as qname, cigar, rest
@@ -199 +234 @@
     if (offset)
       memmove(bam1_seq(bam_)+offset, bam1_seq(bam_),
-              bam_->data_len - 4*core().n_cigar - core().l_qname);
+              data_size() - 4*core().n_cigar - core().l_qname);
 
     // copy new cigar
@@ -205 +240 @@
       bam1_cigar(bam_)[i] = c[i];
 
-    bam_->data_len += offset;
+    data_size() += offset;
     core().n_cigar = c.size();
-    assert(bam_->data_len <= bam_->m_data);
+    assert(data_size() <= bam_->m_data);
   }
 
@@ -219 +254 @@
 
 
+  int& BamRead::data_size(void)
+  {
+#if YAT_HAVE_HTSLIB
+    return bam_->l_data;
+#else
+    return bam_->data_len;
+#endif
+  }
+
+
+  const int& BamRead::data_size(void) const
+  {
+#if YAT_HAVE_HTSLIB
+    return bam_->l_data;
+#else
+    return bam_->data_len;
+#endif
+  }
+
+
   int32_t BamRead::end(void) const
   {
@@ -230 +285 @@
     return bam_->core.flag;
   }
+
 
   const char* BamRead::name(void) const
@@ -239 +295 @@
     int offset = n.size() + 1 - core().l_qname;
     assert(bam_);
-    reserve(bam_->data_len + offset);
+    reserve(data_size() + offset);
     // move remaining data
     if (offset)
       memmove(bam_->data + core().l_qname + offset,
               bam_->data + core().l_qname,
-              bam_->data_len - core().l_qname);
+              data_size() - core().l_qname);
     core().l_qname += offset;
     // copy new name
     memcpy(bam1_qname(bam_), n.c_str(), n.size()+1);
-    bam_->data_len += offset;
-    assert(bam_->data_len <= bam_->m_data);
+    data_size() += offset;
+    assert(data_size() <= bam_->m_data);
   }
 
@@ -278 +334 @@
   {
     std::string result(sequence_length(), ' ');
-#ifdef HAVE_BAM_NT16_REV_TABLE
+#if HAVE_BAM_NT16_REV_TABLE
     for (size_t i=0; i<result.size(); ++i)
       result[i] = bam_nt16_rev_table[sequence(i)];
@@ -298 +354 @@
 
     int aux_pos = bam1_seq(bam_) - bam_->data + (seq.size()+1)/2 + qual.size();
-    int len = aux_pos + bam_->l_aux;
+    int len = aux_pos + aux_size();
     reserve(len);
     // move aux to its new location
     memmove(static_cast<void*>(bam_->data+aux_pos),
-            static_cast<const void*>(bam1_aux(bam_)), bam_->l_aux);
+            static_cast<const void*>(aux()), aux_size());
     // copy sequence
     for (size_t i=0; i<seq.size(); ++i)
+#if YAT_HAVE_HTSLIB
+      sequence(i, seq_nt16_table[static_cast<size_t>(seq[i])]);
+#else
       sequence(i, bam_nt16_table[static_cast<size_t>(seq[i])]);
+#endif
     core().l_qseq = seq.size();
     // copy quality
@@ -311 +371 @@
       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);
+    data_size() = len;
+    assert(data_size() <= bam_->m_data);
   }
 

trunk/yat/omic/BamRead.h

@@ -24 +24 @@
 
 #include "config_bam.h"
-#include YAT_SAM_HEADER
 
 #include "yat/utility/Aligner.h"
@@ -31 +30 @@
 #include "yat/utility/deprecate.h"
 
+#include YAT_SAM_HEADER
+
 #include <functional>
 #include <string>
@@ -309 +310 @@
     friend class BamReadIterator;
     uint32_t calend(const bam1_core_t *c, const uint32_t *cigar) const;
+    // access data length, current length of data
+    int& data_size(void);
+    const int& data_size(void) const;
   };