Changeset 3173

Timestamp:

Mar 8, 2014, 2:48:13 PM (9 years ago)

Author:

Peter

Message:

a first version of BamPairIterator?. refs #784

Location:

trunk

Files:

• test/Makefile.am (modified) (1 diff)
• yat/omic/BamPairIterator.h (added)
• yat/omic/Makefile.am (modified) (2 diffs)
• yat/omic/algorithm.h (modified) (2 diffs)

trunk/test/Makefile.am

@@ -46 +46 @@
   test/bam_header.test \
   test/bam_header2.test \
+  test/bam_pair_iterator.test \
   test/codon.test test/commandline.test \
   test/concept.test \

trunk/yat/omic/Makefile.am

@@ -1 +1 @@
 ## $Id$
 
-# Copyright (C) 2010, 2011, 2012, 2013 Peter Johansson
+# Copyright (C) 2010, 2011, 2012, 2013, 2014 Peter Johansson
 #
 # This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -36 +36 @@
 nobase_include_HEADERS += $(srcdir)/yat/omic/BamFile.h
 nobase_include_HEADERS += $(srcdir)/yat/omic/BamHeader.h
+nobase_include_HEADERS += $(srcdir)/yat/omic/BamPairIterator.h
 nobase_include_HEADERS += $(srcdir)/yat/omic/BamRead.h
 nobase_include_HEADERS += $(srcdir)/yat/omic/BamReadFilter.h

trunk/yat/omic/algorithm.h

@@ -23 +23 @@
 */
 
+#include "BamPairIterator.h"
 #include "BamRead.h"
 #include "GenomicPosition.h"
@@ -59 +60 @@
   void bam_pair_analyse(Iterator first, Iterator last, Visitor& visitor)
   {
-    BOOST_CONCEPT_ASSERT((boost::InputIterator<Iterator>));
-
-    using boost::Convertible;
-    typedef typename std::iterator_traits<Iterator>::reference reference_type;
-    BOOST_CONCEPT_ASSERT((Convertible<reference_type, BamRead>));
-
-    // map in which key is mate's position
-    typedef std::multimap<GenomicPosition, BamRead> Map;
-    Map reads;
-    std::map<std::string, BamRead> siam_reads;
-    GenomicPosition siam_position(0,0);
-    for (; first!=last; ++first) {
-      GenomicPosition position(first->core().tid, first->pos());
-      GenomicPosition mate_position(first->core().mtid, first->mpos());
-
-      if (position != siam_position) {
-        siam_position = position;
-        siam_reads.clear();
-      }
-
-      // have not seen the mate yet
-      if (mate_position > position)
-        reads.insert(std::make_pair(mate_position, *first));
-      else if (mate_position < position) {
-        // erase all reads with mate position less than current
-        // position (assuming input range is sorted)
-        Map::iterator lower = reads.lower_bound(position);
-        reads.erase(reads.begin(), lower);
-        // iterator may have been invalidated
-        lower = reads.begin();
-        // find mate and let visitor act on the pair
-        for (; lower!=reads.end() && lower->first == position; ++lower) {
-          if (same_query_name(lower->second, *first)) {
-            visitor(lower->second, *first);
-            break;
-          }
-        }
-      }
-      else { // read and mate have same position
-        // find mate
-        std::map<std::string, BamRead>::iterator
-          mate = siam_reads.lower_bound(first->name());
-        if (mate!=siam_reads.end() && same_query_name(mate->second, *first)) {
-          visitor(mate->second, *first);
-        }
-        else
-          // insert with hint
-          siam_reads.insert(mate, std::make_pair(first->name(), *first));
-      }
-    }
+    BamPairIterator<Iterator> iter(first, last);
+    BamPairIterator<Iterator> end(last, last);
+    for (; iter!=end; ++iter)
+      visitor(iter->first, iter->second);
   }