From 4343155a017a8bc502af3b84cea966077ac4dde0 Mon Sep 17 00:00:00 2001 From: Steve Mussmann Date: Mon, 30 Jul 2018 16:07:05 -0600 Subject: [PATCH] added f2, f3, f4, and fst calculations per locus --- fnStats.cpp | 156 ++++++++++++++++++++++++++++++++++++++++++++++++---- fnStats.h | 6 ++ fn_main.cpp | 1 + 3 files changed, 151 insertions(+), 12 deletions(-) diff --git a/fnStats.cpp b/fnStats.cpp index e51860e..0b1df0c 100644 --- a/fnStats.cpp +++ b/fnStats.cpp @@ -4,6 +4,7 @@ //#include //#include //#include +#include #include #include //#include @@ -101,9 +102,9 @@ void fnStats::calcAFreqs(fnFiles &f) } it++; } - std::cout << i << std::endl; - std::cout << "anc = " << Afreqs[i]["anc"] << std::endl; - std::cout << "der = " << Afreqs[i]["der"] << std::endl; + //std::cout << i << std::endl; + //std::cout << "anc = " << Afreqs[i]["anc"] << std::endl; + //std::cout << "der = " << Afreqs[i]["der"] << std::endl; } } @@ -149,9 +150,9 @@ void fnStats::calcBFreqs(fnFiles &f) } it++; } - std::cout << i << std::endl; - std::cout << "anc = " << Bfreqs[i]["anc"] << std::endl; - std::cout << "der = " << Bfreqs[i]["der"] << std::endl; + //std::cout << i << std::endl; + //std::cout << "anc = " << Bfreqs[i]["anc"] << std::endl; + //std::cout << "der = " << Bfreqs[i]["der"] << std::endl; } } @@ -197,9 +198,9 @@ void fnStats::calcCFreqs(fnFiles &f) } it++; } - std::cout << i << std::endl; - std::cout << "anc = " << Cfreqs[i]["anc"] << std::endl; - std::cout << "der = " << Cfreqs[i]["der"] << std::endl; + //std::cout << i << std::endl; + //std::cout << "anc = " << Cfreqs[i]["anc"] << std::endl; + //std::cout << "der = " << Cfreqs[i]["der"] << std::endl; } } @@ -245,8 +246,139 @@ void fnStats::calcDFreqs(fnFiles &f) } it++; } - std::cout << i << std::endl; - std::cout << "anc = " << Dfreqs[i]["anc"] << std::endl; - std::cout << "der = " << Dfreqs[i]["der"] << std::endl; + //std::cout << i << std::endl; + //std::cout << "anc = " << Dfreqs[i]["anc"] << std::endl; + //std::cout << "der = " << Dfreqs[i]["der"] << std::endl; + } +} + +void fnStats::calcFstats(fnFiles &f) +{ + for(unsigned int i=0; i al = f.getALocus(i); + std::unordered_map bl = f.getBLocus(i); + std::unordered_map cl = f.getCLocus(i); + std::unordered_map dl = f.getOutgroupLocus(i); + + //initialize iterators + std::unordered_map::iterator ait = al.begin(); + std::unordered_map::iterator bit = bl.begin(); + std::unordered_map::iterator cit = cl.begin(); + std::unordered_map::iterator dit = dl.begin(); + + std::vector avec; + std::vector bvec; + std::vector cvec; + std::vector dvec; + + int atotal=0; + while(ait != al.end()) + { + atotal += ait->second; + avec.push_back(ait->second); + ait++; + } + + int btotal=0; + while(bit != bl.end()) + { + btotal += bit->second; + bvec.push_back(bit->second); + bit++; + } + + int ctotal=0; + while(cit != cl.end()) + { + ctotal += cit->second; + cvec.push_back(cit->second); + cit++; + } + + int dtotal=0; + while(dit != dl.end()) + { + dtotal += dit->second; + dvec.push_back(dit->second); + dit++; + } + + double ahz = hz(avec, atotal); + double bhz = hz(bvec, btotal); + double chz = hz(cvec, ctotal); + double dhz = hz(dvec, dtotal); + + double f2i = f2(Afreqs[i]["anc"], Bfreqs[i]["anc"], ahz, bhz, atotal, btotal); + double f3i = f3(Afreqs[i]["anc"], Bfreqs[i]["anc"], Cfreqs[i]["anc"], chz, ctotal); + double f4i = f4(Afreqs[i]["anc"], Bfreqs[i]["anc"], Cfreqs[i]["anc"], Dfreqs[i]["anc"]); + double fsti = fst(Afreqs[i]["der"], Bfreqs[i]["der"]); + + //std::cout << std::fixed; + std::cout << "len(a) = " << avec.size() << std::endl; + std::cout << "len(b) = " << bvec.size() << std::endl; + std::cout << "len(c) = " << cvec.size() << std::endl; + std::cout << "ahz = " << ahz << std::endl; + std::cout << "bhz = " << bhz << std::endl; + std::cout << "chz = " << chz << std::endl; + std::cout << "f2 = " << f2i << std::endl; + std::cout << "f3 = " << f3i << std::endl; + //std::cout << std::setprecision(10) << "f4 = " << f4i << std::endl; + std::cout << "f4 = " << f4i << std::endl; + std::cout << "fst = " << fsti << std::endl; + std::cout << std::endl; + + } +} + +double fnStats::hz(std::vector &v, int t) +{ + if(v.size() == 1){ + return 0.0; + } + else + { + return ((double)v[0]*(double)v[1])/((double)t*((double)t-1.0)); + } +} + + +double fnStats::f2(double freqa, double freqb, double ha, double hb, int tota, int totb) +{ + double x = freqa-freqb; + return (x*x)-(ha/(double)tota)-(hb/(double)totb); +} + +double fnStats::f3(double freqa, double freqb, double freqc, double hc, int totc) +{ + double x = freqc-freqa; + double y = freqc-freqb; + double xy = (x*y); + return xy-(hc/(double)totc); +} + +double fnStats::f4(double freqa, double freqb, double freqc, double freqd) +{ + double x = freqa-freqb; + double y = freqc-freqd; + return x*y; +} + +double fnStats::fst(double freqa, double freqb) +{ + double n = (freqa-freqb); + n = n*n; + double b = 1.0-freqb; + double a = 1.0-freqa; + b = freqa*b; + a = freqb*a; + if(a+b == 0) + { + return 0.0; + } + else + { + return n/(a+b); } } diff --git a/fnStats.h b/fnStats.h index 90e2ec4..bac5816 100644 --- a/fnStats.h +++ b/fnStats.h @@ -14,6 +14,7 @@ class fnStats { fnStats(int l); void calcAllFreqs(fnFiles &f); void findAncestral(fnFiles &f); + void calcFstats(fnFiles &f); private: //record alleles as anc (ancestral) and derived (der) std::vector > Afreqs; @@ -26,6 +27,11 @@ class fnStats { void calcBFreqs(fnFiles &f); void calcCFreqs(fnFiles &f); void calcDFreqs(fnFiles &f); + double hz(std::vector &v, int t); + double f2(double freqa, double freqb, double ha, double hb, int tota, int totb); + double f3(double freqa, double freqb, double freqc, double hc, int totc); + double f4(double freqa, double freqb, double freqc, double freqd); + double fst(double freqa, double freqb); }; #endif diff --git a/fn_main.cpp b/fn_main.cpp index ba978d3..8fe17f2 100644 --- a/fn_main.cpp +++ b/fn_main.cpp @@ -28,6 +28,7 @@ int main(int argc, char** argv) { fnStats fs(f.getLength()); fs.findAncestral(f); fs.calcAllFreqs(f); + fs.calcFstats(f); return 0; }