Changeset 3173


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

a first version of BamPairIterator?. refs #784

Location:
trunk
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/Makefile.am

    r3170 r3173  
    4646  test/bam_header.test \
    4747  test/bam_header2.test \
     48  test/bam_pair_iterator.test \
    4849  test/codon.test test/commandline.test \
    4950  test/concept.test \
  • trunk/yat/omic/Makefile.am

    r2994 r3173  
    11## $Id$
    22
    3 # Copyright (C) 2010, 2011, 2012, 2013 Peter Johansson
     3# Copyright (C) 2010, 2011, 2012, 2013, 2014 Peter Johansson
    44#
    55# This file is part of the yat library, http://dev.thep.lu.se/yat
     
    3636nobase_include_HEADERS += $(srcdir)/yat/omic/BamFile.h
    3737nobase_include_HEADERS += $(srcdir)/yat/omic/BamHeader.h
     38nobase_include_HEADERS += $(srcdir)/yat/omic/BamPairIterator.h
    3839nobase_include_HEADERS += $(srcdir)/yat/omic/BamRead.h
    3940nobase_include_HEADERS += $(srcdir)/yat/omic/BamReadFilter.h
  • 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.