-
Notifications
You must be signed in to change notification settings - Fork 0
/
drift.py
56 lines (43 loc) · 1.59 KB
/
drift.py
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
# Functions for doing drift-related calculations.
from __future__ import division, print_function
import numpy as np
import pdb
###########################################################################
def estimate_covariance(freq_data, T):
"""
Using the notation of the Berg and Coop paper
"""
G=np.array(freq_data, dtype=float)
G=G.T
M,K = G.shape
eps=np.mean(G, axis=0)
var2=np.expand_dims(1/eps*(1-eps),axis=0)
TGs=T.dot(G*var2)
F=TGs.dot(TGs.T)/(K-1)
return F
###########################################################################
def Fst_matrix(freq, pops):
"""
Method by Nick Patterson. Should be equivalent to the inbreed=True
Option in Eigenstrat/smartpca
"""
npops=len(pops)
Fst=np.zeros( (npops, npops), dtype="float")
for i in range(npops-1):
for k in range(i+1, npops):
s=freq["total"][:,i]
t=freq["total"][:,k]
u=freq["count"][:,i]
v=freq["count"][:,k]
EX=np.mean(np.square(np.true_divide(u,s)-np.true_divide(v,t)))
Eh1=np.mean(np.true_divide(u*(s-u), s*(s-1)))
Eh2=np.mean(np.true_divide(v*(t-v), t*(t-1)))
Eh1s=np.mean(np.true_divide(u*(s-u), s*s*(s-1)))
Eh2t=np.mean(np.true_divide(v*(t-v), t*t*(t-1)))
Nhat=EX-Eh1s-Eh2t
Dhat=Nhat+Eh1+Eh2
# if Nhat/Dhat<0:
# pdb.set_trace()
Fst[i,k]=Fst[k,i]=Nhat/Dhat
return {"Fst":Fst, "pops":pops}
###########################################################################