nucleotide_index.hpp

Go to the documentation of this file.

#ifndef NUCLEOTIDE_INDEX_HPP_
#define NUCLEOTIDE_INDEX_HPP_
#include "boost/algorithm/string/predicate_facade.hpp"
#include "boost/assert.hpp"
#if defined(NDEBUG)
//fix boost::multi_array warnings
//boost/multi_array/base.hpp:503:13: warning: unused variable 'bound_adjustment'
#define BOOST_DISABLE_ASSERTS
#endif
#include "boost/multi_array.hpp"
#include "vdj_pipe/detail/disable_warnings.hpp"
#include "vdj_pipe/exception.hpp"

#define VDJ_PIPE_N_NUCLEOTIDES 16

namespace vdj_pipe{

enum Nucleotide {
   Adenine     = 0,
   Cytosine    = 1,
   Guanine     = 2,
   Thymine     = 3,
   Any         = 4,
   Uracil      = 5,  
   Purine      = 6,  
   Pyrimidine  = 7,  
   Ketone      = 8,  
   Amine       = 9,  
   Strong      = 10, 
   Weak        = 11, 
   not_A       = 12, 
   not_C       = 13, 
   not_G       = 14, 
   not_T       = 15  
};

inline Nucleotide nucleotide_index(const char c) {
   switch (c) {
   case 'A':
   case 'a':
      return Adenine;
   case 'C':
   case 'c':
      return Cytosine;
   case 'G':
   case 'g':
      return Guanine;
   case 'T':
   case 't':
      return Thymine;
   case 'U':
   case 'u':
      return Uracil;
   case 'R':
   case 'r':
      return Purine;
   case 'Y':
   case 'y':
      return Pyrimidine;
   case 'K':
   case 'k':
      return Ketone;
   case 'M':
   case 'm':
      return Amine;
   case 'S':
   case 's':
      return Strong;
   case 'W':
   case 'w':
      return Weak;
   case 'B':
   case 'b':
      return not_A;
   case 'D':
   case 'd':
      return not_C;
   case 'H':
   case 'h':
      return not_G;
   case 'V':
   case 'v':
      return not_T;
   case 'N':
   case 'n':
      return Any;
   }
   BOOST_THROW_EXCEPTION(
            base_exception()
            << base_exception::msg_t("invalid character")
            << base_exception::str1_t(sanitize(c))
   );
   return Adenine; //never comes to this
}

inline bool is_ambiguous(const char c) {
   switch (c) {
   case 'A':
   case 'a':
   case 'C':
   case 'c':
   case 'G':
   case 'g':
   case 'T':
   case 't':
   case 'U':
   case 'u':
      return false;
   case 'R':
   case 'r':
   case 'Y':
   case 'y':
   case 'K':
   case 'k':
   case 'M':
   case 'm':
   case 'S':
   case 's':
   case 'W':
   case 'w':
   case 'B':
   case 'b':
   case 'D':
   case 'd':
   case 'H':
   case 'h':
   case 'V':
   case 'v':
   case 'N':
   case 'n':
      return true;
   }
   BOOST_THROW_EXCEPTION(
            base_exception()
            << base_exception::msg_t("invalid character")
            << base_exception::str1_t(sanitize(c))
   );
   return 0; //never comes to this
}

inline Nucleotide complement(const Nucleotide n) {
   BOOST_ASSERT((int)n < VDJ_PIPE_N_NUCLEOTIDES);
   switch (n) {
   case Adenine: return Thymine;
   case Cytosine: return Guanine;
   case Guanine: return Cytosine;
   case Thymine: return Adenine;
   case Uracil: return Adenine;
   case Purine: return Pyrimidine;
   case Pyrimidine: return Purine;
   case Ketone: return Amine;
   case Amine: return Ketone;
   case Strong: return Weak;
   case Weak: return Strong;
   case not_A: return not_T;
   case not_C: return not_G;
   case not_G: return not_C;
   case not_T: return not_A;
   case Any: return Any;
   }
   return Any; //never comes to this
}

inline char complement(const char c) {
   switch(c) {
   case 'A': return 'T';
   case 'a': return 't';
   case 'T': return 'A';
   case 't': return 'a';
   case 'C': return 'G';
   case 'c': return 'g';
   case 'G': return 'C';
   case 'g': return 'c';
   case 'U': return 'A';
   case 'u': return 'a';
   case 'R': return 'Y';
   case 'r': return 'y';
   case 'Y': return 'R';
   case 'y': return 'r';
   case 'K': return 'M';
   case 'k': return 'm';
   case 'M': return 'K';
   case 'm': return 'k';
   case 'S': return 'W';
   case 's': return 'w';
   case 'W': return 'S';
   case 'w': return 's';
   case 'N': return 'N';
   case 'n': return 'n';
   default:
      BOOST_THROW_EXCEPTION(
               base_exception()
               << base_exception::msg_t("invalid nucleotide symbol")
               << base_exception::str1_t(sanitize(c))
      );
   }
   return 0; //never comes to this
}

inline char to_capital(const Nucleotide n) {
   BOOST_ASSERT((int)n < VDJ_PIPE_N_NUCLEOTIDES);
   switch (n) {
   case Adenine: return 'A';
   case Cytosine: return 'C';
   case Guanine: return 'G';
   case Thymine: return 'T';
   case Uracil: return 'U';
   case Purine: return 'R';
   case Pyrimidine: return 'Y';
   case Ketone: return 'K';
   case Amine: return 'M';
   case Strong: return 'S';
   case Weak: return 'W';
   case not_A: return 'B';
   case not_C: return 'D';
   case not_G: return 'H';
   case not_T: return 'V';
   case Any: return 'N';
   }
   return 'X'; //never comes to this
}

inline char to_small(const Nucleotide n) {
   BOOST_ASSERT((int)n < VDJ_PIPE_N_NUCLEOTIDES);
   switch (n) {
   case Adenine: return 'a';
   case Cytosine: return 'c';
   case Guanine: return 'g';
   case Thymine: return 't';
   case Uracil: return 'u';
   case Purine: return 'r';
   case Pyrimidine: return 'y';
   case Ketone: return 'k';
   case Amine: return 'm';
   case Strong: return 's';
   case Weak: return 'w';
   case not_A: return 'b';
   case not_C: return 'd';
   case not_G: return 'h';
   case not_T: return 'v';
   case Any: return 'n';
   }
   return 'n'; //never comes to this
}

typedef boost::multi_array<int,2> scoring_matrix_t;

namespace detail{

template<int A, int I, int P, int U>
inline scoring_matrix_t make_sm() {
   scoring_matrix_t
   m(boost::extents[VDJ_PIPE_N_NUCLEOTIDES][VDJ_PIPE_N_NUCLEOTIDES]);

   static const int p = VDJ_PIPE_N_NUCLEOTIDES * VDJ_PIPE_N_NUCLEOTIDES;

   const int n[p] =
   {
            //                   AG CT GT AC CG AT ^A ^C ^G ^T
            // A  C  G  T  N  U  R  Y  K  M  S  W  B  D  H  V
      /* A */  A, I, I, I, U, I, P, I, I, P, I, P, I, P, P, P, // A
      /* C */  I, A, I, I, U, I, I, P, I, P, P, I, P, I, P, P, // C
      /* G */  I, I, A, I, U, I, P, I, P, I, P, I, P, P, I, P, // G
      /* T */  I, I, I, A, U, A, I, P, P, I, I, P, P, P, P, I, // T
      /* N */  U, U, U, U, U, U, U, U, U, U, U, U, U, U, U, U, // N
      /* U */  I, I, I, A, U, A, I, P, P, I, I, P, P, P, P, I, // U
      /* R */  P, I, P, I, U, I, P, I, P, P, P, P, P, P, P, P, // AG
      /* Y */  I, P, I, P, U, P, I, P, P, P, P, P, P, P, P, P, // CT
      /* K */  I, I, P, P, U, P, P, P, P, I, P, P, P, P, P, P, // GT
      /* M */  P, P, I, I, U, I, P, P, I, P, P, P, P, P, P, P, // AC
      /* S */  I, P, P, I, U, I, P, P, P, P, P, I, P, P, P, P, // CG
      /* W */  P, I, I, P, U, P, P, P, P, P, I, P, P, P, P, P, // AT
      /* B */  I, P, P, P, U, P, P, P, P, P, P, P, P, P, P, P, // CGT
      /* D */  P, I, P, P, U, P, P, P, P, P, P, P, P, P, P, P, // AGT
      /* H */  P, P, I, P, U, P, P, P, P, P, P, P, P, P, P, P, // ACT
      /* V */  P, P, P, I, U, I, P, P, P, P, P, P, P, P, P, P  // ACG
            // A  C  G  T  N  U  R  Y  K  M  S  W  B  D  H  V
            //                   AG CT GT AC CG AT ^A ^C ^G ^T
   };
   std::copy(n, n + p, m.data());
   return m;
}

}//namespace detail

template<int Match, int Mismatch, int Approximate, int Uncertain>
inline scoring_matrix_t const& scoring_matrix() {
   static const scoring_matrix_t m =
            detail::make_sm<Match,Mismatch,Approximate,Uncertain>();
   return m;
}

inline int identity(
         const Nucleotide n1,
         const Nucleotide n2,
         scoring_matrix_t const& sm = scoring_matrix<2,-2,1,0>()
) {
   BOOST_ASSERT(
            sm.shape()[0] == VDJ_PIPE_N_NUCLEOTIDES &&
            sm.shape()[1] == VDJ_PIPE_N_NUCLEOTIDES
   );
   BOOST_ASSERT((int)n1 < VDJ_PIPE_N_NUCLEOTIDES);
   BOOST_ASSERT((int)n2 < VDJ_PIPE_N_NUCLEOTIDES);
   BOOST_ASSERT(sm[n1][n2] == sm[n2][n1]);
   return sm[n1][n2];
}

struct Is_ambiguous : public boost::algorithm::predicate_facade<Is_ambiguous> {
   Is_ambiguous() {}
   typedef bool result_type;
   bool operator() (const char c) const {return vdj_pipe::is_ambiguous(c);}
};

}//namespace vdj_pipe
#endif /* NUCLEOTIDE_INDEX_HPP_ */