algorithm.h

#ifndef theplu_yat_omic_algorithm
#define theplu_yat_omic_algorithm

// $Id: algorithm.h 2905 2012-12-14 00:03:10Z peter $
//
// Copyright (C) 2012 Peter Johansson
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.

#include "BamRead.h"
#include "GenomicPosition.h"

#include <boost/concept/assert.hpp>

#include <iterator>
#include <map>
#include <string>
#include <utility>

namespace theplu {
namespace yat {
namespace omic {

  template<class Iterator, class Visitor>
  void bam_pair_analyse(Iterator first, Iterator last, Visitor& visitor)
  {
    BOOST_CONCEPT_ASSERT((boost::InputIterator<Iterator>));

    using boost::Convertible;
    typedef typename std::iterator_traits<Iterator>::reference reference_type;
    BOOST_CONCEPT_ASSERT((Convertible<reference_type, BamRead>));

    // map in which key is mate's position
    typedef std::multimap<GenomicPosition, BamRead> Map;
    Map reads;
    std::map<std::string, BamRead> siam_reads;
    GenomicPosition siam_position(0,0);
    for (; first!=last; ++first) {
      GenomicPosition position(first->core().tid, first->pos());
      GenomicPosition mate_position(first->core().mtid, first->mpos());

      if (position != siam_position) {
        siam_position = position;
        siam_reads.clear();
      }

      // have not seen the mate yet
      if (mate_position > position)
        reads.insert(std::make_pair(mate_position, *first));
      else if (mate_position < position) {
        // erase all reads with mate position less than current
        // position (assuming input range is sorted)
        Map::iterator lower = reads.lower_bound(position);
        reads.erase(reads.begin(), lower);
        // iterator may have been invalidated
        lower = reads.begin();
        // find mate and let visitor act on the pair
        for (; lower!=reads.end() && lower->first == position; ++lower) {
          if (same_query_name(lower->second, *first)) {
            visitor(lower->second, *first);
            break;
          }
        }
      }
      else { // read and mate have same position
        // find mate
        std::map<std::string, BamRead>::iterator
          mate = siam_reads.lower_bound(first->name());
        if (mate!=siam_reads.end() && same_query_name(mate->second, *first)) {
          visitor(mate->second, *first);
        }
        else
          // insert with hint
          siam_reads.insert(mate, std::make_pair(first->name(), *first));
      }
    }
  }

}}}
#endif