gdst.hpp

Go to the documentation of this file.

#ifndef GDST_HPP_
#define GDST_HPP_
#include "boost/assert.hpp"
#include "boost/foreach.hpp"
#include "boost/shared_ptr.hpp"

#include "vdj_pipe/config.hpp"
#include "vdj_pipe/detail/disable_warnings.hpp"
#include "vdj_pipe/detail/id_map.hpp"
#include "vdj_pipe/detail/string_ref.hpp"
#include "vdj_pipe/exception.hpp"
#include "vdj_pipe/gdst/common_subsequence.hpp"
#include "vdj_pipe/gdst/depth_first_iterator.hpp"
#include "vdj_pipe/gdst/nodes.hpp"
#include "vdj_pipe/gdst/pattern_iterator.hpp"
#include "vdj_pipe/nucleotide_index.hpp"
#include "vdj_pipe/object_ids.hpp"
#include "vdj_pipe/sequence_entry.hpp"

namespace vdj_pipe{ namespace gdst{
class Gdst_stats;
class Common_substrings;
class Ukkonen_inserter;

struct VDJ_PIPE_DECL Match {
   explicit Match(
            const Branch_id bid = Branch_id(),
            const unsigned min_d = 0,
            const unsigned i = 0
   )
   :curr_(bid),
    min_d_(min_d),
    i_(i),
    at_pattern_end_(),
    pn_(),
    at_tree_end_(),
    next_(),
    mismatch_(),
    en_()
   {}

   void next(const boost::string_ref patt, Gdst const& st);
   bool check(const boost::string_ref patt, Gdst const& st);

   Branch_id curr_; 
   unsigned min_d_; 
   unsigned i_; 
   bool at_pattern_end_; 
   Nucleotide pn_; 
   bool at_tree_end_; 
   Branch_id next_; 
   bool mismatch_; 
   Nucleotide en_; 
};

class VDJ_PIPE_DECL Gdst {
   typedef detail::Id_map<Branch_id, Branch> branch_map;
   typedef detail::Id_map<Leaf_id, Leaf> leaf_map;
   typedef detail::Id_map<Children_id, Children> children_map;
   typedef boost::string_ref seq_type;
   friend class Gdst_stats;
   friend class Common_substrings;
   friend class Ukkonen_inserter;

public:
   struct Err : public base_exception{};

   typedef detail::Id_map<Seq_id, Seq_entry> seq_map;
   typedef boost::shared_ptr<seq_map> seq_map_ptr;
   typedef Seq_pos<Seq_id> match_type;
   typedef std::vector<match_type> match_vector;

   Gdst()
   : ss_(new seq_map(Seq_id(1))),
     bm_(Branch_id(1)),
     cm_(Children_id(1)),
     lm_(Leaf_id(1)),
     root_(bm_.insert(Branch(0)))
   {}

   explicit Gdst(seq_map_ptr ss)
   : ss_(ss),
     bm_(Branch_id(1)),
     cm_(Children_id(1)),
     lm_(Leaf_id(1)),
     root_(bm_.insert(Branch(0)))
   {}

   std::size_t size() const {return bm_.size();}
   Depth_iter depth_first() const {return Depth_iter(*this, root_);}
   Branch_id root() const {return root_;}
   seq_map const& sequence_map() const {return *ss_;}
   seq_map& sequence_map() {return *ss_;}

   Branch const& operator[](const Branch_id bid) const {
      BOOST_ASSERT(bid);
      return bm_[bid];
   }

   Leaf const& operator[](const Leaf_id lid) const {
      BOOST_ASSERT(lid);
      return lm_[lid];
   }

   Seq_entry const& operator[](const Seq_id sid) const {
      return (*ss_)[sid];
   }

   Leaf const& leaf(const Branch_id bid) const {
      BOOST_ASSERT(bid);
      static const Leaf empty = Leaf();
      if( const Leaf_id lid = bm_[bid].leaf_ ) return lm_[lid];
      return empty;
   }

   const seq_type suffix(const Seq_id sid, const unsigned len) const {
      BOOST_ASSERT(sid);
      BOOST_ASSERT((*ss_)[sid].size());
      BOOST_ASSERT((*ss_)[sid].size() >= len);
      const seq_type seq = (*ss_)[sid].sequence();
      return seq.substr(seq.length() - len);
   }

   const seq_type suffix(const Leaf_id lid, const unsigned len) const {
      BOOST_ASSERT(lid);
      BOOST_ASSERT(lm_[lid].size());
      return suffix(lm_[lid].front(), len);
   }

   const seq_type suffix(Branch const& b) const {
      if( b.leaf_ ) return suffix(b.leaf_, b.n_);
      BOOST_ASSERT(b.sb_);
      Branch const& suffix_branch = bm_[b.sb_];
      BOOST_ASSERT(suffix_branch.leaf_);
      return suffix(suffix_branch.leaf_, suffix_branch.n_);
   }

   const seq_type suffix(const Branch_id bid) const {
      BOOST_ASSERT(bid);
      return suffix(bm_[bid]);
   }

   Nucleotide letter(
            const Branch_id bid,
            const unsigned i
   ) const {
      BOOST_ASSERT(i < bm_[bid].n_);
      const seq_type seq = suffix(bid);
      return nucleotide_index(seq[i]);
   }

   Nucleotide letter(
            const Branch_id bid,
            const Nucleotide n,
            const unsigned i
   ) const {
      BOOST_ASSERT(bid);
      BOOST_ASSERT(bm_[bid].n_ < i);
      const Branch_id c = child(bid, n);
      BOOST_ASSERT(c);
      BOOST_ASSERT(i < bm_[c].n_);
      const seq_type seq = suffix(c);
      return nucleotide_index(seq[i]);
   }

   Branch_id child(const Branch_id bid1, const Nucleotide n) const {
      BOOST_ASSERT(bid1);
      if( n >= 4 ) BOOST_THROW_EXCEPTION(
               Err()
               << Err::msg_t("invalid nucleotide (only ACGT supported)")
               << Err::str1_t(sanitize(to_capital(n)))
      );
      Branch const& b1 = bm_[bid1];
      return b1.c_ ? cm_[b1.c_][n] : Branch_id();
   }

   Branch_id child(const Branch_id bid1, const seq_type seq) const {
      BOOST_ASSERT(bid1);
      Branch const& b1 = bm_[bid1];
      if( ! b1.c_ ) return Branch_id();
      BOOST_ASSERT(seq.length() > b1.n_);
      const Nucleotide n = nucleotide_index(seq[b1.n_]);
      return cm_[b1.c_][n];
   }

   void child(
            const Branch_id bid1,
            const Nucleotide n,
            const Branch_id bid2
   ) {
      BOOST_ASSERT(bid1);
      BOOST_ASSERT(bid2);
      BOOST_ASSERT(bm_[bid1].n_ < bm_[bid2].n_);
      Branch& b1 = bm_[bid1];
      if( ! b1.c_ ) b1.c_ = cm_.insert(Children());
      cm_[b1.c_][n] = bid2;
   }

   Common_subseq find_longest(
            const seq_type seq,
            std::size_t min_len = 0
   ) const;

   void find_overlaping(
            const boost::string_ref seq,
            detail::Vector_set<Seq_id>& vs,
            std::size_t min_overlap = 0
   ) const;

   Match find(
            const seq_type seq,
            const Branch_id bid,
            unsigned min_d = 0
   ) const {
      Match m(bid, min_d);
      for( ; m.check(seq, *this); m.next(seq, *this) );
      return m;
   }

   void insert(const Seq_id sid);

private:
   seq_map_ptr ss_;
   branch_map bm_;
   children_map cm_;
   leaf_map lm_;
   Branch_id root_;

   void collect_sequences(
            const gdst::Branch_id nid,
            detail::Vector_set<Seq_id>& vs
   ) const {
      for( gdst::Depth_iter di(*this, nid); ! di.at_end(); di.next() ) {
         gdst::Leaf const& l = leaf(di.id());
         vs.insert(l.begin(), l.end());
      }
   }

   void suffix_link(const Branch_id bid1, const Branch_id bid2) {
      BOOST_ASSERT(bid1);
      BOOST_ASSERT(bid2);
      BOOST_ASSERT( ! bm_[bid1].sl_);
      BOOST_ASSERT( bm_[bid1].n_ == bm_[bid2].n_ + 1);
      bm_[bid1].sl_ = bid2;
   }

   Branch_id suffix_link(const Branch_id bid) const {
      BOOST_ASSERT(bid);
      return bm_[bid].sl_;
   }

   unsigned edge_length(const Branch_id bid, const Nucleotide n) const {
      const Branch_id c = child(bid, n);
      BOOST_ASSERT(c);
      BOOST_ASSERT(bm_[bid].n_ < bm_[c].n_);
      return bm_[c].n_ - bm_[bid].n_;
   }

   unsigned edge_length(const Branch_id bid, const seq_type suff) const {
      BOOST_ASSERT(bid);
      BOOST_ASSERT(bm_[bid].n_ < suff.size());
      return edge_length(bid, nucleotide_index(suff[bm_[bid].n_]));
   }

   void check_edge(
            Branch_id& an,
            unsigned& aei,
            unsigned& al,
            const boost::string_ref seq
   ) const {
      BOOST_ASSERT(an);
      for( ; ; ) {
         BOOST_ASSERT(aei < seq.size());
         const Nucleotide ae = nucleotide_index(seq[aei]);
         const Branch_id c = child(an, ae);
         BOOST_ASSERT(c);
         BOOST_ASSERT(bm_[an].n_ < bm_[c].n_);
         const unsigned l = bm_[c].n_ - bm_[an].n_;
         if( al < l ) return;
         an = c;
         aei += l;
         al -= l;
      }
   }

   void add_to_leaf(const Branch_id bid, const Seq_id sid) {
      Branch& b = bm_[bid];
      if( ! b.leaf_ ) {
         b.leaf_ = lm_.insert(Leaf());
      }
      lm_[b.leaf_].insert(sid);
   }

   Branch_id leaf_from_branch(
            const Branch_id bid1,
            const Nucleotide n1,
            const unsigned i,
            const Seq_id sid
   ) {
      const Branch_id bid2 = bm_.insert(Branch(i));
      add_to_leaf(bid2, sid);
      bm_[bid2].sb_ = bid2;  //its own suffix branch
      child(bid1, n1, bid2);
      return bid2;
   }

   Branch_id split_edge(
            const Branch_id bid1,
            const Nucleotide n1,
            const unsigned i,
            const Nucleotide n2
   ) {
      BOOST_ASSERT(bm_[bid1].n_ < i && "new child index should be greater");
      const Branch_id bid2 = bm_.insert(Branch(i));
      const Branch_id bid3 = child(bid1, n1);
      BOOST_ASSERT(i < bm_[bid3].n_ && "old child index should be greater");
      child(bid1, n1, bid2);
      child(bid2, n2, bid3);
      return bid2;
   }
};

}//namespace gdst
}//namespace vdj_pipe
#endif /* GDST_HPP_ */