Mercurial > repos > miller-lab > genome_diversity
comparison genome_diversity/src/Fst_lib.c @ 14:8ae67e9fb6ff
Uploaded Miller Lab Devshed version a51c894f5bed again [possible toolshed.g2 bug]
author | miller-lab |
---|---|
date | Fri, 28 Sep 2012 11:35:56 -0400 |
parents | 2c498d40ecde |
children | f04f40a36cc8 |
comparison
equal
deleted
inserted
replaced
13:fdb4240fb565 | 14:8ae67e9fb6ff |
---|---|
1 // procedure to compute either Wright's Fst or an unbiased estimator of if | |
2 | |
3 #include "lib.h" | |
4 // Wright's Fst | |
5 static double Wright(double f1, double f2) { | |
6 double | |
7 f, // frequency in the pooled population | |
8 H_ave, // average of HWE heterogosity in the two populations | |
9 H_all; // HWE heterozygosity in the pooled popuations | |
10 | |
11 H_ave = f1*(1.0 - f1) + f2*(1.0 - f2); | |
12 f = (f1 + f2)/2.0; | |
13 if (f == 0.0 || f == 1.0) | |
14 return 0.0; | |
15 H_all = 2.0*f*(1.0 - f); | |
16 return (H_all - H_ave) / H_all; | |
17 } | |
18 | |
19 /* unbiased estimator of Fst from: | |
20 Weir, B.S. and Cockerham, C.C. 1984. Estimating F-statistics for the | |
21 analysis of population structure. Evolution 38: 1358–1370. | |
22 as interpreted by: | |
23 Akey, J.M., Zhang, G., Zhang, K., Jin, L., and Shriver, M.D. 2002. | |
24 Interrogating a high-density SNP map for signatures of natural | |
25 selection. Genome Res. 12: 1805–1814. | |
26 */ | |
27 static double Weir(int n1, double p1, int n2, double p2) { | |
28 double F, p_bar, nc, MSP, MSG, N = n1 + n2; | |
29 | |
30 if (p1 == p2) | |
31 return 0.0; | |
32 MSG = (n1*p1*(1.0-p1) + n2*p2*(1.0-p2))/(N-1.0); | |
33 p_bar = (n1*p1 + n2*p2)/N; | |
34 MSP = n1*(p1-p_bar)*(p1-p_bar) + n2*(p2-p_bar)*(p2-p_bar); | |
35 nc = N - (double)(n1*n1 + n2*n2)/N; | |
36 F = (MSP - MSG) / (MSP + (nc-1)*MSG); | |
37 if (F < 0.0) | |
38 F = 0.0; | |
39 return F; | |
40 } | |
41 | |
42 double Fst(int nA1, int na1, int nA2, int na2, int unbiased) { | |
43 double p1, p2; | |
44 | |
45 p1 = (double)nA1 / (double)(nA1+na1); | |
46 p2 = (double)nA2 / (double)(nA2+na2); | |
47 | |
48 return (unbiased ? Weir(nA1+na1, p1, nA2+na2, p2) : Wright(p1, p2)); | |
49 } |