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

Last change on this file since 2904 was 2904, checked in by Peter, 10 years ago

boost concept check. refs #729

  • 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 2904 2012-12-14 00:00:22Z 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    using yat::omic::GenomicPosition;
64    // map in which key is mate's position
65    typedef std::multimap<GenomicPosition, BamRead> Map;
66    Map reads;
67    std::map<std::string, BamRead> siam_reads;
68    GenomicPosition siam_position(0,0);
69    for (; first!=last; ++first) {
70      GenomicPosition position(first->core().tid, first->pos());
71      GenomicPosition mate_position(first->core().mtid, first->mpos());
72
73      if (position != siam_position) {
74        siam_position = position;
75        siam_reads.clear();
76      }
77
78      // have not seen the mate yet
79      if (mate_position > position)
80        reads.insert(std::make_pair(mate_position, *first));
81      else if (mate_position < position) {
82        // erase all reads with mate position less than current
83        // position (assuming input range is sorted)
84        Map::iterator lower = reads.lower_bound(position);
85        reads.erase(reads.begin(), lower);
86        // iterator may have been invalidated
87        lower = reads.begin();
88        // find mate and let visitor act on the pair
89        for (; lower!=reads.end() && lower->first == position; ++lower) {
90          if (same_query_name(lower->second, *first)) {
91            visitor(lower->second, *first);
92            break;
93          }
94        }
95      }
96      else { // read and mate have same position
97        // find mate
98        std::map<std::string, BamRead>::iterator
99          mate = siam_reads.lower_bound(first->name());
100        if (mate!=siam_reads.end() && same_query_name(mate->second, *first)) {
101          visitor(mate->second, *first);
102        }
103        else
104          // insert with hint
105          siam_reads.insert(mate, std::make_pair(first->name(), *first));
106      }
107    }
108  }
109
110}}}
111#endif
Note: See TracBrowser for help on using the repository browser.