This post covers how to run pairwise f2-statistics and FST in AdmixPy. They are simple to interpret, and are also useful computationally. Once f2 blocks have been computed and cached, many downstream analyses can reuse them without repeatedly reading and converting the original genotype data.

What does f2 measure?

The f2-statistic quantifies allele-frequency differentiation between two populations, AA and BB, and is defined as:

f2(A,B)=E[(pA−pB)2] f_2(A, B) = E[(p_A - p_B)^2]

where pAp_A and pBp_B are the allele frequencies of populations AA and BB at a SNP, and the squared allele-frequency differences are averaged across SNPs.

This makes f2 a compact, distance-like summary of population differentiation. Because of sampling variation, bias-corrected estimates can sometimes be slightly negative.

Bias-corrected f2

A basic calculation of f2 would simply average the squared allele-frequency differences:

(pA−pB)2 (p_A - p_B)^2

But allele frequencies are estimated from finite samples. If a population has only one or a few observed alleles at a SNP, the observed allele frequency can be different from the true population frequency because of sampling noise. Without correction, this sampling noise inflates f2.

AdmixPy uses an allele-count-aware finite-sample correction:

f^2(A,B)=E[(pA−pB)2−pA(1−pA)cA−1−pB(1−pB)cB−1], \hat f_2(A,B)=E\left[ (p_A-p_B)^2 -\frac{p_A(1-p_A)}{c_A-1} -\frac{p_B(1-p_B)}{c_B-1} \right],

Here, cAc_A and cBc_B are the observed allele counts for populations AA and BB. With the bias correction default, SNP values with cA≤1c_A\leq1 or cB≤1c_B\leq1 are excluded.

With bias correction enabled, AdmixPy requires at least two observed alleles in each population at a SNP. Values with (c<2) are excluded with a warning. Setting apply_corr=False retains them, but explicitly returns a sampling-biased estimate.

Diploid and pseudohaploid samples

Finite-sample correction is based on the number of observed alleles. Ploidy therefore directly affects estimates for populations represented by a single sample or only a few samples.

A diploid individual contributes two alleles per SNP:

c=2 c = 2

Its observed allele frequency can be:

p∈{0,0.5,1} p \in \{0, 0.5, 1\}

A pseudohaploid individual contributes only one allele per SNP:

c=1 c = 1

Its observed allele frequency can only be:

p∈{0,1} p \in \{0, 1\}

For larger population groups, the effect is reduced because the allele-frequency estimate is averaged across more samples. For example, a pseudohaploid population represented by many individuals still has many observed alleles in total, even though each individual contributes only one allele per SNP. In that case, pseudohaploid calling is still taken into account.

f2 and FST are population-level summaries. Bias-corrected f2 and FST can be calculated for a single diploid sample, although the estimates may be noisy. A population represented by one pseudohaploid sample has only one observed allele per called SNP, so bias-corrected estimates cannot be calculated for those SNPs.

Running pairwise f2 in AdmixPy

For AdmixPy setup instructions, refer to the previous post: Introducing AdmixPy: f-statistics, qpAdm, and qpWave in Python. After activating the virtual environment, start a Python REPL by entering python in the terminal.

Then import AdmixPy:

import admixpy as ap

Since we imported admixpy as ap, pairwise f2-statistics can now be run with ap.f2().

For this example, I will compare a Late Neolithic population from Shandong against a set of modern East Asian populations:

# Adjust to your dataset prefix or path
prefix = "ho_v66"

pops = ["Han", "Japanese", "Korean", "Mongol", "Tibetan"]

Now run:

ap.f2(prefix, pop1="China_Shandong_Dinggong_LN", pop2=pops)

The result:

                         pop1      pop2         est           se       n
0  China_Shandong_Dinggong_LN       Han  0.00220414  0.000125263  206614
1  China_Shandong_Dinggong_LN  Japanese  0.00413708  0.000136411  206614
2  China_Shandong_Dinggong_LN    Korean  0.00460585  0.000190963  206614
3  China_Shandong_Dinggong_LN    Mongol  0.00462775  0.000142786  206614
4  China_Shandong_Dinggong_LN   Tibetan  0.00379421  0.000132498  206614

est is calculated by weighting each block estimate by the number of usable SNPs for that population pair. n is the total number of contributing SNPs, and se is the block-jackknife standard error, calculated by leaving out one block at a time.

The interpretation is simple: lower f2 means the two populations are more similar in allele-frequency space, while higher f2 means they are more differentiated. In this example, the Han average has the lowest f2 value with China_Shandong_Dinggong_LN, suggesting that among the listed comparison populations, Han are closest to this Shandong Late Neolithic group in allele-frequency space.

This should not be interpreted as an ancestry model. The result only says that, among this set of populations, Han show the lowest pairwise allele-frequency differentiation from China_Shandong_Dinggong_LN.

Because China_Shandong_Dinggong_LN contains multiple samples, it provides enough observed alleles for bias correction at more SNPs than a single pseudohaploid sample population would.

Effect of the finite-sample correction

For comparison, the same analysis can be run without the finite-sample correction (apply_corr=False):

                         pop1      pop2        est           se       n
0  China_Shandong_Dinggong_LN       Han  0.0226191  0.000127416  206955
1  China_Shandong_Dinggong_LN  Japanese  0.0251914  0.000140491  206955
2  China_Shandong_Dinggong_LN    Korean  0.0334269  0.000201275  206955
3  China_Shandong_Dinggong_LN    Mongol  0.0255762  0.000140268  206955
4  China_Shandong_Dinggong_LN   Tibetan  0.0244494  0.000134095  206955

Without the correction, the estimates are approximately 5.5 to 10 times higher. Han still has the lowest f2, but Korean and Mongol switch ranks.

Blocks and standard errors

Like other ADMIXTOOLS-style statistics, AdmixPy estimates uncertainty using SNP blocks. SNPs are grouped by chromosome and genetic position. AdmixPy computes f2 separately for each block, then combines the blocks into the final estimate.

The standard error is calculated using a leave-one-block-out jackknife. In each jackknife replicate, one block is left out and the statistic is recomputed from the remaining blocks. The variation across these leave-one-block-out estimates gives the standard error.

If standard errors were computed as if every SNP were independent, they would usually be too small. Block jackknifing gives a more realistic estimate of uncertainty by accounting for correlation among nearby SNPs.

Using get_f2() and the f2 cache

The main helper is get_f2(). It accepts genotype data, an in-memory F2Blocks object, or an on-disk f2 cache.

For example, you can compute f2 blocks from a TGENO prefix:

blocks = ap.get_f2(prefix, pops=["Chimp", "Sardinian", "Orcadian", "Norwegian", "Mongol"])

AdmixPy automatically detects pseudohaploid samples by default and adjusts their allele counts when computing allele frequencies. In most cases this does not need to be changed. It can be changed by adding the adjust_pseudohaploid parameter; False treats all samples as diploid, while an integer changes how many SNPs are checked during pseudohaploid detection.

If you want to reuse these blocks later, you can write them to disk with write_f2():

ap.write_f2(blocks, "f2_cache")

This creates a folder called f2_cache. Inside it, the block-level pairwise estimates, per-pair usable-SNP counts, and block lengths are saved. Then you can load the folder again with get_f2():

blocks = ap.get_f2("f2_cache")

The loaded blocks object can be passed to functions like ap.f2(), ap.f3(), or ap.f4().

Cached f2 blocks can be reused for f3, f4, qpWave, and qpAdm, but they cannot reproduce direct-genotype allsnps=True calculations when different statistics require different SNP intersections.

The default maxmiss=0 keeps only SNPs observed in every selected population, which can be restrictive for groups with high missingess. With finite-sample correction, SNPs with fewer than two allele observations are then excluded pairwise. Setting maxmiss=1 allows pairs to use their available SNP overlap, improving estimates for high-missingness populations. remove_na=False preserves blocks lacking data for some pairs instead of discarding them globally.

With many aDNA populations, setting remove_na=False - Trueis the default - is often needed for building a usable cache. When coverage varies across populations, blocks can lack estimates for some population pairs. Keeping those blocks nevertheless preserves the available pairwise data instead of discarding the entire block.

Running FST in AdmixPy

AdmixPy also provides pairwise Hudson-style FST through ap.fst():

ap.fst(prefix, pop1="China_Shandong_Dinggong_LN", pop2=pops)

Result:

                         pop1      pop2         est           se       n
0  China_Shandong_Dinggong_LN       Han  0.00863094  0.000492654  273728
1  China_Shandong_Dinggong_LN  Japanese   0.0162544  0.000530885  273728
2  China_Shandong_Dinggong_LN    Korean    0.017967  0.000742235  273728
3  China_Shandong_Dinggong_LN    Mongol   0.0176329  0.000540258  273728
4  China_Shandong_Dinggong_LN   Tibetan   0.0147032  0.000511363  273728

For each block (b), AdmixPy calculates Hudson-style FST as:

FST,b=∑s∈b[(pA−pB)2−pA(1−pA)cA−1−pB(1−pB)cB−1]∑s∈b[(pA−pB)2+pA(1−pA)+pB(1−pB)] F_{ST,b}= \frac{ \sum_{s\in b}\left[ (p_A-p_B)^2- \frac{p_A(1-p_A)}{c_A-1}- \frac{p_B(1-p_B)}{c_B-1} \right] }{ \sum_{s\in b}\left[ (p_A-p_B)^2+p_A(1-p_A)+p_B(1-p_B) \right] }

With the bias corrected default, SNPs with cA≤1c_A\leq1 or cB≤1c_B\leq1 are excluded. Setting apply_corr=False removes the two correction terms.

By default, fst_aggregation="block_ratios" averages these block-level estimates using their usable SNP counts. fst_aggregation="pooled_components" instead pools the numerators and denominators across blocks before taking the ratio. The results can differ when genetic variation varies among blocks.

FST is closely related to f2, but it is normalized by the amount of genetic variation in the compared populations. AdmixPy uses a Hudson-style estimator: the numerator is the same finite-sample-corrected allele-frequency difference used for f2, while the denominator is the expected pairwise allele difference between the two populations.

By default, f2 uses only sites polymorphic across populations, while FST keeps them. To compare f2 and FST on the same sites, pass poly_only=True to both calls or poly_only=False.

f2 and FST often give a similar ordering of comparison populations. In this example, Han again has the lowest value relative to China_Shandong_Dinggong_LN among the listed populations. The difference is scale and interpretation: f2 is an unnormalized measure of allele-frequency differentiation, while FST asks how large that differentiation is relative to the total allele-frequency variation available at the SNPs being compared.

Summary

At the pairwise level, f2 and FST summarize how differentiated two populations are in allele-frequency space. Lower estimates generally mean more similar allele frequencies; higher estimates mean more differentiation.

AdmixPy uses a finite-sample correction based on observed allele counts, which is relevant for single-sample populations and pseudohaploid ancient DNA. It also computes statistics by genomic blocks and uses a block jackknife for standard errors.

A practical advantage is caching. By storing block-level f2 values in an F2Blocks object or an on-disk f2 cache, AdmixPy makes repeated f2, f3, f4, qpWave, and qpAdm analyses much faster. Instead of repeatedly reading and converting genotype data, you can compute the shared foundation once and reuse it across many tests.

For exploratory population-genetic analysis, this makes f2 both a useful statistic and a useful workflow primitive.