-
Notifications
You must be signed in to change notification settings - Fork 2
/
FstAlleleCounter.cpp
96 lines (76 loc) · 2.56 KB
/
FstAlleleCounter.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
/*
* FstAlleleCounter
* Date: Aug-15-2012
* Author : Gabriel Renaud [email protected]
*
*/
#include "FstAlleleCounter.h"
FstAlleleCounter::FstAlleleCounter(){
reinitializedCounters();
}
void FstAlleleCounter::reinitializedCounters(){
counterSame =0;
counterCommon =0;
counterReference =0;
counterSample =0;
}
double FstAlleleCounter::lowerConf (const unsigned int shortBranch,const unsigned int commonBranch) const{
double pHat = double(shortBranch)/ ( double(shortBranch)+double(commonBranch));
//1.0 = 85%,
//1.6 = 95%
double z = 1.6;
double n=shortBranch+commonBranch;
//taken from http://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval
double numerator= ( pHat + (z*z/(2.0*n)) - z* sqrt( (pHat*(1.0-pHat)+(z*z/(4.0*n)) )/n ));
double denominator = (1.0+(z*z/n));
// cout<<"pHat "<<pHat<<endl;
// // cout<<(z*z)<<endl;
// // cout<<(2.0*n)<<endl;
// cout<<((z*z)/(2.0*n))<<endl;
// cout<<numerator<<endl;
// cout<<denominator<<endl;
// cout<<"####"<<endl;
return ( numerator / denominator );
}
string FstAlleleCounter::getHeader(string prefixToAdd){
return
prefixToAdd+"noMut\t"+
prefixToAdd+"common\t"+
prefixToAdd+"ind1Spec\t"+
prefixToAdd+"ind2Spec\t"+
prefixToAdd+"divInd1M\t"+
prefixToAdd+"divInd1L\t"+
prefixToAdd+"divInd1H\t"+
prefixToAdd+"divInd2M\t"+
prefixToAdd+"divInd2L\t"+
prefixToAdd+"divInd2H";
}
double FstAlleleCounter::highConf (const unsigned int shortBranch,const unsigned int commonBranch) const {
double pHat = double(shortBranch)/ ( double(shortBranch)+double(commonBranch));
//1.0 = 85%,
//1.6 = 95%
double z = 1.6;
double n=shortBranch+commonBranch;
//taken from http://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval
double numerator= ( pHat + (z*z/(2.0*n)) + z* sqrt( (pHat*(1.0-pHat)+(z*z/(4.0*n)) )/n ));
double denominator = (1.0+(z*z/n));
return ( numerator / denominator );
}
// int main (int argc, char *argv[]) {
// FstAlleleCounter al;
// cout<<al.lowerConf(46,396)<<endl;
// }
double FstAlleleCounter::fstRefSam () const {
if( (counterReference + counterCommon) != 0){
return double(counterReference)/double(counterReference+counterCommon);
}else{
return std::numeric_limits<double>::infinity();
}
}
double FstAlleleCounter::fstSamRef () const {
if( (counterSample + counterCommon) != 0){
return double(counterSample)/double(counterSample+counterCommon);
}else{
return std::numeric_limits<double>::infinity();
}
}