25 #include <Sequence/PolySNPimpl.hpp> 43 : PolySNP(data, false, 0, true)
48 rep->_NumPoly = NumPoly();
51 PolySIM::~PolySIM(
void) {}
54 PolySIM::ThetaPi(
void)
const 60 std::lock_guard<std::mutex> lock(rep->instance_lock);
61 if (rep->_know_pi ==
false)
63 double Pi = 0.0, nsam = double(rep->_nsam);
64 for (std::vector<stateCounter>::const_iterator i
65 = rep->_counts.begin();
66 i != rep->_counts.end(); ++i)
68 Pi += (2.0 * i->one * (nsam - i->one))
69 / (nsam * (nsam - 1.));
78 PolySIM::ThetaW(
void)
const 85 return (rep->_NumPoly > 0) ? (double(rep->_NumPoly) / a_sub_n()) : 0.;
89 PolySIM::ThetaH(
void)
const 96 double H = 0.0, nsam = double(rep->_nsam);
98 for (std::vector<stateCounter>::const_iterator i
99 = rep->_counts.begin();
100 i != rep->_counts.end(); ++i)
103 ? (2.0 * i->one * i->one) / (nsam * (nsam - 1.0))
110 PolySIM::ThetaL(
void)
const 117 if (rep->_NumPoly == 0)
119 const char state =
'1';
123 double thetal = 0.0, nsam = double(rep->_nsam);
125 for (site = 0; site < rep->_data->numsites(); ++site)
127 for (seq = 0, num_changes = 0; seq < unsigned(nsam); ++seq)
130 += (*rep->_data)[seq][site] == state ? 1 : 0;
132 if (num_changes < nsam)
134 double nc = double(num_changes);
135 thetal += nc / (nsam - 1.0);
142 PolySIM::TajimasD(
void)
const 144 if (rep->_NumPoly == 0)
145 return std::numeric_limits<double>::quiet_NaN();
146 double Pi = ThetaPi();
148 return ((Pi - W) / Dnominator());
152 PolySIM::Hprime(
const bool &likeThorntonAndolfatto)
const 158 if (rep->_NumPoly == 0)
159 return std::numeric_limits<double>::quiet_NaN();
160 double pi = ThetaPi();
161 double theta = ThetaW();
162 double thetal = ThetaL();
163 double a = a_sub_n();
164 double b = b_sub_n();
165 double b_n_plus1 = b_sub_n_plus1();
166 double S = rep->_NumPoly;
167 double thetasq = (likeThorntonAndolfatto ==
false)
168 ? S * (S - 1) / (a * a + b)
171 double vThetal = (rep->_nsam * theta) / (2.0 * (rep->_nsam - 1.0))
172 + (2.0 * pow(rep->_nsam / (rep->_nsam - 1.0), 2.0)
178 = (3.0 * rep->_nsam * (rep->_nsam + 1.0) * theta
179 + 2.0 * (rep->_nsam * rep->_nsam + rep->_nsam + 3.0) * thetasq)
180 / (9 * rep->_nsam * (rep->_nsam - 1.0));
182 double cov = (rep->_nsam + 1.0) * theta / (3.0 * (rep->_nsam - 1.0))
183 + (7.0 * rep->_nsam * rep->_nsam + 3.0 * rep->_nsam - 2.0
184 - 4.0 * rep->_nsam * (rep->_nsam + 1.0) * b_n_plus1)
185 * thetasq / (2.0 * pow((rep->_nsam - 1.0), 2.0));
188 Hpr /= pow((vThetal + vPi - 2.0 * cov), 0.5);
193 PolySIM::Dnominator(
void)
const 195 if (rep->_NumPoly == 0)
196 return std::numeric_limits<double>::quiet_NaN();
197 double S = rep->_NumPoly;
198 double a1, a2, b1, b2, c1, c2, e1, e2;
202 b1 = (rep->_nsam + 1.0) / (3.0 * (rep->_nsam - 1.0));
203 b2 = (2.0 * (pow(rep->_nsam, 2.0) + rep->_nsam + 3.0))
204 / (9.0 * rep->_nsam * (rep->_nsam - 1.0));
206 c2 = b2 - (rep->_nsam + 2.0) / (a1 * rep->_nsam) + a2 / pow(a1, 2.0);
208 e2 = c2 / (pow(a1, 2.0) + a2);
209 double denominator = pow((e1 * S + e2 * S * (S - 1.0)), 0.5);
210 return (denominator);
214 PolySIM::HudsonsHaplotypeTest(
const int &subsize,
const int &subss)
const 226 int *subslist, i, seq;
231 subslist =
new int[subsize];
233 for (i = 0; i < subsize; ++i)
239 = poly(subslist,
int(rep->_nsites), subsize, subss, &seq);
242 = nextsample(subslist, subsize,
int(rep->_nsam), seq);
262 PolySIM::poly(
int *subslist,
const int &ss,
const int &subsize,
263 const int &subss,
int *seq)
const 273 int npoly = 0, i, j, *polyvector;
275 polyvector =
new int[ss];
276 for (i = 0; i < ss; ++i)
279 for (i = 1; i < subsize; ++i)
280 for (j = 0; j < ss; ++j)
281 if ((*rep->_data)[PolyTable::size_type(subslist[i])]
282 [std::string::size_type(j)]
283 != (*rep->_data)[PolyTable::size_type(subslist[0])]
284 [std::string::size_type(j)])
286 if (polyvector[j] == 0)
290 if (npoly == subss + 1)
299 PolySIM::nextsample(
int *subslist,
const int &subsize,
const int &nsam,
312 if (subslist[seq] > nsam - subsize + seq)
319 return (nextsample(subslist, subsize, nsam, seq));
322 for (i = seq + 1; i < subsize; i++)
323 subslist[i] = subslist[i - 1] + 1;
328 PolySIM::FuLiD(
void)
const 330 if (rep->_NumPoly == 0)
331 return std::numeric_limits<double>::quiet_NaN();
332 double ExternalMutations = double(NumExternalMutations());
333 double NumMut = double(NumMutations());
334 double a = a_sub_n();
335 double b = b_sub_n();
336 double c = c_sub_n();
338 double vD = 1.0 + (pow(a, 2.0) / (b + pow(a, 2.0))
339 * (c - (rep->_nsam + 1.0) / (rep->_nsam - 1.0)));
340 double uD = a - 1.0 - vD;
341 double D = NumMut - a * double(ExternalMutations);
342 D /= pow((uD * NumMut + vD * pow(NumMut, 2.0)), 0.5);
347 PolySIM::FuLiF(
void)
const 349 if (rep->_NumPoly == 0)
350 return std::numeric_limits<double>::quiet_NaN();
351 double Pi = ThetaPi();
352 double NumMut = double(NumMutations());
353 double ExternalMutations = double(NumExternalMutations());
354 double a = a_sub_n();
355 double a_n_plus1 = a_sub_n_plus1();
356 double b = b_sub_n();
357 double c = c_sub_n();
359 + 2.0 * (pow(rep->_nsam, 2.0) + rep->_nsam + 3.0)
360 / (9.0 * rep->_nsam * (double(rep->_nsam - 1.0)));
361 vF -= (2.0 / (rep->_nsam - 1.0));
362 vF /= (pow(a, 2.0) + b);
365 = 1.0 + (rep->_nsam + 1.0) / (3.0 * (double(rep->_nsam - 1.0)));
366 uF -= 4.0 * ((rep->_nsam + 1.0) / (pow(rep->_nsam - 1.0, 2.0)))
367 * (a_n_plus1 - 2.0 * rep->_nsam / (rep->_nsam + 1.0));
371 double F = Pi - ExternalMutations;
372 F /= pow(uF * NumMut + vF * pow(NumMut, 2.0), 0.5);
377 PolySIM::FuLiDStar(
void)
const 379 if (rep->_NumPoly == 0)
380 return std::numeric_limits<double>::quiet_NaN();
381 double Singletons = double(NumSingletons());
382 double NumMut = double(NumMutations());
384 double a = a_sub_n();
385 double b = b_sub_n();
386 double d = d_sub_n();
388 double vD = pow(rep->_nsam / (rep->_nsam - 1.0), 2.0) * b;
389 vD += pow(a, 2.0) * d;
390 vD -= 2.0 * (rep->_nsam * a * (a + 1.0))
391 / (pow(
double(rep->_nsam - 1.0), 2.0));
392 vD /= (pow(a, 2.0) + b);
394 double uD = (rep->_nsam / (rep->_nsam - 1.0))
395 * (a - (rep->_nsam / (rep->_nsam - 1.0)))
398 double DStar = (rep->_nsam / (rep->_nsam - 1.0)) * NumMut
399 - a *
double(Singletons);
400 DStar /= pow(uD * NumMut + vD * pow(NumMut, 2.0), 0.5);
405 PolySIM::FuLiFStar(
void)
const 407 if (rep->_NumPoly == 0)
408 return std::numeric_limits<double>::quiet_NaN();
409 double Singletons = double(NumSingletons());
410 double Pi = ThetaPi();
411 double NumMut = double(NumMutations());
413 double a = a_sub_n();
414 double a_n_plus1 = a_sub_n_plus1();
415 double b = b_sub_n();
418 double vF = 2.0 * pow(rep->_nsam, 3.0) + 110.0 * pow(rep->_nsam, 2.0)
419 - 255.0 * rep->_nsam + 153.0;
420 vF /= (9.0 * pow(rep->_nsam, 2.0) * (rep->_nsam - 1.0));
421 vF += (((2.0 * (rep->_nsam - 1.0) * a) / pow(rep->_nsam, 2.0))
422 - (8.0 * b / rep->_nsam));
423 vF /= (pow(a, 2.0) + b);
425 double uF = (4.0 * pow(rep->_nsam, 2.0) + 19.0 * rep->_nsam + 3.0
426 - 12.0 * (rep->_nsam + 1.0) * a_n_plus1);
427 uF /= (3.0 * rep->_nsam * (rep->_nsam - 1.0));
431 = Pi - (((rep->_nsam - 1.0) / rep->_nsam)) * double(Singletons);
432 FStar /= pow((uF * NumMut + vF * pow(NumMut, 2.0)), 0.5);
437 PolySIM::NumMutations(
void)
const 442 return rep->_NumPoly;
448 PolySIM::NumSingletons(
void)
const 456 if (rep->_counted_singletons ==
false)
458 unsigned singletonCount = 0;
460 for (std::vector<stateCounter>::const_iterator i
461 = rep->_counts.begin();
462 i != rep->_counts.end(); ++i)
465 += (i->one == 1 || i->zero == 1) ? 1u : 0u;
467 rep->_counted_singletons =
true;
468 rep->_singletons = singletonCount;
469 return (singletonCount);
472 return rep->_singletons;
474 return rep->_singletons;
478 PolySIM::NumExternalMutations(
void)
const 489 unsigned numExternal = 0;
490 for (std::vector<stateCounter>::const_iterator i
491 = rep->_counts.begin();
492 i != rep->_counts.end(); ++i)
499 return (numExternal);
503 PolySIM::Minrec(
void)
const 511 if (rep->_NumPoly < 2)
513 unsigned a, b, e, numgametes, Rmin = 0, x = 0;
516 if (rep->_NumPoly < 2 || x >= (rep->_NumPoly - 1))
518 for (a = x + 1; a < rep->_nsites; ++a)
520 for (b = (flag ==
false) ? x : a - 1; b < a; ++b)
524 for (e = 0; e < rep->_nsam; ++e)
526 if ((*rep->_data)[e][b] ==
'0' 527 && (*rep->_data)[e][a] ==
'0')
533 for (e = 0; e < rep->_nsam; ++e)
535 if ((*rep->_data)[e][b] ==
'0' 536 && (*rep->_data)[e][a] ==
'1')
542 for (e = 0; e < rep->_nsam; ++e)
544 if ((*rep->_data)[e][b] ==
'1' 545 && (*rep->_data)[e][a] ==
'0')
551 for (e = 0; e < rep->_nsam; ++e)
553 if ((*rep->_data)[e][b] ==
'1' 554 && (*rep->_data)[e][a] ==
'1')
576 PolySIM::WallStats(
void)
const 578 std::lock_guard<std::mutex> lock(rep->instance_lock);
579 if (!rep->_calculated_wall_stats)
581 unsigned S = rep->_NumPoly;
584 unsigned n00, n01, n10, n11;
585 unsigned nhap_curr, nhap_left;
586 unsigned n0site1, n0site2;
590 for (
unsigned site1 = 0; site1 < rep->_nsites - 1;
595 for (
unsigned site2 = site1 + 1;
596 site2 < rep->_nsites; ++site2)
599 n00 = n01 = n10 = n11 = n0site1
601 for (
unsigned i = 0; i < rep->_nsam;
605 = (*rep->_data)[i][site2];
607 (*rep->_data)[i][site1])
639 && n0site1 < rep->_nsam)
641 && n0site2 < rep->_nsam))
643 nhap_curr += (n00 > 0) ? 1 : 0;
644 nhap_curr += (n01 > 0) ? 1 : 0;
645 nhap_curr += (n10 > 0) ? 1 : 0;
646 nhap_curr += (n11 > 0) ? 1 : 0;
651 ++rep->_walls_Bprime;
658 ++rep->_walls_Bprime;
665 nhap_left = nhap_curr;
669 = double(rep->_walls_Bprime) / (double(S - 1));
671 = (rep->_walls_Bprime + double(A)) / (
double(S));
676 = std::numeric_limits<double>::quiet_NaN();
677 rep->_walls_Bprime = 0;
679 = std::numeric_limits<double>::quiet_NaN();
682 rep->_calculated_wall_stats =
true;
686 PolySIM::WallsB(
void)
const 689 return rep->_walls_B;
692 PolySIM::WallsBprime(
void)
const 695 return rep->_walls_Bprime;
698 PolySIM::WallsQ(
void)
const 701 return rep->_walls_Q;
double thetal(const AlleleCountMatrix &ac, const std::int8_t refstate)
Zeng et al. .
Methods dealing with recombination.
The namespace in which this library resides.
const unsigned SEQMAXUNSIGNED
declaration of Sequence::PolySIM, a class to analyze coalescent simulation data
namespace Sequence::Recombination
Declaration of Sequence::SimData, a class representing polymorphism data from coalescent simulations ...
Data from coalescent simulations.