sequence_fls.hpp

Go to the documentation of this file.

#ifndef SEQUENCE_FLS_HPP_
#define SEQUENCE_FLS_HPP_
#include <algorithm>
#include <ostream>
#include "boost/assert.hpp"
#include "boost/cstdint.hpp"
#include "boost/integer/static_log2.hpp"
#include "vdj_pipe/nucleotide_index.hpp"
#include "vdj_pipe/exception.hpp"
#include "vdj_pipe/detail/string_ref.hpp"

namespace vdj_pipe{

template<unsigned N = 4, typename S = boost::uint_least32_t> class Seq_fls {
   static const unsigned BITS_PER_CHAR = boost::static_log2<N>::value;
   static const unsigned MAX_CHARS = sizeof(S) * 8 / BITS_PER_CHAR;
   static const unsigned MASK = ~(~0 << BITS_PER_CHAR);

public:
   static const unsigned alphabet_size = N;
   static std::size_t length() {return MAX_CHARS;}
   typedef Seq_fls self_type;

   Seq_fls() {}

   explicit Seq_fls(const boost::string_ref seq) : val_(0) {
      assign(seq);
   }

   Seq_fls(self_type const& ss) : val_(ss.val_) {}
   std::size_t size() const {return MAX_CHARS;}
   self_type& operator=(self_type const& ss) {
      val_ = ss.val_;
      return *this;
   }

   self_type& operator=(const boost::string_ref seq) {
      val_ = 0;
      assign(seq);
      return *this;
   }

   bool operator==(self_type const& sf) const {return val_ == sf.val_;}
   bool operator!=(self_type const& sf) const {return !(*this == sf);}

   Nucleotide operator[](const std::size_t i) const {
      BOOST_ASSERT(i < MAX_CHARS);
      return Nucleotide((val_>> i * BITS_PER_CHAR) & MASK);
   }

   void set(const unsigned i, const Nucleotide n) {
      BOOST_ASSERT(i < MAX_CHARS);
      BOOST_ASSERT(n < alphabet_size);
      val_ &= ~(MASK << i * BITS_PER_CHAR);
      val_ |= n << i * BITS_PER_CHAR;
   }

   std::string to_string(std::size_t n = self_type::length()) const {
      std::string s(n, 'X');
      n = std::min(self_type::length(), n);
      for(std::size_t i = 0; i != n; ++i) {
         s[i] = to_capital(this->operator[](i));
      }
      return s;
   }

   std::ostream& print(
            std::ostream& os,
            std::size_t n = self_type::length()
   ) const {
      n = std::min(self_type::length(), n);
      for(std::size_t i = 0; i != n; ++i) {
         os << to_capital(this->operator[](i));
      }
      return os;
   }

private:
   S val_;

   void assign(const boost::string_ref seq) {
      if( seq.size() > size() ) BOOST_THROW_EXCEPTION(
         base_exception()
         << base_exception::msg_t("sequence too long")
         << base_exception::int1_t(seq.size())
         << base_exception::int2_t(size())
      );
      for(unsigned i = seq.size(); i != 0; --i) {
         const Nucleotide n = nucleotide_index(seq[i-1]);
         if( n >= alphabet_size ) {
            std::string c("\"");
            c += seq[i-1];
            c += '\"';
            BOOST_THROW_EXCEPTION(
               base_exception()
               << base_exception::msg_t("character not supported")
               << base_exception::str1_t(sanitize(seq[i-1]))
            );
         }
         val_ <<= BITS_PER_CHAR;
         val_ |= n;
      }
   }

   friend inline std::size_t hash_value(self_type const& seq) {return seq.val_;}

   friend std::ostream& operator<<(std::ostream& os, self_type const& seq) {
      return seq.print(os);
   }

   friend bool assign(self_type& seq1, const boost::string_ref seq2) {
      if( seq2.size() > length() ) return false;
      seq1.val_ = 0;
      for(unsigned i = seq2.size(); i != 0; --i) {
         const Nucleotide n = nucleotide_index(seq2[i-1]);
         if( n >= alphabet_size ) return false;
         seq1.val_ <<= BITS_PER_CHAR;
         seq1.val_ |= n;
      }
      return true;
   }
};

template<class Seq> inline int
identity(
         Seq const& s1,
         const boost::string_ref s2,
         scoring_matrix_t const& sm,
         const std::size_t = 0
) {
   int i = 0;
   for(std::size_t n = 0; n != s1.size() && n != s2.size(); ++n) {
      const Nucleotide n1 = s1[n];
      const Nucleotide n2 = nucleotide_index(s2[n]);
      i += identity(n1, n2, sm);
   }
   return i;
}

template<unsigned N, typename S> inline int
identity(
         const Seq_fls<N,S> seq1,
         const Seq_fls<N,S> seq2,
         scoring_matrix_t const& sm,
         const unsigned len = N
) {
   int i = 0;
   for( unsigned j = 0; j != len; ++j ) {
      i += identity(seq1[j], seq2[j], sm);
   }
   return i;
}

}//namespace vdj_pipe
#endif /* SEQUENCE_FLS_HPP_ */