Ignore:
Timestamp:
Mar 8, 2014, 2:48:13 PM (8 years ago)
Author:
Peter
Message:

a first version of BamPairIterator?. refs #784

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/yat/omic/algorithm.h

    r3114 r3173  
    2323*/
    2424
     25#include "BamPairIterator.h"
    2526#include "BamRead.h"
    2627#include "GenomicPosition.h"
     
    5960  void bam_pair_analyse(Iterator first, Iterator last, Visitor& visitor)
    6061  {
    61     BOOST_CONCEPT_ASSERT((boost::InputIterator<Iterator>));
    62 
    63     using boost::Convertible;
    64     typedef typename std::iterator_traits<Iterator>::reference reference_type;
    65     BOOST_CONCEPT_ASSERT((Convertible<reference_type, BamRead>));
    66 
    67     // map in which key is mate's position
    68     typedef std::multimap<GenomicPosition, BamRead> Map;
    69     Map reads;
    70     std::map<std::string, BamRead> siam_reads;
    71     GenomicPosition siam_position(0,0);
    72     for (; first!=last; ++first) {
    73       GenomicPosition position(first->core().tid, first->pos());
    74       GenomicPosition mate_position(first->core().mtid, first->mpos());
    75 
    76       if (position != siam_position) {
    77         siam_position = position;
    78         siam_reads.clear();
    79       }
    80 
    81       // have not seen the mate yet
    82       if (mate_position > position)
    83         reads.insert(std::make_pair(mate_position, *first));
    84       else if (mate_position < position) {
    85         // erase all reads with mate position less than current
    86         // position (assuming input range is sorted)
    87         Map::iterator lower = reads.lower_bound(position);
    88         reads.erase(reads.begin(), lower);
    89         // iterator may have been invalidated
    90         lower = reads.begin();
    91         // find mate and let visitor act on the pair
    92         for (; lower!=reads.end() && lower->first == position; ++lower) {
    93           if (same_query_name(lower->second, *first)) {
    94             visitor(lower->second, *first);
    95             break;
    96           }
    97         }
    98       }
    99       else { // read and mate have same position
    100         // find mate
    101         std::map<std::string, BamRead>::iterator
    102           mate = siam_reads.lower_bound(first->name());
    103         if (mate!=siam_reads.end() && same_query_name(mate->second, *first)) {
    104           visitor(mate->second, *first);
    105         }
    106         else
    107           // insert with hint
    108           siam_reads.insert(mate, std::make_pair(first->name(), *first));
    109       }
    110     }
     62    BamPairIterator<Iterator> iter(first, last);
     63    BamPairIterator<Iterator> end(last, last);
     64    for (; iter!=end; ++iter)
     65      visitor(iter->first, iter->second);
    11166  }
    11267
Note: See TracChangeset for help on using the changeset viewer.