source: trunk/yat/omic/BamPairIterator.h @ 3193

#ifndef theplu_yat_omic_bam_pair_iterator
#define theplu_yat_omic_bam_pair_iterator

// $Id: BamPairIterator.h 3193 2014-04-22 05:45:36Z peter $

/*
  Copyright (C) 2014 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library 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.

  The yat library 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 <iostream>

#include "BamPair.h"
#include "BamRead.h"

#include <boost/concept/assert.hpp>
#include <boost/iterator/iterator_facade.hpp>
#include <boost/shared_ptr.hpp>

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

namespace theplu {
namespace yat {
namespace omic {

  /**
     \since New in yat 0.12

     Type Requirements:
     - \c Iterator must be an \input_iterator
     - \c Iterator 's \c reference type must be convertible to BamRead

     This iterator works on a sorted range [begin, end) and provies a
     covenient way to access the pairs rather than the reads
     individually. When BamPairIterator is incremented, rather than
     just step to next element in [begin, end), BamPairIterator keeps
     stepping until it finds a read whose mate it has already seen, and
     operator* return the pair of BamReads.

     Note that BamPairIterator is a single-pass iterator, i.e., once it
     is incremented the behaviour of its copies is undefined.

     \note \c operator-> is not working with Boost 1.48 (or
     older). See <a
     href="https://svn.boost.org/trac/boost/ticket/5697">Boost Bug
     #5697</a>
   */
  template<typename Base>
  class BamPairIterator
    : public boost::iterator_facade<
    BamPairIterator<Base>, const BamPair, std::input_iterator_tag,
    const BamPairProxy
    >
  {
  public:
    /**
       \brief default constructor
     */
    BamPairIterator(void);

    /**
       Creates an iterator that will work on [begin, end).

       \see bam_pair_iterator

       \note Input range \c [\a first, \a last \c ) must be sorted or
       behaviour is undefined.
     */
    explicit BamPairIterator(Base begin, Base end);
  private:
    Base iter_;
    Base end_;
    const BamRead* mate_;
    boost::shared_ptr<std::map<std::string, BamRead> > siam_reads_;
    typedef std::pair<int32_t, int32_t> Position;
    boost::shared_ptr<std::multimap<Position, BamRead> > reads_;
    friend class boost::iterator_core_access;

    const BamPairProxy dereference(void) const;
    bool equal(const BamPairIterator& other) const;
    void increment(void);
    void find_next(void);
  };


  /**
     Convenient function to create BamPairIterator

     \relates BamPairIterator
   */
  template<typename Base>
  BamPairIterator<Base> bam_pair_iterator(Base base, Base end)
  { return BamPairIterator<Base>(base, end); }


  /**
     Convenient function to create BamPairIterator that works on empty
     range. Equivalent to bam_pair_iterator(end, end).

     \relates BamPairIterator
  */
  template<typename Base>
  BamPairIterator<Base> bam_pair_iterator(Base end)
  { return BamPairIterator<Base>(end, end); }


  //             template implementations
  ///////////////////////////////////////////////////////////////

  template<typename Base>
  BamPairIterator<Base>::BamPairIterator(Base base, Base end)
    : iter_(base), end_(end),
      siam_reads_(new std::map<std::string, BamRead>),
      reads_(new std::multimap<Position, BamRead>)
  {
    BOOST_CONCEPT_ASSERT((boost::InputIterator<Base>));
    using boost::Convertible;
    typedef typename std::iterator_traits<Base>::reference reference_type;
    BOOST_CONCEPT_ASSERT((Convertible<reference_type, BamRead>));
    find_next();
  }


  template<typename Base>
  BamPairIterator<Base>::BamPairIterator(void)
  {
    BOOST_CONCEPT_ASSERT((boost::InputIterator<Base>));
    using boost::Convertible;
    typedef typename std::iterator_traits<Base>::reference reference_type;
    BOOST_CONCEPT_ASSERT((Convertible<reference_type, BamRead>));
  }


  template<typename Base>
  const BamPairProxy
  BamPairIterator<Base>::dereference(void) const
  {
    return BamPairProxy(mate_, &*iter_);
  }


  template<typename Base>
  bool BamPairIterator<Base>::equal(const BamPairIterator& other) const
  {
    return iter_ == other.iter_;
  }


  template<typename Base>
  void BamPairIterator<Base>::increment(void)
  {
    ++iter_;
    find_next();
  }


  template<typename Base>
  void BamPairIterator<Base>::find_next(void)
  {
    BamLessPos less;
    for (; iter_!=end_; ++iter_) {
      Position position(iter_->tid(), iter_->pos());
      Position mate_position(iter_->mtid(), iter_->mpos());
      // clear siam reads if iter is more advanced than siam reads
      if (siam_reads_->size() && less(siam_reads_->begin()->second, *iter_))
        siam_reads_->clear();

      // have not seen mate yet
      if (mate_position > position)
        reads_->insert(std::make_pair(mate_position, *iter_));
      else if (position > mate_position) {
        std::multimap<Position, BamRead>::iterator
          lower = reads_->lower_bound(position);
        // erase all reads with mate position less than current
        // position (assuming input range is sorted)
        reads_->erase(reads_->begin(), lower);

        // find mate and assign pair
        for (; lower!=reads_->end() && lower->first == position; ++lower)
          if (same_query_name(lower->second, *iter_)) {
            mate_ = &lower->second;
            return;
          }
      }
      else { // siam read, i.e., iter_ and mate have same position
        // check if we have already seen mate and stored it in map
        std::map<std::string, BamRead>::iterator
          mate = siam_reads_->lower_bound(iter_->name());
        if (mate!=siam_reads_->end() && same_query_name(mate->second, *iter_)) {
          mate_ = &mate->second;
          return;
        }
        else
          // insert with hint
          siam_reads_->insert(mate, std::make_pair(iter_->name(), *iter_));
      }
    }
  }

}}}
#endif