42 seqlen (seqa->length())
50 throw std::runtime_error (
"Sequence::Kimura80::Kimura80(): constructor called with two sequence objects of unequal length");
64 unsigned ungapped_sites = 0;
65 for (i = 0; i < seqlen; ++i)
70 if (std::toupper((*seq1)[i]) !=
71 std::toupper((*seq2)[i]))
73 auto type =
TsTv ((*seq1)[i], (*seq2)[i]);
75 if (type == Mutations::Ts)
79 else if (type == Mutations::Tv)
88 sites_compared = (ungapped_sites < seqlen) ? ungapped_sites : seqlen;
90 P = double (num_Ts) / double (sites_compared);
91 Q = double (num_Tv) / double (sites_compared);
94 if (fabs(1.0 - 2.0 * P - Q) > DBL_EPSILON)
96 divergence = -1.0 * 0.5 * std::log ((1.0 - 2.0 * P - Q)
97 * pow ((1 - 2.0 * Q), 0.5));
104 if (divergence <= 0.0-DBL_EPSILON)
117 if (!std::isfinite(divergence))
118 return std::numeric_limits<double>::quiet_NaN();
128 return (sites_compared);
Kimura80(const Sequence::Seq *seqa, const Sequence::Seq *seqb)
Mutations TsTv(const char &i, const char &j)
declaration of Sequence::Kimura80
class Sequence::Seq, an abstract base class for Sequences
Abstract interface to sequence objects.
The namespace in which this library resides.
size_type length(void) const
bool NotAGap(const char &c)
delcaration of routines for comparing DNA sequences This file declares a set of functions useful for ...
Definition of enumeration types.