Changeset 3306

Timestamp:

Aug 21, 2014, 6:37:21 AM (9 years ago)

Author:

Peter

Message:

merge ptach release 0.12.1 into trunk

Location:

trunk

Files:

• . (modified) (1 prop)
• NEWS (modified) (1 diff)
• configure.ac (modified) (1 diff)
• m4/version.m4 (modified) (1 diff)
• m4/yat_check_libbam.m4 (modified) (3 diffs)
• test/bam.cc (modified) (2 diffs)
• yat/omic/BamRead.cc (modified) (3 diffs)
• yat/omic/BamRead.h (modified) (1 diff)
• yat/utility/Cigar.h (modified) (1 diff)

trunk/NEWS

@@ -12 +12 @@
 
 yat 0.12.x series from http://dev.thep.lu.se/yat/svn/branches/0.12-stable
+
+version 0.12.1 (released 21 August 2014)
+  - bam_calend is buggy in older versions of samtools; configure
+    script checks that bam_calend works as expected; if not use its
+    own implementation (see bug #810)
+
+  A complete list of closed tickets can be found here [[br]]
+  http://dev.thep.lu.se/yat/query?status=closed&milestone=yat+0.12.1
 
 version 0.12 (released 25 July 2014)

trunk/configure.ac

@@ -383 +383 @@
   # check if global variable bam_nt16_rev_table is available
   YAT_BAM_NT16_REV_TABLE
+  # check if function bam_calend works. Result of this test can be
+  # overridden by setting variable 'yat_cv_func_bam_calend'
+  YAT_FUNC_BAM_CALEND([$srcdir/test/data/foo.sam])
   AC_PATH_PROG([SAMTOOLS], [samtools], [false])
   AC_ARG_VAR([SAMTOOLS], [Tools for alignment in the SAM format])

trunk/m4/version.m4

@@ -77 +77 @@
 # yat-0.11.2 8:2:0
 # yat-0.12   9:0:0
+# yat-0.12.1 9:1:0
 #
 # *Accidently, the libtool number was not updated for yat 0.5

trunk/m4/yat_check_libbam.m4

@@ -1 +1 @@
 ## $Id$
 #
-# serial 4 (yat 0.11)
+# serial 5 (yat 0.12.1)
 #
 #
-#   Copyright (C) 2012, 2013 Peter Johansson
+#   Copyright (C) 2012, 2013, 2014 Peter Johansson
 #
 #   This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -107 +107 @@
 
 
+# YAT_FUNC_BAM_CALEND(SAM_FILE)
+# ===================
+# Check if there is a bam_calend that works
+AC_DEFUN([YAT_FUNC_BAM_CALEND],
+[
+  AC_CACHE_CHECK([for working bam_calend],
+    [yat_cv_func_bam_calend],
+    [AC_RUN_IFELSE(
+       [AC_LANG_PROGRAM(
+          [_YAT_BAM_INCLUDES],
+          [
+          samfile_t* samfile = samopen("$1", "r", NULL);
+          if (!samfile)
+            return 1;
+          bam1_t* read = bam_init1();
+          while (samread(samfile, read) >= -1) {
+            bool have_match = false;
+            if (read->core.n_cigar < 1)
+              continue;
+            uint32_t end = bam_calend(&read->core, bam1_cigar(read));
+            // replace BAM_CMATCH with BAM_CEQUAL in cigar
+            for (size_t k=0; k<read->core.n_cigar; ++k) {
+              uint32_t& element = bam1_cigar(read)@<:@k@:>@;
+              uint32_t op = element & BAM_CIGAR_MASK;
+              if (op == BAM_CMATCH) {
+                have_match = true;
+                uint32_t oplen = element >> BAM_CIGAR_SHIFT;
+                element = oplen<<BAM_CIGAR_SHIFT|BAM_CEQUAL;
+              }
+            }
+            // if no match element found, next read please
+            if (!have_match)
+              continue;
+            // bam_calend works on BAM_CEQUAL
+            if (bam_calend(&read->core, bam1_cigar(read)) == end)
+              return 0;
+            // bam_calend is not working as expected
+            return 2;
+          }
+          // end of data, fail
+          return 3;
+          ])
+       ],
+       [yat_cv_func_bam_calend=yes],
+       [yat_cv_func_bam_calend=no],
+       [yat_cv_func_bam_calend="guessing no"])
+    ])
+  AS_IF([test x"$yat_cv_func_bam_calend" = x"yes"],[
+    AC_DEFINE([HAVE_BAM_CALEND], [1], [Define to 1 if bam_calend is working])
+  ])
+])
+
+
 # _YAT_BAM_INCLUDES
 # =================
@@ -114 +167 @@
 @%:@if YAT_HAVE_BAM_H
 @%:@ include <bam.h>
+@%:@ include <sam.h>
 @%:@elif YAT_HAVE_BAM_BAM_H
 @%:@ include <bam/bam.h>
+@%:@ include <bam/sam.h>
 @%:@elif YAT_HAVE_SAMTOOLS_BAM_H
 @%:@ include <samtools/bam.h>
+@%:@ include <samtools/sam.h>
 @%:@endif
 ]) # _YAT_BAM_HEADER

trunk/test/bam.cc

@@ -120 +120 @@
 
 
+void test_end(test::Suite& suite, BamRead bam)
+{
+  suite.out() << "test end:\n";
+  suite.out() << bam.pos() << "-" << bam.end() << " " << bam.cigar_str()
+              << "\n";
+  uint32_t end = bam.end();
+  if (!suite.add(end==24093)) {
+    suite.err() << "error: incorrect end: " << end << "\n";
+  }
+
+  // change the cigar to include a DIFF and EQUAL
+  utility::Aligner::Cigar cig;
+  cig.push_back(BAM_CEQUAL, 55);
+  cig.push_back(BAM_CDIFF, 1);
+  cig.push_back(BAM_CEQUAL, 44);
+  bam.cigar(cig);
+  suite.out() << bam.pos() << "-" << bam.end() << " " << bam.cigar_str()
+              << "\n";
+  end = bam.end();
+  if (!suite.add(end==24093)) {
+    suite.err() << "error: incorrect end: " << end << "\n";
+  }
+}
+
+
 void test_sequence(test::Suite& suite, const BamRead& b)
 {
@@ -234 +259 @@
   in.read(bam);
   // test that functions are implemented
-  bam.end();
+  test_end(suite, bam);
   bam.flag();
   bam.mpos();

trunk/yat/omic/BamRead.cc

@@ -107 +107 @@
     assert(bam_);
     return bam_->l_aux;
+  }
+
+
+  uint32_t BamRead::calend(const bam1_core_t *c, const uint32_t *cigar) const
+  {
+#ifdef HAVE_BAM_CALEND
+    return bam_calend(&core(), bam1_cigar(bam_));
+#endif
+    int k, end = c->pos;
+    for (k = 0; k < c->n_cigar; ++k) {
+      int op  = bam_cigar_op(cigar[k]);
+      int len = bam_cigar_oplen(cigar[k]);
+      if (op == BAM_CBACK) { // move backward
+        int l, u, v;
+        if (k == c->n_cigar - 1) break; // skip trailing 'B'
+        for (l = k - 1, u = v = 0; l >= 0; --l) {
+          int op1  = bam_cigar_op(cigar[l]);
+          int len1 = bam_cigar_oplen(cigar[l]);
+          if (bam_cigar_type(op1)&1) { // consume query
+            if (u + len1 >= len) { // stop
+              if (bam_cigar_type(op1)&2) v += len - u;
+              break;
+            } else u += len1;
+          }
+          if (bam_cigar_type(op1)&2) v += len1;
+        }
+        end = l < 0? c->pos : end - v;
+      } else if (bam_cigar_type(op)&2) end += bam_cigar_oplen(cigar[k]);
+    }
+    return end;
   }
 
@@ -186 +216 @@
   bam1_core_t& BamRead::core(void)
   { return bam_->core; }
+
+
+  int32_t BamRead::end(void) const
+  {
+    assert(bam_);
+    return this->calend(&core(), bam1_cigar(bam_));
+  }
 
 
@@ -328 +365 @@
 
 
-  int32_t BamRead::end(void) const
-  {
-    assert(bam_);
-    return bam_calend(&core(), bam1_cigar(bam_));
-  }
-
-
   void BamRead::swap(BamRead& other)
   {

trunk/yat/omic/BamRead.h

@@ -309 +309 @@
     friend class OutBamFile;
     friend class BamReadIterator;
+    uint32_t calend(const bam1_core_t *c, const uint32_t *cigar) const;
   };
 

trunk/yat/utility/Cigar.h

@@ -56 +56 @@
 /// \abstract CIGAR: X = mismatch
 #define BAM_CDIFF       8
-#define BAM_CBACK       9
 #endif // end of YAT_HAVE_LIBBAM
+
+// BAM_CBACK is not defined in old versions of bam.h
+#ifndef BAM_CBACK
+#define BAM_CBACK 9
+#endif
 
 // backport #defines from samtools 0.1.19