sequence_map_aligned.hpp

Go to the documentation of this file.

#ifndef SEQUENCE_MAP_ALIGNED_HPP_
#define SEQUENCE_MAP_ALIGNED_HPP_
#include <limits>
#include <stdint.h>
#include <vector>

#include "boost/assert.hpp"
#include "boost/foreach.hpp"

#include "external/ssw_cpp.h"
#include "vdj_pipe/detail/id_bimap.hpp"
#include "vdj_pipe/detail/id_iterator.hpp"
#include "vdj_pipe/exception.hpp"
#include "vdj_pipe/gdst/gdst.hpp"
#include "vdj_pipe/interval_iterator.hpp"
#include "vdj_pipe/parser_fasta.hpp"
#include "vdj_pipe/sequence_interval.hpp"
#include "vdj_pipe/sequence_transform.hpp"

namespace vdj_pipe{

class Seq_map_aligned {
   typedef gdst::Gdst suffix_tree;
   typedef suffix_tree::seq_map seq_map;
   typedef detail::Id_bimap<Seq_id, std::string> id_map;
public:
   typedef Id_iterator<Seq_id> iterator;
   typedef iterator const_iterator;
   struct Err : public base_exception{};

   struct Match {
      Match()
      : id_(),
        score_(std::numeric_limits<int>::min()),
        n_mismatches_(std::numeric_limits<unsigned>::max()),
        identity_(0.0),
        si_(sequence_interval_invalid()),
        cigar_()
      {}

      Seq_id id_;
      int score_;
      unsigned n_mismatches_;
      double identity_;
      sequence_interval si_;
      std::vector<uint32_t> cigar_;
   };

   Seq_map_aligned()
   : st_(),
     idm_(Seq_id(1))
   {}

   iterator begin() const {return iterator(st_.sequence_map().min_id());}
   iterator end() const {return ++iterator(st_.sequence_map().max_id());}
   Seq_id max_id() const {return st_.sequence_map().max_id();}
   std::size_t size() const {return st_.sequence_map().size();}

   void load_file(std::string const& fasta, const bool reverse) {
      for(Parser_fasta pf((File_input(fasta))); pf.has_next(); pf.next_record() ) {
         const Seq_record sr = pf.get_record();
         insert(sr.id_, sr.seq_, reverse);
      }
   }

   void insert(std::string const& id, std::string seq, const bool reverse) {
      check(seq);
      std::pair<Seq_id,bool> p = idm_.insert(id);
      if( ! p.second ) BOOST_THROW_EXCEPTION(
               Err()
               << Err::msg_t("duplicate sequence")
               << Err::str1_t(sanitize(id))
               << Err::str2_t(sanitize(seq))
      );
      if( reverse ) {
         st_.sequence_map().insert(p.first) = Seq_entry(complement(seq));
      } else {
         st_.sequence_map().insert(p.first) = Seq_entry(seq);
      }

      st_.insert(p.first);
   }

   Seq_entry const& seq(const Seq_id sid) const {
      BOOST_ASSERT(st_.sequence_map().find(sid));
      return st_.sequence_map()[sid];
   }

   std::string const& seq_id(const Seq_id sid) const {
      BOOST_ASSERT(st_.sequence_map().find(sid));
      return idm_[sid];
   }

   Match best_match(
            const boost::string_ref seq,
            const std::size_t min_match,
            const bool track_mismatches
   ) const {
      detail::Vector_set<Seq_id> vs = find_overlaping(seq, min_match);
      Match m;
      BOOST_FOREACH(const Seq_id sid, vs) {
         StripedSmithWaterman::Aligner aligner;
         StripedSmithWaterman::Filter filter;
         StripedSmithWaterman::Alignment al;
         Seq_entry const& se = st_.sequence_map()[sid];
         aligner.Align(seq.data(), se.sequence().data(), se.size(), filter, &al);
         if(
                  (
                           track_mismatches &&
                           ( m.n_mismatches_ > (unsigned)al.mismatches )
                  ) ||
                  ( ! track_mismatches && ( m.score_ < al.sw_score ) )
         ) {
            m.id_ = sid;
            m.score_ = al.sw_score;
            m.n_mismatches_ = (unsigned)al.mismatches;
            m.cigar_ = al.cigar;
            m.si_ = sequence_interval(al.query_begin, al.query_end);
            double len =
                     0.5 *
                     (al.ref_end - al.ref_begin + al.query_end - al.query_begin);
            m.identity_ = (len - al.mismatches) / len;
         }
      }
      return m;
   }

   detail::Vector_set<Seq_id> find_overlaping(
            const boost::string_ref seq,
            const std::size_t min_match
   ) const {
      detail::Vector_set<Seq_id> vs;
      for(
               Unambiguous_interval_iter uii(seq, min_match);
               uii.has_subseq();
               uii.next()
      ) {
         st_.find_overlaping(uii.subseq(), vs, min_match);
      }
      return vs;
   }

private:
   suffix_tree st_;
   id_map idm_;

   static void check(const boost::string_ref seq) {
      for(std::size_t n = 0; n != seq.size(); ++n) {
         if( is_ambiguous(seq[n]) ) BOOST_THROW_EXCEPTION(
                  Err()
                  << Err::msg_t("undefined nucleotide")
                  << Err::str1_t(sanitize(seq, 500))
                  << Err::str2_t(sanitize(seq[n]))
                  << Err::int1_t(n)
         );
      }
   }
};


}//namespace vdj_pipe
#endif /* SEQUENCE_MAP_ALIGNED_HPP_ */