me_ambiguous.hpp

Go to the documentation of this file.

#ifndef ME_AMBIGUOUS_HPP_
#define ME_AMBIGUOUS_HPP_
#include <limits>
#include "boost/foreach.hpp"
#include "boost/utility/string_ref.hpp"
#include "me_relative_interval.hpp"
#include "sequence_map_full.hpp"
#include "sequence_map_types.hpp"
#include "me_types.hpp"
#include "vdj_pipe/exception.hpp"
#include "vdj_pipe/object_ids.hpp"
#include "vdj_pipe/sequence_record.hpp"
#include "vdj_pipe/value_map.hpp"

namespace vdj_pipe{ namespace match{

class Match_element_ambiguous {
public:
   struct Err : public base_exception{};
   typedef boost::string_ref sequence;
   typedef Qual_record::quality quality;

   Match_element_ambiguous(
            Value_map const& vm,
            std::string const& match_value_name,
            std::string const& score_value_name,
            Relative_interval const& ri,
            detail::name_seq_vector const& sv,
            const int min_score,
            const bool ignore_dups,
            const bool require_best,
            const bool track_mismatches
   )
   : vm_(vm),
     name_val_id_(vm_.insert_new_name(match_value_name)),
     score_val_id_(vm_.insert_new_name(score_value_name)),
     sm_(ignore_dups),
     ri_(ri),
     scoring_matrix_(
              track_mismatches ?
                       &scoring_matrix<0,-1,0,0>() :
                       &scoring_matrix<2,-2,1,0>()
     ),
     min_score_(track_mismatches ? - min_score : min_score),
     require_best_(require_best)
   {
      if( sv.empty() ) BOOST_THROW_EXCEPTION(
               Err()
               << Err::msg_t("no sequences provided")
      );

      if( ! ri_.length_defined() ) {
         ri_ = Relative_interval(ri_.position(), sv.front().second.size());
      }

      BOOST_FOREACH(detail::match_seq const& ms, sv) {
         if( ms.first.empty() ) {
            sm_.insert(ms.second, ms.second);
         } else {
            sm_.insert(ms.first, ms.second);
         }
      }
   }

   sequence_interval operator() (
            const sequence seq,
            quality const& qual,
            sequence_interval const& si
   ) {
      sequence_interval res = sequence_interval_invalid();
      if( ! is_valid(si) ) return res;
      const sequence_interval si1 = ri_(si, seq.size());
      if( ! is_valid(si1) || (std::size_t)width(si1) < sm_.min_len() ) {
         return res;
      }

      std::size_t best_pos = 0;
      Seq_map_full::match_type m;

      for(
            std::size_t n = si1.lower(), max = si1.upper() - sm_.min_len() + 1;
            n != max;
            ++n
      ) {
         const sequence s = seq.substr(n);
         if( sm_.find_closest(s, m, *scoring_matrix_) ) {
            best_pos = n;
         }
      }

      if( m.has_score() && score_val_id_ ) {
         vm_[score_val_id_] = (long)m.score1();
      }

      if( m.is_acceptable(min_score_, require_best_) ) {
         if( name_val_id_ ) vm_[name_val_id_] = sm_.name(m.id1());
         res.assign(best_pos, best_pos + sm_.seq(m.id1()).size());
      }

      return res;
   }

   void finish() {}

private:
   Value_map vm_;
   Val_id name_val_id_;
   Val_id score_val_id_;
   Seq_map_full sm_;
   Relative_interval ri_;
   scoring_matrix_t const* scoring_matrix_;
   int min_score_;
   bool require_best_;
};

}//namespace match
}//namespace vdj_pipe
#endif /* ME_AMBIGUOUS_HPP_ */