35 #ifndef DOXYGEN_SKIP //doxygen should skip this 36 class shortestPath::shortestPathImpl
43 std::vector<std::string> _path;
44 std::pair<double,shortestPath::pathType>process_path(
const std::string &intermediate,
46 std::pair<double,shortestPath::pathType>process_path2(
const std::string &intermediate1,
47 const std::string &intermediate2,
49 shortestPathImpl(
const std::string &codon1,
50 const std::string &codon2,
54 shortestPath::shortestPathImpl::shortestPathImpl(
const std::string &codon1,
55 const std::string &codon2,
61 if(codon1.length()!=3 || codon2.length() != 3)
63 throw std::runtime_error(
"Codons are not both of length 3");
68 _path.push_back(codon1);
75 _type = shortestPath::pathType::NONE;
77 else if (t1 !=
"X" && t2 !=
"X" &&
78 (std::find_if(codon1.begin(),codon1.end(),
79 ambiguousNucleotide()) == codon1.end()) &&
80 (std::find_if(codon2.begin(),codon2.end(),
81 ambiguousNucleotide()) == codon2.end()) )
85 _distance = distances(t1[0],t2[0]);
88 _type = shortestPath::pathType::N;
92 _type = shortestPath::pathType::S;
100 std::vector< std::pair<double, shortestPath::pathType> >
101 paths(2,std::pair<double, shortestPath::pathType>(0.,shortestPath::pathType::AMBIG));
104 paths[0] = process_path(intermediates[0],code);
105 paths[1] = process_path(intermediates[1],code);
106 if (paths[0].first < paths[1].first)
108 _type = paths[0].second;
109 _distance = paths[0].first;
110 _path.push_back(intermediates[0]);
114 _type = paths[1].second;
115 _distance = paths[1].first;
116 _path.push_back(intermediates[1]);
119 else if (ndiffs == 3)
122 std::vector< std::pair<double, shortestPath::pathType> >
123 paths(6,std::pair<double, shortestPath::pathType>(0.,shortestPath::pathType::AMBIG));
127 paths[0] = process_path2(intermediates[0],intermediates[1],code);
128 paths[1] = process_path2(intermediates[0],intermediates[2],code);
129 paths[2] = process_path2(intermediates[3],intermediates[4],code);
130 paths[3] = process_path2(intermediates[3],intermediates[5],code);
131 paths[4] = process_path2(intermediates[6],intermediates[7],code);
132 paths[5] = process_path2(intermediates[6],intermediates[8],code);
135 double shortest_path = paths[0].first;
136 for(
unsigned i = 0 ; i < 6 ; ++i)
138 if (paths[i].first < shortest_path)
141 shortest_path=paths[i].first;
144 _type = paths[shortest].second;
145 _distance = paths[shortest].first;
149 _path.push_back(intermediates[0]);
150 _path.push_back(intermediates[1]);
153 _path.push_back(intermediates[0]);
154 _path.push_back(intermediates[2]);
157 _path.push_back(intermediates[3]);
158 _path.push_back(intermediates[4]);
161 _path.push_back(intermediates[3]);
162 _path.push_back(intermediates[5]);
165 _path.push_back(intermediates[6]);
166 _path.push_back(intermediates[7]);
169 _path.push_back(intermediates[6]);
170 _path.push_back(intermediates[8]);
177 _type = shortestPath::pathType::AMBIG;
179 _path.push_back(codon2);
182 std::pair<double,shortestPath::pathType>
183 shortestPath::shortestPathImpl::process_path(
const std::string &intermediate,
189 double d = ( distances(t1[0],tint[0]) + distances(tint[0],t2[0]) );
190 if ( (t1 != tint) && (tint != t2) )
192 t = shortestPath::pathType::NN;
194 else if ( ((t1 != tint) && (tint == t2)) ||
195 ((t1 == tint) && (tint != t2)) )
197 t = shortestPath::pathType::SN;
199 else if ( t1==tint && tint==t2 )
201 t = shortestPath::pathType::SS;
203 return std::make_pair(d,t);
206 std::pair<double,shortestPath::pathType>
207 shortestPath::shortestPathImpl::process_path2(
const std::string &intermediate1,
208 const std::string &intermediate2,
212 intermediate1.end(),code);
214 intermediate2.end(),code);
216 double d = (distances(t1[0],tint1[0]) +
217 distances(tint1[0],tint2[0]) +
218 distances(tint2[0],t2[0]));
220 bool a = (t1==tint1),
224 if ( !a && !b && !c )
226 t = shortestPath::pathType::NNN;
228 else if ( (a && b && !c) ||
232 t = shortestPath::pathType::SSN;
234 else if ( (!a && !b && c) ||
238 t = shortestPath::pathType::SNN;
241 else if (a && b && c)
243 t = shortestPath::pathType::SSS;
245 return std::make_pair(d,t);
250 const std::string &codon2,
265 impl = std::unique_ptr<shortestPath::shortestPathImpl>(
new shortestPath::shortestPathImpl(codon1,codon2,code));
267 catch(std::runtime_error &e)
271 catch(std::exception &e)
273 throw (std::runtime_error(e.what()));
277 throw (std::runtime_error(
"caught exception of unknown type"));
281 shortestPath::~shortestPath()
301 return impl->_path.begin();
311 return impl->_path.end();
319 return impl->_distance;
323 const std::string &codon2,
353 std::pair<unsigned,unsigned> rv;
360 case shortestPath::pathType::S :
362 rv = std::make_pair(1,0);
365 case shortestPath::pathType::N :
367 rv = std::make_pair(0,1);
370 case shortestPath::pathType::SS :
372 rv = std::make_pair(2,0);
375 case shortestPath::pathType::SN :
377 rv = std::make_pair(1,1);
380 case shortestPath::pathType::NN :
382 rv = std::make_pair(0,2);
385 case shortestPath::pathType::SSS :
387 rv = std::make_pair(3,0);
390 case shortestPath::pathType::SSN :
392 rv = std::make_pair(2,1);
395 case shortestPath::pathType::SNN :
397 rv = std::make_pair(1,2);
400 case shortestPath::pathType::NNN :
402 rv = std::make_pair(0,3);
405 case shortestPath::pathType::NONE :
407 rv = std::make_pair(0,0);
410 case shortestPath::pathType::AMBIG :
418 catch (std::runtime_error &e)
425 std::pair<unsigned,shortestPath::pathType>
diffType(
const std::string &codon1,
426 const std::string &codon2,
452 if (type == shortestPath::pathType::AMBIG)
456 unsigned site=0,ndiffs=0;
457 for(
unsigned i = 0 ; i < 3 ; ++i)
459 if (codon1[i] != codon2[i])
469 return std::make_pair(site,type);
471 catch (std::runtime_error &e)
477 std::tuple<shortestPath::pathType,shortestPath::pathType,shortestPath::pathType>
479 const std::string &codon2,
499 std::string t1(codon1),t2(codon2);
500 for(
unsigned i = 0; i < 3 ; ++i)
503 std::swap(t1[i],t2[i]);
504 std::pair<unsigned,shortestPath::pathType> dt1 =
diffType(t1,codon1,code);
505 std::pair<unsigned,shortestPath::pathType> dt2 =
diffType(t2,codon2,code);
509 p[i] = shortestPath::pathType::AMBIG;
511 std::swap(t1[i],t2[i]);
513 return std::make_tuple(p[0],p[1],p[2]);
515 catch(std::runtime_error &e)
std::string Translate(std::string::const_iterator beg, std::string::const_iterator end, Sequence::GeneticCodes genetic_code=GeneticCodes::UNIVERSAL, const char &gapchar='-')
std::pair< unsigned, shortestPath::pathType > diffType(const std::string &codon1, const std::string &codon2, const Sequence::GeneticCodes &code=GeneticCodes::UNIVERSAL)
const_iterator begin() const
double path_distance() const
The namespace in which this library resides.
int NumDiffs(const std::string &seq1, const std::string &seq2, const bool &skip_missing=true, const bool &nucleic_acid=true)
shortestPath(const std::string &codon1, const std::string &codon2, const Sequence::GeneticCodes &code=GeneticCodes::UNIVERSAL)
Inter2_t Intermediates2(const std::string &codon1, const std::string &codon2)
Calculate the intermediate codons between a pair of codons diverged at 2 positions.
declares Sequence::Translate,a function to translate CDS sequences into peptide sequences ...
Grantham's distances (Sequence::Grantham)
const unsigned SEQMAXUNSIGNED
declarations of Sequence::Intermediates2 and Sequence::Intermediates3
Calculate shortest path between 2 codons.
Grantham's distances.
const_iterator end() const
std::pair< unsigned, unsigned > mutsShortestPath(const std::string &codon1, const std::string &codon2, const Sequence::GeneticCodes &code=GeneticCodes::UNIVERSAL)
Inter3_t Intermediates3(const std::string &codon1, const std::string &codon2)
Calculate the intermediate codons between a pair of codons diverged at 3 positions.
delcaration of routines for comparing DNA sequences This file declares a set of functions useful for ...
std::tuple< shortestPath::pathType, shortestPath::pathType, shortestPath::pathType > diffTypeMulti(const std::string &codon1, const std::string &codon2, const Sequence::GeneticCodes &code=GeneticCodes::UNIVERSAL)
Routines to find the shortest distance between any 2 codons, using Grantham's distance. Declares the class Sequence::shortestPath, and the functions Sequence::mutsShortestPath and Sequence::diffType.