find_interval.hpp

Go to the documentation of this file.

#ifndef FIND_INTERVAL_HPP_
#define FIND_INTERVAL_HPP_
#include <utility>
#include <deque>
#include "vdj_pipe/sequence_record.hpp"
#include "vdj_pipe/nucleotide_index.hpp"

namespace vdj_pipe{

inline std::pair<std::size_t,std::size_t> longest_average_interval(
         Qual_record::quality const& q,
         const std::size_t min_q,
         const std::size_t win
) {
   std::pair<std::size_t,std::size_t> p(0,0);
   if( q.size() < win ) return p;
   std::size_t cqs = 0; //current sum of scores
   for( std::size_t i = 0; i != win; ++i ) cqs += q[i];

   std::size_t lo = 0, up = win;
   bool expanding = cqs >= min_q;
   for( ; up != q.size(); ++up ) {
      cqs -= q[up - win];
      cqs += q[up];
      if( cqs < min_q ) {
         if( expanding ) {
            if( p.second < up - lo ) {
               p.first = lo;
               p.second = up - lo;
            }
            expanding = false;
         }
      } else {
         if( ! expanding ) {
            lo = up - win + 1;
            expanding = true;
         }
      }
   }
   if( expanding && p.second < up - lo ) {
      p.first = lo;
      p.second = up - lo;
   }
   return p;
}

inline std::pair<std::size_t,std::size_t> longest_min_interval(
         Qual_record::quality const& q,
         const Qual_record::quality::value_type min_q
) {
   std::pair<std::size_t,std::size_t> p(0,0);
   std::size_t lo = 0, up = 0;
   bool expanding = false;
   for( ; up != q.size(); ++up) {
      if( q[up] < min_q ) {
         if( expanding ) {
            if( p.second < up - lo ) {
               p.first = lo;
               p.second = up - lo;
            }
            expanding = false;
         }
      } else {
         if( ! expanding ) {
            lo = up;
            expanding = true;
         }
      }
   }
   if( expanding && p.second < up - lo ) {
      p.first = lo;
      p.second = up - lo;
   }
   return p;
}

inline std::pair<std::size_t,std::size_t> longest_unambiguous_interval(
         Seq_record::sequence const& seq,
         const std::size_t max
) {
   std::pair<std::size_t,std::size_t> p(0,0);
   std::size_t lo = 0, up = 0;

   if( max ) {
      std::deque<std::size_t> amb_pos;
      for( ; up != seq.size(); ++up) {
         if( is_ambiguous(seq[up]) ) {
            if( amb_pos.size() == max ) {
               if( p.second < up - lo ) {
                  p.first = lo;
                  p.second = up - lo;
               }
               lo = amb_pos.front() + 1;
               amb_pos.pop_front();
            }
            amb_pos.push_back(up);
         }
      }
   } else {
      for( ; up != seq.size(); ++up) {
         if( is_ambiguous(seq[up]) ) {
            if( p.second < up - lo ) {
               p.first = lo;
               p.second = up - lo;
            }
            lo = up + 1;
         }
      }
   }

   if( p.second < up - lo ) {
      p.first = lo;
      p.second = up - lo;
   }
   return p;
}

}//namespace vdj_pipe
#endif /* FIND_INTERVAL_HPP_ */