Fst in tskit (equivalent to Reynolds et al 1983)

[mshpak76]

Does anybody know if the Fst implemented in tskit is formally equivalent to the coefficient of coancestry derived in Reynolds et al 1983 (theta on pg 769, Genetics 105:767-779)? There are different ways of estimating Fst, but when I collect the various x1,x2,n1,n2 terms in

Fst = 1 - 2 * (d(X) + d(Y)) / (d(X) + 2 * d(X, Y) + d(Y))

(see summary functions at the end of https://github.com/tskit-dev/tskit/discussions/new ), I don't quite get the Reynolds et al expression. I don't know if something went wrong with my algebra or if the estimates are in fact not equivalent.

[petrelharp]

I think it is not formally equivalent, even for the biallelic case, as Reynolds et al is defining things in terms of variance components (that's what a and b are), and our definition is in terms of probabilities of idenity, but I'd have to do a fair bit of math to check. We certainly could add the Reynolds estimator as an option (e.g., ts.Fst(..., method="reynolds"), but for genomic data (ie, lots of mostly biallelic snps) the difference is negligible (at least in cases where I've done that check).

[mshpak76]

Do you happen to have a block of python code available that implements Reynolds Fst calculation for ts objects? If not, I have the Reynolds estimator coded in a separate python script, but I'm not quite sure how I would integrate it with tskit tools and apply it to the ts object, apart from writing the haplotypes to a matrix and calculating Reynolds Fst from the frequencies.

[petrelharp]

I don't, sorry.

[mshpak76]

Unfortunately, we will probably need to run Reynolds Fst in order to compare the output of demographic models to our data for consistency. If I have a python script def rey_fst(freq1, freq2, SampleSize1, SampleSize2) that calculates Reynolds Fst at each segregating site given allele frequency and sample size for samples 1,2, what is the most efficient way to extract this information from a ts object to get window Reynolds Fst?

…

On Wed, Jun 15, 2022 at 3:11 PM Peter Ralph ***@***.***> wrote: I *think* it is not formally equivalent, even for the biallelic case, as Reynolds et al is defining things in terms of variance components (that's what a and b are), and our definition is in terms of probabilities of idenity, but I'd have to do a fair bit of math to check. We certainly *could* add the Reynolds estimator as an option (e.g., ts.Fst(..., method="reynolds"), but for genomic data (ie, lots of mostly biallelic snps) the difference is negligible (at least in cases where I've done that check). — Reply to this email directly, view it on GitHub <#2341 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AUZMZCIJ3U5MHBWGUHI4QVTVPI2HHANCNFSM5Y4AF5NA> . You are receiving this because you authored the thread.Message ID: ***@***.***>

-- ======================= Max Shpak, Ph.D. Department of Genetics University of Wisconsin Madison, WI 53706

[petrelharp]

Depends - are you wanting per-SNP Fst or the combined Fst across windows (or the genome)? If the latter, I'd get the joint allele frequency spectrum. If the former, you can get allele frequencies this way or extract them from the genotype matrix.

We would like to have this Fst definition implemented, but don't have the time to do it right now.