source: branches/0.12-stable/yat/omic/BamPairIterator.h @ 3234

Last change on this file since 3234 was 3234, checked in by Peter, 9 years ago

fixes #797

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 6.9 KB
Line 
1#ifndef theplu_yat_omic_bam_pair_iterator
2#define theplu_yat_omic_bam_pair_iterator
3
4// $Id: BamPairIterator.h 3234 2014-05-23 10:08:48Z 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 "BamPair.h"
26#include "BamRead.h"
27
28#include <yat/utility/config_public.h>
29#include <yat/utility/yat_assert.h>
30
31#include <boost/concept/assert.hpp>
32#include <boost/iterator/iterator_concepts.hpp>
33#include <boost/iterator/iterator_facade.hpp>
34#include <boost/shared_ptr.hpp>
35
36#include <iterator>
37#include <map>
38#include <utility>
39
40namespace theplu {
41namespace yat {
42namespace omic {
43
44  /**
45     \since New in yat 0.12
46
47     Type Requirements:
48     - \c Iterator must be an \input_iterator
49     - \c Iterator 's \c reference type must be convertible to BamRead
50
51     This iterator works on a sorted range [begin, end) and provies a
52     covenient way to access the pairs rather than the reads
53     individually. When BamPairIterator is incremented, rather than
54     just step to next element in [begin, end), BamPairIterator keeps
55     stepping until it finds a read whose mate it has already seen, and
56     operator* return the pair of BamReads.
57
58     Note that BamPairIterator is a single-pass iterator, i.e., once it
59     is incremented the behaviour of its copies is undefined.
60   */
61  template<typename Base>
62  class BamPairIterator
63    : public boost::iterator_facade<
64    BamPairIterator<Base>, const BamPair, std::input_iterator_tag,
65    const BamPairProxy
66    >
67  {
68  public:
69    /**
70       \brief default constructor
71     */
72    BamPairIterator(void);
73
74    /**
75       Creates an iterator that will work on [begin, end).
76
77       \see bam_pair_iterator
78
79       \note Input range \c [\a first, \a last \c ) must be sorted or
80       behaviour is undefined.
81     */
82    explicit BamPairIterator(Base begin, Base end);
83
84#ifndef YAT_HAVE_BOOST_ITERATOR_FACADE_PROXY_PTR
85    /**
86       This is workaround implementation of operator-> when
87       implementation in base class (boost) does not compile. This
88       implementation may be slow, so when using old boost it often
89       prefereble to use operator*.
90     */
91    typename BamPairIterator::pointer operator->(void) const
92    {
93      // older versions of boost mandates pointer to be
94      // value_type*. To accomplish that we implement this slow thing
95      // which keeps a copy of a value_type as member.
96      // See https://svn.boost.org/trac/boost/ticket/5697 for discussion.
97
98      // Possibly stupid to assign each time, but why bother optimize
99      // for old bugs in boost. Users who are eager for speed, should
100      // either upgrade their boost or use (*iterator).function()
101      // rather than iterator->function().
102      YAT_ASSERT(mate_);
103      YAT_ASSERT(iter_ != end_);
104      dummie_.first() = *mate_;
105      dummie_.second() = *iter_;
106      return &dummie_;
107    }
108  private:
109    mutable BamPair dummie_;
110#endif
111
112  private:
113    Base iter_;
114    Base end_;
115    const BamRead* mate_;
116    boost::shared_ptr<std::map<std::string, BamRead> > siam_reads_;
117    typedef std::pair<int32_t, int32_t> Position;
118    boost::shared_ptr<std::multimap<Position, BamRead> > reads_;
119    friend class boost::iterator_core_access;
120
121    const BamPairProxy dereference(void) const;
122    bool equal(const BamPairIterator& other) const;
123    void increment(void);
124    void find_next(void);
125  };
126
127
128  /**
129     Convenient function to create BamPairIterator
130
131     \relates BamPairIterator
132   */
133  template<typename Base>
134  BamPairIterator<Base> bam_pair_iterator(Base base, Base end)
135  { return BamPairIterator<Base>(base, end); }
136
137
138  /**
139     Convenient function to create BamPairIterator that works on empty
140     range. Equivalent to bam_pair_iterator(end, end).
141
142     \relates BamPairIterator
143  */
144  template<typename Base>
145  BamPairIterator<Base> bam_pair_iterator(Base end)
146  { return BamPairIterator<Base>(end, end); }
147
148
149  //             template implementations
150  ///////////////////////////////////////////////////////////////
151
152  template<typename Base>
153  BamPairIterator<Base>::BamPairIterator(Base base, Base end)
154    : iter_(base), end_(end),
155      siam_reads_(new std::map<std::string, BamRead>),
156      reads_(new std::multimap<Position, BamRead>)
157  {
158    BOOST_CONCEPT_ASSERT((boost::InputIterator<Base>));
159    using boost::Convertible;
160    typedef typename std::iterator_traits<Base>::reference reference_type;
161    BOOST_CONCEPT_ASSERT((Convertible<reference_type, BamRead>));
162    find_next();
163  }
164
165
166  template<typename Base>
167  BamPairIterator<Base>::BamPairIterator(void)
168  {
169    BOOST_CONCEPT_ASSERT((boost::InputIterator<Base>));
170    using boost::Convertible;
171    typedef typename std::iterator_traits<Base>::reference reference_type;
172    BOOST_CONCEPT_ASSERT((Convertible<reference_type, BamRead>));
173  }
174
175
176  template<typename Base>
177  const BamPairProxy
178  BamPairIterator<Base>::dereference(void) const
179  {
180    return BamPairProxy(mate_, &*iter_);
181  }
182
183
184  template<typename Base>
185  bool BamPairIterator<Base>::equal(const BamPairIterator& other) const
186  {
187    return iter_ == other.iter_;
188  }
189
190
191  template<typename Base>
192  void BamPairIterator<Base>::increment(void)
193  {
194    ++iter_;
195    find_next();
196  }
197
198
199  template<typename Base>
200  void BamPairIterator<Base>::find_next(void)
201  {
202    BamLessPos less;
203    for (; iter_!=end_; ++iter_) {
204      Position position(iter_->tid(), iter_->pos());
205      Position mate_position(iter_->mtid(), iter_->mpos());
206      // clear siam reads if iter is more advanced than siam reads
207      if (siam_reads_->size() && less(siam_reads_->begin()->second, *iter_))
208        siam_reads_->clear();
209
210      // have not seen mate yet
211      if (mate_position > position)
212        reads_->insert(std::make_pair(mate_position, *iter_));
213      else if (position > mate_position) {
214        std::multimap<Position, BamRead>::iterator
215          lower = reads_->lower_bound(position);
216        // erase all reads with mate position less than current
217        // position (assuming input range is sorted)
218        reads_->erase(reads_->begin(), lower);
219
220        // find mate and assign pair
221        for (; lower!=reads_->end() && lower->first == position; ++lower)
222          if (same_query_name(lower->second, *iter_)) {
223            mate_ = &lower->second;
224            return;
225          }
226      }
227      else { // siam read, i.e., iter_ and mate have same position
228        // check if we have already seen mate and stored it in map
229        std::map<std::string, BamRead>::iterator
230          mate = siam_reads_->lower_bound(iter_->name());
231        if (mate!=siam_reads_->end() && same_query_name(mate->second, *iter_)) {
232          mate_ = &mate->second;
233          return;
234        }
235        else
236          // insert with hint
237          siam_reads_->insert(mate, std::make_pair(iter_->name(), *iter_));
238      }
239    }
240  }
241
242}}}
243#endif
Note: See TracBrowser for help on using the repository browser.