source: trunk/yat/omic/algorithm.h @ 2905

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

remove superfluous 'using'

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 3.5 KB
Line 
1#ifndef theplu_yat_omic_algorithm
2#define theplu_yat_omic_algorithm
3
4// $Id: algorithm.h 2905 2012-12-14 00:03:10Z peter $
5//
6// Copyright (C) 2012 Peter Johansson
7//
8// This program is free software; you can redistribute it and/or modify
9// it under the terms of the GNU General Public License as published by
10// the Free Software Foundation; either version 3 of the License, or
11// (at your option) any later version.
12//
13// This program is distributed in the hope that it will be useful, but
14// WITHOUT ANY WARRANTY; without even the implied warranty of
15// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16// General Public License for more details.
17//
18// You should have received a copy of the GNU General Public License
19// along with this program. If not, see <http://www.gnu.org/licenses/>.
20
21#include "BamRead.h"
22#include "GenomicPosition.h"
23
24#include <boost/concept/assert.hpp>
25
26#include <iterator>
27#include <map>
28#include <string>
29#include <utility>
30
31namespace theplu {
32namespace yat {
33namespace omic {
34
35  /**
36     \c bam_pair_analyse performs an operation on bam read pairs as
37     defined by \a visitor. The function iterates over sorted input
38     range of reads; if read is first read, it is cached for later
39     use; if read is second read and mate is present in cache, \a
40     visitor operates on pair, i.e., \c Visitor (\c mate, \c read) is
41     called.
42
43     Type Requirements:
44     - \c Iterator must be an \input_iterator
45     - \c Iterator 's \c reference type must be convertible to BamRead
46     - \c Visitor must have an \c operator()(BamRead, BamRead) (or any
47       \c const or reference combination)
48
49     \note Input range \c [\a first, \a last \c ) must be sorted or behaviour is
50     undefined.
51
52     \since New in yat 0.10
53   */
54  template<class Iterator, class Visitor>
55  void bam_pair_analyse(Iterator first, Iterator last, Visitor& visitor)
56  {
57    BOOST_CONCEPT_ASSERT((boost::InputIterator<Iterator>));
58
59    using boost::Convertible;
60    typedef typename std::iterator_traits<Iterator>::reference reference_type;
61    BOOST_CONCEPT_ASSERT((Convertible<reference_type, BamRead>));
62
63    // map in which key is mate's position
64    typedef std::multimap<GenomicPosition, BamRead> Map;
65    Map reads;
66    std::map<std::string, BamRead> siam_reads;
67    GenomicPosition siam_position(0,0);
68    for (; first!=last; ++first) {
69      GenomicPosition position(first->core().tid, first->pos());
70      GenomicPosition mate_position(first->core().mtid, first->mpos());
71
72      if (position != siam_position) {
73        siam_position = position;
74        siam_reads.clear();
75      }
76
77      // have not seen the mate yet
78      if (mate_position > position)
79        reads.insert(std::make_pair(mate_position, *first));
80      else if (mate_position < position) {
81        // erase all reads with mate position less than current
82        // position (assuming input range is sorted)
83        Map::iterator lower = reads.lower_bound(position);
84        reads.erase(reads.begin(), lower);
85        // iterator may have been invalidated
86        lower = reads.begin();
87        // find mate and let visitor act on the pair
88        for (; lower!=reads.end() && lower->first == position; ++lower) {
89          if (same_query_name(lower->second, *first)) {
90            visitor(lower->second, *first);
91            break;
92          }
93        }
94      }
95      else { // read and mate have same position
96        // find mate
97        std::map<std::string, BamRead>::iterator
98          mate = siam_reads.lower_bound(first->name());
99        if (mate!=siam_reads.end() && same_query_name(mate->second, *first)) {
100          visitor(mate->second, *first);
101        }
102        else
103          // insert with hint
104          siam_reads.insert(mate, std::make_pair(first->name(), *first));
105      }
106    }
107  }
108
109}}}
110#endif
Note: See TracBrowser for help on using the repository browser.