Open atongsa opened 1 year ago
A pull request would be very welcome if you'd like to contribute @atongsa! If you're interested, we'd be happy to provide some concrete pointers on where to start. (The development docs should get you started, in the first instance)
thank you, sir. I am learning the development docs and will try to add a PBS() function to tskit.
i imagine the ts.PBS() logic is like,
1) read vcf and 3 pops
2) use the ts.fst() to get 3 pair fst and then calculate the PBS for all 3 branch
but i dont find ways to read vcf into tree sequence (just write_vcf) in tskit document, can you give me some prompt
As it's a "derived" statistic, I would imagine the implementation would look something like that of Fst. However, it may be that it's not as complicated as that.
Perhaps you could try implementing this as a standalone function based on some simple simulations, and paste your code in here?
We wouldn't be working with VCF here, just tskit tree sequence objects.
but i dont find ways to read vcf into tree sequence (just write_vcf) in tskit document, can you give me some prompt
Just FYI @atongsa, there is a long and complicated "inference" procedure to try and construct a (sometimes) reasonable tree sequence from the much more limited genetic variation data present in a VCF (e.g. via tsinfer). Moreover, even once you manage to infer a tree sequence, to calculate branch length stats like Fst, you need to date the tree sequence too (e.g. via tsdate).
So to implement stats calculations in tskit, you'll want to work with directly simulated tree sequences, not VCF data.
i write a demo function, and i will improve it later for sites and windows calculation (because i don't understand some function in trees.py, i need some time to learn), and then try to make a tskit.pbs() pull request
def pbs(ts, test, ref1, ref2):
"""
Calculates the population branch statistic (PBS) between a test population and two reference populations.
Parameters
----------
ts: tree sequence
test : numpy.ndarray, nodes as the focus population
ref1 : numpy.ndarray, nodes as the far away population
ref2 : numpy.ndarray, nodes as the sister group of focus population
Returns
-------
PBS values of the focus population
"""
# Calculate pairwise Fst values
fst_test_ref1 = ts.Fst([test, ref1])
fst_test_ref2 = ts.Fst([test, ref2])
fst_ref1_ref2 = ts.Fst([ref1, ref2])
# Calculate PBS
pbs = (-math.log(1-fst_test_ref1, 10) + -math.log(1-fst_test_ref2, 10) - -math.log(1-fst_ref1_ref2, 10)) / 2
return pbs
the test code below get a return value 0.1198192924722683/2
import math
import tskit
import numpy as np
ts = msprime.sim_ancestry(3, recombination_rate=0.1, sequence_length=10, random_seed=42)
mts = msprime.sim_mutations(ts, rate=0.1, random_seed=42)
p1 = mts.samples()[[0,5]]
p2 = mts.samples()[[2,3]]
p3 = mts.samples()[[1,4]]
pbs(mts,p1,p2,p3)
For reference, here's the definition of PBS:
The code looks right, although I think you're missing a factor of 2? And, it could be slightly more concise if you used the indexes
argument to Fst
.
(It's also interesting to note that PBS is framed almost in terms of just a branch length; it's worth thinking about what a tree-based analogous thing would be, although I'm not suggesting that for this issue.)
Working through the motivation and implementation in #2777 here: so, the statistic is defined in terms of some variables called T
; these are -log(1-Fst)
, which are "divergence time in units scaled by the population size"; later on in the document it says that the log transform "places branches of different magnitudes on the same scale". In the paper they say they computed Fst using the method of Reynolds; OF COURSE, Reynolds et al do not actually use the notation "Fst" anywhere. However, they discuss estimators of the 'coancestry coefficient', which is defined to be the probability of identity between two alleles drawn from the same population, and provide a ratio of variances estimate for this; this is (one of the usual definitions of) Fst. And, they have this:
The main question I have at this point is whether we should really be using T = -log(1-Fst)
in our definition of PBS or if we should be using T = divergence
.
Now, our Fst is implemented as
Fst = 1 - 2 * (d(X) + d(Y)) / (d(X) + 2 * d(X, Y) + d(Y))
and so letting d(X+Y) = (d(X) + 2 * d(X, Y) + d(Y))
,
-log(1 - Fst) = -log(2 * (d(X) + d(Y)) / d(X+Y))
= log(d(X+Y)) - log(d(X) + d(Y)) - log(2)
and so plugging this in to the definition of PBS we'd get that
PBS(T,H,D) = (
log(d(T+H)) - log(d(T) + d(H)) - log(2)
+ log(d(T+D)) - log(d(T) + d(D)) - log(2)
- log(d(H+D)) + log(d(H) + d(D)) + log(2)
) / 2
= log(
( d(T+H) * d(T+D) * (d(H) + d(D)) )
/ ( 2 * (d(T) + d(H)) * (d(T) + d(D)) * d(H+D) )
) / 2
I'm not seeing any particular simplification or insight here, just writing it down.
However, notice that since (as defined here) Fst is a unitless quantity, PBS is as well, even though it's defined with some variables named T
which suggest "time". In the last expression we see this because it's the log of a ratio of three products of times. This is quite different to what we'd get if we used divergence directly for T
.
Okay - I think defining PBS in terms of Fst is what we ought to do - at least, it agrees with the literature. It would be interesting to look at the other statistic we'd get by just putting in divergences (I think that's been called a "y-statistic"?), but that's a separate topic.
PBS (population branch statistics) [^1] is derived from Fst, which can detect positive selection signals in admixted populations .
the formula is
[^1]: Yi, X. et al. Sequencing of 50 Human Exomes Reveals Adaptation to High Altitude. Science 329, 75–78 (2010).