51 struct Comeron95::Com95impl
53 double Qs, Bs, Qa, Ba, A2S, A4, As, A2V, A0, Aa,
54 q0, q2S, q2V, q4, p0, p2S, p2V, p4,Ka,Ks;
55 std::unique_ptr<RedundancyCom95> sitesObj;
58 const WeightingScheme2 *_weights2,
59 const WeightingScheme3 *_weights3,
61 void omega (
const Sites * s,
64 Qs(0.), Bs(0.), Qa(0.), Ba(0.), A2S(0.), A4(0.), As(0.), A2V(0.), A0(0.), Aa(0.),
65 q0(0.), q2S(0.), q2V(0.), q4(0.), p0(0.), p2S(0.), p2V(0.), p4(0.),Ka(0.),Ks(0.),
66 sitesObj(
std::unique_ptr<RedundancyCom95>(new RedundancyCom95(__code))),
113 double ka (
void)
const;
114 double ks (
void)
const;
115 double ratio (
void)
const;
116 double P0 (
void)
const;
117 double P2S (
void)
const;
118 double P2V (
void)
const;
119 double P4 (
void)
const;
120 double Q0 (
void)
const;
121 double Q2S (
void)
const;
122 double Q2V (
void)
const;
123 double Q4 (
void)
const;
124 double as (
void)
const;
125 double aa (
void)
const;
126 double bs (
void)
const;
127 double ba (
void)
const;
128 double L0 (
const Sites *)
const;
129 double L2S (
const Sites * )
const;
130 double L2V (
const Sites *)
const;
131 double L4 (
const Sites *)
const;
134 Comeron95::Comeron95(
GeneticCodes code ) : impl(
std::unique_ptr<Com95impl>(new Com95impl(code)))
138 Comeron95::~Comeron95(){}
163 return this->
operator()(seqa,seqb,&w2,&w3,max);
172 Sites s(seqa,seqb,(*impl->sitesObj),maxdiffs);
173 impl->diverge(seqa,seqb,weights2,weights3,maxdiffs);
174 impl->omega(&s,seqa,seqb);
175 return Com95_t({{impl->ka(),
197 void Comeron95::Com95impl::omega (
const Sites * s,
212 Qs = (q2V + q4) / (s->
L2V() + s->
L4());
214 if (!std::isfinite (Qs))
217 Bs = (-0.5) * log (1.0 - (2.0 * Qs));
219 Qa = (q0 + q2S) / (s->
L0() + s->
L2S());
221 if (!std::isfinite (Qa))
224 Ba = (-0.5) * std::log (1.0 - (2.0 * Qa));
225 if (!std::isfinite (Ba))
231 std::unique_ptr<Kimura80> K80(
new Kimura80 (&seqobj1, &seqobj2));
237 double P2S_site = p2S / s->
L2S();
238 double P2V_site = p2V / s->
L2V();
239 double P0_site = p0 / s->
L0();
240 double Q0_site = q0 / s->
L0();
241 double P4_site = p4 / s->
L4();
242 double Q4_site = q4 / s->
L4();
244 log1 = std::log (1.0 - (2.0 * P2S_site) - Qa);
245 log2 = std::log (1.0 - (2.0 * Qa));
247 if (!std::isfinite (log1))
249 if (!std::isfinite (log2))
252 A2S = (-0.5) * log1 + (0.25) * log2;
254 log1 = std::log (1.0 - (2.0 * P4_site) - Q4_site);
255 log2 = std::log (1.0 - (2.0 * Q4_site));
257 if (!std::isfinite (log1))
259 if (!std::isfinite (log2))
262 A4 = (-0.5) * log1 + (0.25) * log2;
264 As = (s->
L2S() * A2S + s->
L4() * A4) / (s->
L2S() +
267 log1 = std::log (1.0 - (2.0 * P2V_site) - Qs);
268 log2 = std::log (1.0 - (2.0 * Qs));
270 if (!std::isfinite (log1))
272 if (!std::isfinite (log2))
275 A2V = (-0.5) * log1 + (0.25) * log2;
277 log1 = std::log (1.0 - (2.0 * P0_site) - Q0_site);
278 log2 = std::log (1.0 - (2.0 * Q0_site));
279 if (!std::isfinite (log1))
281 if (!std::isfinite (log2))
284 A0 = (-0.5) * log1 + (0.25) * log2;
286 Aa = (s->
L2V() * A2V + s->
L0() * A0) / (s->
L2V() +
301 if (!std::isfinite (Ks))
302 Ks = std::numeric_limits<double>::quiet_NaN();
303 if (!std::isfinite (Ka))
304 Ka = std::numeric_limits<double>::quiet_NaN();
307 void Comeron95::Com95impl::diverge (
const Sequence::Seq & seq1,
319 size_t length = seq1.
length();
320 q0= q2S= q2V= q4= p0= p2S= p2V= p4 = 0.;
322 std::string codon1, codon2;
326 for (i = 0; i <= (length - 3); i += 3)
328 for (j = 0; j <= 2; ++j)
331 codon1[j] = char(std::toupper(seq1[i + j]));
332 codon2[j] = char(std::toupper(seq2[i + j]));
341 int ndiff =
NumDiffs (codon1, codon2);
346 Single(*sitesObj.get(), codon1, codon2);
348 p2S += Single.
P2S ();
349 p2V += Single.
P2V ();
352 q2S += Single.
Q2S ();
353 q2V += Single.
Q2V ();
373 Double(*sitesObj.get(), codon1, codon2,weights2);
375 p2S += Double.
P2S ();
376 p2V += Double.
P2V ();
379 q2S += Double.
Q2S ();
380 q2V += Double.
Q2V ();
383 else if (ndiff == 3 && maxdiffs > 2)
386 Triple(*sitesObj.get(),codon1, codon2,weights3);
388 p2S += Triple.
P2S ();
389 p2V += Triple.
P2V ();
392 q2S += Triple.
Q2S ();
393 q2V += Triple.
Q2V ();
402 if (!std::isfinite (p0))
404 if (!std::isfinite (p2S) )
406 if (!std::isfinite (p2V) )
408 if (!std::isfinite (p4) )
410 if (!std::isfinite (q0) )
412 if (!std::isfinite (q2S) )
414 if (!std::isfinite (q2V) )
416 if (!std::isfinite (q4) )
420 double Comeron95::Com95impl::L0 (
const Sites * sites)
const 427 double Comeron95::Com95impl::L2S (
const Sites * sites)
const 434 double Comeron95::Com95impl::L2V (
const Sites * sites)
const 441 double Comeron95::Com95impl::L4 (
const Sites * sites)
const 449 double Comeron95::Com95impl::as (
void)
const 454 if (!std::isfinite ( As))
455 return std::numeric_limits<double>::quiet_NaN();
458 double Comeron95::Com95impl::aa (
void)
const 463 if (!std::isfinite ( Aa))
464 return std::numeric_limits<double>::quiet_NaN();
467 double Comeron95::Com95impl::bs (
void)
const 472 if (!std::isfinite ( Bs))
473 return std::numeric_limits<double>::quiet_NaN();
476 double Comeron95::Com95impl::ba (
void)
const 481 if (!std::isfinite ( Ba))
482 return std::numeric_limits<double>::quiet_NaN();
486 double Comeron95::Com95impl::ratio(
void)
const 492 if ( !std::isfinite(Ka) || !std::isfinite(Ks) || fabs(Ks-0.) <= DBL_EPSILON)
493 return std::numeric_limits<double>::quiet_NaN();
497 double Comeron95::Com95impl::ka (
void)
const 505 double Comeron95::Com95impl::ks (
void)
const 514 double Comeron95::Com95impl::P0 (
void)
const 521 double Comeron95::Com95impl::P2S (
void)
const 528 double Comeron95::Com95impl::P2V (
void)
const 535 double Comeron95::Com95impl::P4 (
void)
const 542 double Comeron95::Com95impl::Q0 (
void)
const 549 double Comeron95::Com95impl::Q2S (
void)
const 556 double Comeron95::Com95impl::Q2V (
void)
const 563 double Comeron95::Com95impl::Q4 (
void)
const
declaration of Sequence::Kimura80
abstract interface to weighting schemes when codons differ at 2 positions
class Sequence::Seq, an abstract base class for Sequences
Abstract interface to sequence objects.
declaration of Sequence::Comeron95 to calculate Ka/Ks
Kimura's 2-parameter distance.
Com95_t operator()(const Sequence::Seq &seqa, const Sequence::Seq &seqb, int maxdiffs=3)
Tests if a character is in the set A,G,C,T.
used by Sequence::Comeron95, class Sequence::TwoSubs calculates divergence between codons that differ...
class Sequence::Sites calculates the length of a pairwise alignment between coding sequences in terms...
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)
Deal with codons differing at 1 position.
declaration of classes to weight codons by Grantham distance (i.e. for Sequence::Comeron95). Declares Sequence::GranthamWeights2 and Sequence::GranthamWeights3
class Sequence::RedundancyCom95
used by Sequence::Comeron95, class Sequence::ThreeSubs calculates divergence between codons that diff...
Grantham's distances (Sequence::Grantham)
Weights paths by Grantham's distances for codons differing at 3 sites.
used by Sequence::Comeron95, class Sequence::SingleSub calculates divergence between codons that diff...
abstract interface to weighting schemes when codons differ at 3 positions
Calculate length statistics for divergence calculations.
size_type length(void) const
bool Different(const std::string &seq1, const std::string &seq2, const bool &skip_missing=true, const bool &nucleic_acid=true)
Weights paths by Grantham's distances for codons differing at 2 sites.
Deal with codons differing at all 3 positions.
delcaration of routines for comparing DNA sequences This file declares a set of functions useful for ...
Deal with codons differing at 2 positions.