libsequence: RedundancyCom95.cc Source File

00001 /*
00002 
00003 Copyright (C) 2003-2009 Kevin Thornton, krthornt[]@[]uci.edu
00004 
00005 Remove the brackets to email me.
00006 
00007 This file is part of libsequence.
00008 
00009 libsequence is free software: you can redistribute it and/or modify
00010 it under the terms of the GNU General Public License as published by
00011 the Free Software Foundation, either version 3 of the License, or
00012 (at your option) any later version.
00013 
00014 libsequence is distributed in the hope that it will be useful,
00015 but WITHOUT ANY WARRANTY; without even the implied warranty of
00016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00017 GNU General Public License for more details.
00018 
00019 You should have received a copy of the GNU General Public License
00020 long with libsequence.  If not, see <http://www.gnu.org/licenses/>.
00021 
00022 */
00023 
00024 #include <cassert>
00025 #include <map>
00026 #include <cctype>
00027 #include <cstring>
00028 #include <algorithm>
00029 #include <functional>
00030 #include <Sequence/SeqProperties.hpp>
00031 #include <Sequence/SeqEnums.hpp>
00032 #include <Sequence/Translate.hpp>
00033 #include <Sequence/RedundancyCom95.hpp>
00034 
00035 
00036 namespace Sequence
00037 {
00038 #ifndef DOXYGEN_SKIP
00039   struct RedundancyCom95impl
00040   {
00041     explicit RedundancyCom95impl(const GeneticCodes &code);
00042 
00043     const GeneticCodes genetic_code;
00044       
00045     //matrices to be filled on the fly
00046     double firstNon[4][4][4];
00047     double first2S[4][4][4];
00048     double first2V[4][4][4];
00049     double thirdFour[4][4][4];
00050     double thirdNon[4][4][4];
00051     double third2S[4][4][4];
00052     double third2V[4][4][4];
00053     double l0_vals[4][4][4];
00054     double l2S_vals[4][4][4];
00055     double l2V_vals[4][4][4];
00056     double l4_vals[4][4][4];
00057     char intToNuc[4];
00058     //the map maps the nucleotides A,G,C,T onto array indices
00059     std::map<char,int> nucToInt;
00060     void FillFirstPositionCounts ();
00061     void FillThirdPositionCounts ();
00062     void FillLValues ();
00063     void codonPrecondition(const std::string & codon) throw (Sequence::SeqException);
00064   };
00065 
00066   RedundancyCom95impl::RedundancyCom95impl(const GeneticCodes &code) : genetic_code(code)
00067   {
00068     intToNuc[0] = 'A';
00069     intToNuc[1] = 'T';
00070     intToNuc[2] = 'G';
00071     intToNuc[3] = 'C';
00072     //the map is case-insensitive...
00073     nucToInt['A'] = 0;
00074     nucToInt['T'] = 1;
00075     nucToInt['G'] = 2;
00076     nucToInt['C'] = 3;
00077     nucToInt['a'] = 0;
00078     nucToInt['t'] = 1;
00079     nucToInt['g'] = 2;
00080     nucToInt['c'] = 3;
00081     FillFirstPositionCounts ();
00082     FillThirdPositionCounts ();
00083     FillLValues ();
00084   }    
00085 
00086   void RedundancyCom95impl::codonPrecondition(const std::string & codon) throw (Sequence::SeqException)
00093   {
00094     if ( codon.length() != 3 ||
00095          std::find_if(codon.begin(),codon.end(),ambiguousNucleotide()) != codon.end() )
00096       {
00097         throw(SeqException("Sequence::RedundancyCom95 -- precondition failed, invalid codon"));
00098       }
00099   }
00100 
00101   void
00102   RedundancyCom95impl::FillFirstPositionCounts ()
00103   {
00104     int i, j, k, l;
00105     std::string codon,mutation;
00106     int numN, numTsS, numTvS;   //numbers of changes of various types, N = nonsynon, S = synon
00107     int codonState;
00108     int type;           //type refers to Ts or Tv
00109     int numPossChanges, numPossTs, numPossTv;
00110     double FirstN, First2S, First2V;
00111     std::string codonTrans;
00112     std::string mutationTrans;
00113 
00114     bool S, N;
00115     //make every codon
00116     codon.resize(3);
00117     mutation.resize(3);
00118     for (i = Nucleotides (A); i <= Nucleotides (C); ++i)
00119       {
00120         for (j = Nucleotides (A); j <= Nucleotides (C); ++j)
00121           {
00122             for (k = Nucleotides (A); k <= Nucleotides (C); ++k)
00123               {
00124                 numN = numTsS = numTvS = 0;     //initialize
00125                 numPossChanges = numPossTs = numPossTv = 0;
00126                 S = N = 0;
00127                 codon[0] = intToNuc[i];
00128                 codon[1] = intToNuc[j];
00129                 codon[2] = intToNuc[k];
00130 
00131                 //now make all possible mutations at first position
00132                 codonState = i;
00133                 for (l = Nucleotides (A);
00134                      l <= Nucleotides (C); ++l)
00135                   {
00136                     if (l != codonState)
00137                       { //only include mutations
00138                         mutation[0] = intToNuc[l];
00139                         mutation[1] = codon[1];
00140                         mutation[2] = codon[2]; //the mutant codon is now defined
00141                         type = TsTv (intToNuc[codonState], mutation[0]);//find out if change is Ts or Tv
00142 
00143                         codonTrans = Translate (codon.begin(),
00144                                                 codon.end(), genetic_code);
00145                         mutationTrans = Translate (mutation.begin(),
00146                                                    mutation.end(),genetic_code);
00147 
00148                         //if neither codon is a stop
00149                         if ((codonTrans)[0] != '*'
00150                             && (mutationTrans)[0] != '*')
00151                           {
00152                             if(codonTrans == mutationTrans)
00153                               { //synonymous change
00154                                 S = 1;
00155                                 N = 0;
00156                               }
00157                             else
00158                               { //nonsynonymous change
00159                                 S = 0;
00160                                 N = 1;
00161                               }
00162                             if (type == Mutations(Ts))
00163                               ++numPossTs;
00164                             else if (type == Mutations(Tv))
00165                               ++numPossTv;
00166                             if (N == 1) //if change is nonsynonymous
00167                               ++numN;
00168                             else if (S == 1 && type == Mutations(Ts))
00169                               //if change is synonymous transition
00170                               ++numTsS;
00171                             else if (S == 1 && type == Mutations(Tv))
00172                               //if change is synonymous transversion
00173                               ++numTvS;
00174                           }
00175                       }
00176                   }
00177                 //process here
00178                 if ((numPossTs + numPossTv) == 0)
00179                   {     //stop codons
00180                     First2S = 0;
00181                     First2V = 0;
00182                     FirstN = 0;
00183                   }
00184                 else if ((numTsS + numTvS) == 0)        //non-degenerate sites
00185                   {
00186                     First2S = 0;
00187                     First2V = 0;
00188                     FirstN = 1.0;
00189                   }
00190                 else if (numPossTs != 0 && numPossTv != 0)
00191                   {     //special cases
00192                     if (double (numTsS) /
00193                         double (numPossTs) != 1.0
00194                         && double (numTvS) /
00195                         double (numPossTv) != 1.0)
00196                       { //fractional redundancy
00197                         First2S =
00198                           double (numTsS) /
00199                           double (numPossTs +
00200                                   numPossTv);
00201                         First2V =
00202                           double (numTvS) /
00203                           double (numPossTs +
00204                                   numPossTv);
00205                         FirstN = 1.0 - First2S -
00206                           First2V;
00207                       }
00208                     else
00209                       { //odd types of degeneracy
00210                         First2S =
00211                           double (numTsS) /
00212                           double (numPossTs);
00213                         First2V =
00214                           double (numTvS) /
00215                           double (numPossTv);
00216                         FirstN = 1.0 - First2S -
00217                           First2V;
00218                       }
00219                   }
00220                 else
00221                   {
00222                     if (numPossTs > 0)
00223                       First2S =
00224                         double (numTsS) /
00225                         double (numPossTs + numPossTv);
00226                     else
00227                       First2S = 0.0;
00228 
00229                     if (numPossTv > 0)
00230                       First2V =
00231                         double (numTvS) /
00232                         double (numPossTv + numPossTs);
00233                     else
00234                       First2V = 0.0;
00235                     FirstN = 1.0 - First2S - First2V;
00236                   }
00237                 //fill array
00238                 firstNon[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]]=FirstN;
00239                 first2S[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]]=First2S;
00240                 first2V[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]]=First2V;
00241               }
00242           }
00243       }
00244   }
00245 
00246   void
00247   RedundancyCom95impl::FillThirdPositionCounts ()
00248   {
00249     int i, j, k, l;
00250     std::string codon, mutation;        //the two codons
00251     int numN, numTsS, numTvS;   //numbers of changes of various types, N = nonsynon, S = synon
00252     int codonState;
00253     int type;           //type refers to Ts or Tv
00254     int numPossChanges, numPossTs, numPossTv;
00255     double ThirdN, Third2S, Third2V, Third4;
00256     std::string codonTrans;
00257     std::string mutationTrans;
00258     bool S, N;
00259     //make every codon
00260     codon.resize(3);
00261     mutation.resize(3);
00262     for (i = Nucleotides (A); i <= Nucleotides (C); ++i)
00263       {
00264         for (j = Nucleotides (A); j <= Nucleotides (C); ++j)
00265           {
00266             for (k = Nucleotides (A); k <= Nucleotides (C); ++k)
00267               {
00268                 numN = numTsS = numTvS = 0;     //initialize
00269                 numPossChanges = numPossTs = numPossTv = 0;
00270                 S = N = 0;
00271                 codon[0] = intToNuc[i];
00272                 codon[1] = intToNuc[j];
00273                 codon[2] = intToNuc[k];
00274                 //now make all possible mutations at first position
00275                 codonState = k;
00276 
00277                 for (l = Nucleotides (A);
00278                      l <= Nucleotides (C); ++l)
00279                   {
00280                     if (l != codonState)
00281                       { //only include mutations
00282                         mutation[0] = codon[0];
00283                         mutation[1] = codon[1];
00284                         mutation[2] = intToNuc[l];      //the mutant codon is now defined
00285                         type = TsTv (intToNuc[codonState], mutation[2]);        //find out if change is Ts or Tv
00286 
00287                         codonTrans = Translate (codon.begin(),
00288                                                 codon.end(),genetic_code);
00289                         mutationTrans = Translate (mutation.begin(),
00290                                                    mutation.end(),genetic_code);
00291 
00292                         if ((codonTrans)[0] != '*'
00293                             && (mutationTrans)[0] != '*')
00294                           {
00295                             if(codonTrans==mutationTrans)
00296                               {
00297                                 S = 1;
00298                                 N = 0;
00299                               }
00300                             else
00301                               {
00302                                 S = 0;
00303                                 N = 1;
00304                               }
00305                             if (type == Mutations(Ts))
00306                               ++numPossTs;
00307                             else if (type == Mutations(Tv))
00308                               ++numPossTv;
00309                             if (N)
00310                               ++numN;
00311                             else if (S && type == Mutations(Ts))
00312                               ++numTsS;
00313                             else if (S && type == Mutations(Tv))
00314                               ++numTvS;
00315                           }
00316                       }
00317                   }
00318                 //process here
00319                 ThirdN = Third2S = Third2V = Third4 = 0.0;
00320                 if (numPossTs + numPossTv == 0)
00321                   {     //stop codons
00322                     ThirdN = 0.0;
00323                     Third2S = 0.0;
00324                     Third2V = 0.0;
00325                     Third4 = 0.0;
00326                   }
00327                 else if (numTsS + numTvS == 0)
00328                   {     //nondegenerate
00329                     ThirdN = 1.0;
00330                     Third2S = 0.0;
00331                     Third2V = 0.0;
00332                     Third4 = 0.0;
00333                   }
00334                 else if (numTsS + numTvS == 3)
00335                   {     //4fold degenerate
00336                     ThirdN = 0.0;
00337                     Third2S = 0.0;
00338                     Third2V = 0.0;
00339                     Third4 = 1.0;
00340                   }
00341                 else if (numTsS == 0 || numTvS == 0)
00342                   {     //if there is non-degeneracy for one type
00343                     if (numPossTs > 0)
00344                       Third2S =
00345                         double (numTsS) /
00346                         double (numPossTs);
00347                     else
00348                       Third2S = 0.0;
00349                     if (numPossTv > 0)
00350                       Third2V =
00351                         double (numTvS) /
00352                         double (numPossTv);
00353                     else
00354                       Third2V = 0.0;
00355                     ThirdN = 1.0 - Third2S - Third2V;
00356                     Third4 = 0.0;
00357                   }
00358                 else if (numTsS > 0 && numTvS > 0
00359                          && (double (numTsS) /
00360                              double (numPossTs) != 1.0
00361                              || double (numTvS) /
00362                              double (numPossTv) != 1.0))
00363                   {     //fractional degeneracy
00364                     Third2S = 1.0 / 3.0;
00365                     Third2V = 1.0 / 3.0;
00366                     ThirdN = 1.0 / 3.0;
00367                     Third4 = 0.0;
00368                   }
00369                 thirdNon[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]]=ThirdN;
00370                 third2S[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]]=Third2S;
00371                 third2V[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]]=Third2V;
00372                 thirdFour[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]]=Third4;
00373               }
00374           }
00375       }
00376   }
00377 
00378   //the L values are just the sums of the site degeneracy values for each codon
00379   void
00380   RedundancyCom95impl::FillLValues ()
00381   {
00382     int i, j, k;
00383     std::string trans;
00384 
00385     std::string codon(3,'A');
00386     for (i = Nucleotides (A); i <= Nucleotides (C); ++i)
00387       for (j = Nucleotides (A); j <= Nucleotides (C); ++j)
00388         for (k = Nucleotides (A); k <= Nucleotides (C); ++k)
00389           {
00390             codon[0] = intToNuc[i];
00391             codon[1] = intToNuc[j];
00392             codon[2] = intToNuc[k];
00393 
00394             trans =Translate (codon.begin(),codon.end(),  genetic_code);
00395             if(strcmp(trans.c_str(),"*")!=0)
00396               { //if not a stop codon
00397                 l0_vals[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]] = 1. 
00398                   + firstNon[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]] 
00399                   + thirdNon[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]];
00400 
00401                 l2S_vals[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]] = 
00402                   first2S[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]] + 
00403                   third2S[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]];
00404 
00405                 l2V_vals[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]] =
00406                   first2V[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]] + 
00407                   third2V[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]];
00408 
00409                 l4_vals[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]] = 
00410                   thirdFour[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]];                
00411               }
00412             else
00413               {
00414                 l0_vals[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]] = 0.0;
00415                 l2S_vals[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]] = 0.0;
00416                 l2V_vals[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]] = 0.0;
00417                 l4_vals[nucToInt[codon[0]]][nucToInt[codon[1]]][nucToInt[codon[2]]] = 0.0;
00418               }
00419           }
00420   }
00421 #endif
00422   RedundancyCom95::RedundancyCom95 (const GeneticCodes &code):
00423     impl(std::auto_ptr<RedundancyCom95impl>(new RedundancyCom95impl(code)))
00427   {
00428   }
00429     
00430   RedundancyCom95::~RedundancyCom95(void)
00431   {
00432   }
00433 
00434   //the functions below are klugdy, but there is no real loss of efficieny,
00435   //so it is left alone for now
00436   double
00437   RedundancyCom95::FirstNon (const std::string & codon) const throw (Sequence::SeqException)
00443   {
00444     impl->codonPrecondition(codon);
00445     return impl->firstNon[impl->nucToInt[codon[0]]][impl->nucToInt[codon[1]]][impl->nucToInt[codon[2]]];
00446   }
00447 
00448   double
00449   RedundancyCom95::First2S (const std::string & codon) const throw (Sequence::SeqException)
00456   {
00457     impl->codonPrecondition(codon);
00458     return impl->first2S[impl->nucToInt[codon[0]]][impl->nucToInt[codon[1]]][impl->nucToInt[codon[2]]];
00459   }
00460 
00461   double
00462   RedundancyCom95::First2V (const std::string & codon) const throw (Sequence::SeqException)
00469   {
00470     impl->codonPrecondition(codon);
00471     return impl->first2V[impl->nucToInt[codon[0]]][impl->nucToInt[codon[1]]][impl->nucToInt[codon[2]]];
00472   }
00473 
00474   double
00475   RedundancyCom95::ThirdNon (const std::string & codon) const throw (Sequence::SeqException)
00481   {
00482     impl->codonPrecondition(codon);
00483     return impl->thirdNon[impl->nucToInt[codon[0]]][impl->nucToInt[codon[1]]][impl->nucToInt[codon[2]]];
00484   }
00485 
00486   double
00487   RedundancyCom95::ThirdFour (const std::string & codon) const throw (Sequence::SeqException)
00493   {
00494     impl->codonPrecondition(codon);
00495     return impl->thirdFour[impl->nucToInt[codon[0]]][impl->nucToInt[codon[1]]][impl->nucToInt[codon[2]]];
00496   }
00497 
00498   double
00499   RedundancyCom95::Third2S (const std::string & codon) const throw (Sequence::SeqException)
00505   {
00506     impl->codonPrecondition(codon);
00507     return impl->third2S[impl->nucToInt[codon[0]]][impl->nucToInt[codon[1]]][impl->nucToInt[codon[2]]];
00508   }
00509 
00510   double
00511   RedundancyCom95::Third2V (const std::string & codon) const throw (Sequence::SeqException)
00517   {
00518     impl->codonPrecondition(codon);
00519     return impl->third2V[impl->nucToInt[codon[0]]][impl->nucToInt[codon[1]]][impl->nucToInt[codon[2]]];
00520   }
00521 
00522   double
00523   RedundancyCom95::L0_vals (const std::string & codon) const throw (Sequence::SeqException)
00530   {
00531     impl->codonPrecondition(codon);
00532     return impl->l0_vals[impl->nucToInt[codon[0]]][impl->nucToInt[codon[1]]][impl->nucToInt[codon[2]]];
00533   }
00534 
00535   double
00536   RedundancyCom95::L2S_vals (const std::string & codon) const throw (Sequence::SeqException)
00543   {
00544     impl->codonPrecondition(codon);
00545     return impl->l2S_vals[impl->nucToInt[codon[0]]][impl->nucToInt[codon[1]]][impl->nucToInt[codon[2]]];
00546   }
00547 
00548   double
00549   RedundancyCom95::L2V_vals (const std::string & codon) const throw (Sequence::SeqException)
00556   {
00557     impl->codonPrecondition(codon);
00558     return impl->l2V_vals[impl->nucToInt[codon[0]]][impl->nucToInt[codon[1]]][impl->nucToInt[codon[2]]];
00559   }
00560 
00561   double
00562   RedundancyCom95::L4_vals (const std::string & codon) const throw (Sequence::SeqException)
00569   {
00570     impl->codonPrecondition(codon);
00571     return impl->l4_vals[impl->nucToInt[codon[0]]][impl->nucToInt[codon[1]]][impl->nucToInt[codon[2]]];
00572   }
00573 }