source: branches/0.11-stable/yat/omic/algorithm.h @ 3233

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

refs #797. Fix the bug in 0.11.x

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 3.6 KB
Line 
1#ifndef theplu_yat_omic_algorithm
2#define theplu_yat_omic_algorithm
3
4// $Id: algorithm.h 3233 2014-05-23 09:47:19Z peter $
5
6/*
7  Copyright (C) 2012, 2013, 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 "BamRead.h"
26#include "GenomicPosition.h"
27
28#include <boost/concept/assert.hpp>
29#include <boost/iterator/iterator_concepts.hpp>
30
31#include <iterator>
32#include <map>
33#include <string>
34#include <utility>
35
36namespace theplu {
37namespace yat {
38namespace omic {
39
40  /**
41     \c bam_pair_analyse performs an operation on bam read pairs as
42     defined by \a visitor. The function iterates over sorted input
43     range of reads; if read is first read, it is cached for later
44     use; if read is second read and mate is present in cache, \a
45     visitor operates on pair, i.e., \c Visitor (\c mate, \c read) is
46     called.
47
48     Type Requirements:
49     - \c Iterator must be an \input_iterator
50     - \c Iterator 's \c reference type must be convertible to BamRead
51     - \c Visitor must have an \c operator()(BamRead, BamRead) (or any
52       \c const or reference combination)
53
54     \note Input range \c [\a first, \a last \c ) must be sorted or behaviour is
55     undefined.
56
57     \since New in yat 0.10
58   */
59  template<class Iterator, class Visitor>
60  void bam_pair_analyse(Iterator first, Iterator last, Visitor& visitor)
61  {
62    BOOST_CONCEPT_ASSERT((boost::InputIterator<Iterator>));
63
64    using boost::Convertible;
65    typedef typename std::iterator_traits<Iterator>::reference reference_type;
66    BOOST_CONCEPT_ASSERT((Convertible<reference_type, BamRead>));
67
68    // map in which key is mate's position
69    typedef std::multimap<GenomicPosition, BamRead> Map;
70    Map reads;
71    std::map<std::string, BamRead> siam_reads;
72    GenomicPosition siam_position(0,0);
73    for (; first!=last; ++first) {
74      GenomicPosition position(first->core().tid, first->pos());
75      GenomicPosition mate_position(first->core().mtid, first->mpos());
76
77      if (position != siam_position) {
78        siam_position = position;
79        siam_reads.clear();
80      }
81
82      // have not seen the mate yet
83      if (mate_position > position)
84        reads.insert(std::make_pair(mate_position, *first));
85      else if (mate_position < position) {
86        // erase all reads with mate position less than current
87        // position (assuming input range is sorted)
88        Map::iterator lower = reads.lower_bound(position);
89        reads.erase(reads.begin(), lower);
90        // iterator may have been invalidated
91        lower = reads.begin();
92        // find mate and let visitor act on the pair
93        for (; lower!=reads.end() && lower->first == position; ++lower) {
94          if (same_query_name(lower->second, *first)) {
95            visitor(lower->second, *first);
96            break;
97          }
98        }
99      }
100      else { // read and mate have same position
101        // find mate
102        std::map<std::string, BamRead>::iterator
103          mate = siam_reads.lower_bound(first->name());
104        if (mate!=siam_reads.end() && same_query_name(mate->second, *first)) {
105          visitor(mate->second, *first);
106        }
107        else
108          // insert with hint
109          siam_reads.insert(mate, std::make_pair(first->name(), *first));
110      }
111    }
112  }
113
114}}}
115#endif
Note: See TracBrowser for help on using the repository browser.