alg/alphabet.cpp

   1 #include <stdexcept>
   2
   3 #include "mussa_exceptions.hpp"
   4 #include "alg/alphabet.hpp"
   5
   6 // some standard dna alphabets
   7 // Include nl (\012), and cr (\015) to make sequence parsing eol
   8 // convention independent.
   9 const char *Alphabet::reduced_dna_cstr = "AaCcGgTtNn\012\015";
  10 const char *Alphabet::reduced_dna_reverse_cstr = "TtGgCcAaNn\012\015";
  11 const char *Alphabet::reduced_rna_cstr = "AaCcGgUuNn\012\015";
  12 const char *Alphabet::reduced_rna_reverse_cstr = "UuGgCcAaNn\012\015";
  13 const char *Alphabet::reduced_nucleic_cstr = "AaCcGgTtUuNn\012\015";
  14 const char *Alphabet::reduced_nucleic_reverse_cstr = "TtGgCcAaAaNn\012\015";
  15 //! this is the general iupac alphabet for nucleotides
  16 const char *Alphabet::dna_cstr =
  17   "AaCcGgTtRrYySsWwKkMmBbDdHhVvNnXx-~.?\012\015";
  18 const char *Alphabet::dna_reverse_cstr =
  19   "TtGgCcAaYyRrSsWwMmKkVvHhDdBbNnXx-~.?\012\015";
  20 const char *Alphabet::rna_cstr =
  21   "AaCcGgUuRrYySsWwKkMmBbDdHhVvNnXx-~.?\012\015";
  22 const char *Alphabet::rna_reverse_cstr =
  23   "UuGgCcAaYyRrSsWwMmKkVvHhDdBbNnXx-~.?\012\015";
  24 const char *Alphabet::nucleic_cstr =
  25   "AaCcGgTtUuRrYySsWwKkMmBbDdHhVvNnXx-~.?\012\015";
  26 const char *Alphabet::nucleic_reverse_cstr =
  27   "TtGgCcAaAaYyRrSsWwMmKkVvHhDdBbNnXx-~.?\012\015";
  28 //! the protein alphabet
  29 const char *Alphabet::protein_cstr =
  30   "AaCcDdEeFfGgHhIiKkLlMmNnPpQqRrSsTtVvWwYy\012\015";
  31 const char *Alphabet::empty_cstr = "";
  32
  33 const Alphabet& Alphabet::reduced_dna_alphabet() {
  34   static Alphabet *a = new Alphabet(reduced_dna_cstr, reduced_dna_reverse_cstr);
  35   return *a;
  36 }
  37 const Alphabet& Alphabet::reduced_rna_alphabet() {
  38   static Alphabet *a = new Alphabet(reduced_rna_cstr, reduced_rna_reverse_cstr);
  39   return *a;
  40 }
  41 const Alphabet& Alphabet::reduced_nucleic_alphabet() {
  42   static Alphabet *a = new Alphabet(reduced_nucleic_cstr, reduced_nucleic_reverse_cstr);
  43   return *a;
  44 }
  45 const Alphabet& Alphabet::dna_alphabet() {
  46   static Alphabet *a = new Alphabet(dna_cstr, dna_reverse_cstr);
  47   return *a;
  48 }
  49 const Alphabet& Alphabet::rna_alphabet() {
  50   static Alphabet *a = new Alphabet(rna_cstr, rna_reverse_cstr);
  51   return *a;
  52 }
  53 const Alphabet& Alphabet::nucleic_alphabet() {
  54   static Alphabet *a = new Alphabet(nucleic_cstr, nucleic_reverse_cstr);
  55   return *a;
  56 }
  57 const Alphabet& Alphabet::protein_alphabet() {
  58   static Alphabet *a = new Alphabet(protein_cstr, protein_cstr);
  59   return *a;
  60 }
  61 const Alphabet& Alphabet::empty_alphabet() {
  62   static Alphabet *a = new Alphabet(empty_cstr, empty_cstr);
  63   return *a;
  64 }
  65
  66 Alphabet::Alphabet(const char *a, const char *reverse_alphabet) :
  67   alphabet(a),
  68   complement_map(create_complement_map(reverse_alphabet))
  69 {
  70   alphabet_set.insert(alphabet.begin(), alphabet.end());
  71 }
  72
  73 void Alphabet::assign(const Alphabet& a)
  74 {
  75   alphabet = a.alphabet;
  76   complement_map = a.complement_map;
  77   alphabet_set.clear();
  78   alphabet_set.insert(alphabet.begin(), alphabet.end());
  79 }
  80
  81 bool operator==(const Alphabet& x, const Alphabet& y)
  82 {
  83   return x.alphabet == y.alphabet;
  84 }
  85
  86 std::ostream& operator<<(std::ostream& out, const Alphabet& a)
  87 {
  88   out << a.alphabet;
  89 }
  90
  91 bool Alphabet::exists(const char c) const
  92 {
  93   return (alphabet_set.find(c) != alphabet_set.end());
  94 }
  95
  96 const Alphabet& Alphabet::get_alphabet(AlphabetRef alpha)
  97 {
  98   switch (alpha) {
  99     case ::reduced_dna_alphabet:
 100       return Alphabet::reduced_dna_alphabet();
 101     case ::reduced_rna_alphabet:
 102       return Alphabet::reduced_rna_alphabet();
 103     case ::reduced_nucleic_alphabet:
 104       return Alphabet::reduced_nucleic_alphabet();
 105     case ::dna_alphabet:
 106       return Alphabet::dna_alphabet();
 107     case ::rna_alphabet:
 108       return Alphabet::rna_alphabet();
 109     case ::nucleic_alphabet:
 110       return Alphabet::nucleic_alphabet();
 111     case ::protein_alphabet:
 112       return Alphabet::protein_alphabet();
 113     case ::empty_alphabet:
 114       return Alphabet::empty_alphabet();
 115     default:
 116       throw std::runtime_error("unrecognized alphabet type");
 117       break;
 118   }
 119 }
 120
 121 std::string Alphabet::create_complement_map(const std::string& ra) const
 122 {
 123   if (ra.size() == 0) {
 124     // reverse complement is undefined
 125     throw mussa_error("no reverse complement");
 126   }
 127   std::string rc_map(256,'~');
 128
 129   for(int i = 0; alphabet[i] and ra[i]; ++i) {
 130     rc_map[ alphabet[i] ] = ra[i];
 131   }
 132   return rc_map;
 133 }
 134
 135 boost::shared_ptr<std::string> Alphabet::reverse_complement(const std::string &seq) const
 136 {
 137   boost::shared_ptr<std::string> new_seq(new std::string());
 138
 139   new_seq->reserve(seq.size());
 140
 141   for(std::string::const_reverse_iterator seq_i = seq.rbegin();
 142       seq_i != seq.rend();
 143       ++seq_i)
 144   {
 145     new_seq->append(1, complement_map[ *seq_i ]);
 146   }
 147   return new_seq;
 148 }