source: trunk/yat/omic/BamPairIterator.h @ 3176

Last change on this file since 3176 was 3176, checked in by Peter, 8 years ago

closes #784. New class BamPairProxy? to finalize implementation of BamPairIterator?.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 5.8 KB
Line 
1#ifndef theplu_yat_omic_bam_pair_iterator
2#define theplu_yat_omic_bam_pair_iterator
3
4// $Id: BamPairIterator.h 3176 2014-03-15 09:10:37Z peter $
5
6/*
7  Copyright (C) 2014 Peter Johansson
8
9  This file is part of the yat library, http://dev.thep.lu.se/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 3 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with this program. If not, see <http://www.gnu.org/licenses/>.
23*/
24
25#include <iostream>
26
27#include "BamPair.h"
28#include "BamRead.h"
29
30#include <boost/concept/assert.hpp>
31#include <boost/iterator/iterator_facade.hpp>
32#include <boost/shared_ptr.hpp>
33
34#include <iterator>
35#include <map>
36#include <utility>
37
38namespace theplu {
39namespace yat {
40namespace omic {
41
42  /**
43     \since New in yat 0.12
44
45     Type Requirements:
46     - \c Iterator must be an \input_iterator
47     - \c Iterator 's \c reference type must be convertible to BamRead
48
49     This iterator works on a sorted range [begin, end) and provies a
50     covenient way to access the pairs rather than the reads
51     individually. When BamPairIterator is incremented, rather than
52     just step to next element in [begin, end), BamPairIterator keeps
53     stepping until it finds a read whose mate it has already seen, and
54     operator* return the pair of BamReads.
55
56     Note that BamPairIterator is a single-pass iterator, i.e., once it
57     is incremented the behaviour of its copies is undefined.
58   */
59  template<typename Base>
60  class BamPairIterator
61    : public boost::iterator_facade<
62    BamPairIterator<Base>, const BamPair, std::input_iterator_tag,
63    const BamPairProxy
64    >
65  {
66  public:
67    /**
68       \brief default constructor
69     */
70    BamPairIterator(void);
71
72    /**
73       Creates an iterator that will work on [begin, end).
74
75       \see bam_pair_iterator
76
77       \note Input range \c [\a first, \a last \c ) must be sorted or
78       behaviour is undefined.
79     */
80    explicit BamPairIterator(Base begin, Base end);
81  private:
82    Base iter_;
83    Base end_;
84    const BamRead* mate_;
85    boost::shared_ptr<std::map<std::string, BamRead> > siam_reads_;
86    typedef std::pair<int32_t, int32_t> Position;
87    boost::shared_ptr<std::multimap<Position, BamRead> > reads_;
88    friend class boost::iterator_core_access;
89
90    const BamPairProxy dereference(void) const;
91    bool equal(const BamPairIterator& other) const;
92    void increment(void);
93    void find_next(void);
94  };
95
96
97  /**
98     Convenient function to create BamPairIterator
99
100     \relates BamPairIterator
101   */
102  template<typename Base>
103  BamPairIterator<Base> bam_pair_iterator(Base base, Base end)
104  { return BamPairIterator<Base>(base, end); }
105
106
107  /**
108     Convenient function to create BamPairIterator that works on empty
109     range. Equivalent to bam_pair_iterator(end, end).
110
111     \relates BamPairIterator
112  */
113  template<typename Base>
114  BamPairIterator<Base> bam_pair_iterator(Base end)
115  { return BamPairIterator<Base>(end, end); }
116
117
118  //             template implementations
119  ///////////////////////////////////////////////////////////////
120
121  template<typename Base>
122  BamPairIterator<Base>::BamPairIterator(Base base, Base end)
123    : iter_(base), end_(end),
124      siam_reads_(new std::map<std::string, BamRead>),
125      reads_(new std::multimap<Position, BamRead>)
126  {
127    BOOST_CONCEPT_ASSERT((boost::InputIterator<Base>));
128    using boost::Convertible;
129    typedef typename std::iterator_traits<Base>::reference reference_type;
130    BOOST_CONCEPT_ASSERT((Convertible<reference_type, BamRead>));
131    find_next();
132  }
133
134
135  template<typename Base>
136  BamPairIterator<Base>::BamPairIterator(void)
137  {
138    BOOST_CONCEPT_ASSERT((boost::InputIterator<Base>));
139    using boost::Convertible;
140    typedef typename std::iterator_traits<Base>::reference reference_type;
141    BOOST_CONCEPT_ASSERT((Convertible<reference_type, BamRead>));
142  }
143
144
145  template<typename Base>
146  const BamPairProxy
147  BamPairIterator<Base>::dereference(void) const
148  {
149    return BamPairProxy(mate_, &*iter_);
150  }
151
152
153  template<typename Base>
154  bool BamPairIterator<Base>::equal(const BamPairIterator& other) const
155  {
156    return iter_ == other.iter_;
157  }
158
159
160  template<typename Base>
161  void BamPairIterator<Base>::increment(void)
162  {
163    ++iter_;
164    find_next();
165  }
166
167
168  template<typename Base>
169  void BamPairIterator<Base>::find_next(void)
170  {
171    BamLessPos less;
172    for (; iter_!=end_; ++iter_) {
173      Position position(iter_->tid(), iter_->pos());
174      Position mate_position(iter_->mtid(), iter_->mpos());
175      // clear siam reads if iter is more advanced than siam reads
176      if (siam_reads_->size() && less(siam_reads_->begin()->second, *iter_))
177        siam_reads_->clear();
178
179      // have not seen mate yet
180      if (mate_position > position)
181        reads_->insert(std::make_pair(mate_position, *iter_));
182      else if (position > mate_position) {
183        std::multimap<Position, BamRead>::iterator
184          lower = reads_->lower_bound(position);
185        // erase all reads with mate position less than current
186        // position (assuming input range is sorted)
187        reads_->erase(reads_->begin(), lower);
188
189        // find mate and assign pair
190        for (; lower!=reads_->end() && lower->first == position; ++lower)
191          if (same_query_name(lower->second, *iter_)) {
192            mate_ = &lower->second;
193            return;
194          }
195      }
196      else { // siam read, i.e., iter_ and mate have same position
197        // check if we have already seen mate and stored it in map
198        std::map<std::string, BamRead>::iterator
199          mate = siam_reads_->lower_bound(iter_->name());
200        if (mate!=siam_reads_->end() && same_query_name(mate->second, *iter_)) {
201          mate_ = &mate->second;
202          return;
203        }
204        else
205          // insert with hint
206          siam_reads_->insert(mate, std::make_pair(iter_->name(), *iter_));
207      }
208    }
209  }
210
211}}}
212#endif
Note: See TracBrowser for help on using the repository browser.