Interpreting f4-Statistics with AdmixPy

f4-statistics can be used to measure the covariance between allele-frequency differences across two pairs of populations. A value around zero implies that the two population pairs are symmetrically related, while a non-zero value implies excess allele sharing between one of the cross-pairs.

Theory and formula

An f4-statistic is the average, across SNPs, of the product of the allele-frequency differences between two pairs of populations:

Multiplying out results in:

In the formula, $A$, $B$, $C$, and $D$ are populations: $p_{A,i}$, $p_{B,i}$, $p_{C,i}$, and $p_{D,i}$ are their respective allele frequencies at SNP $i$, and $\mathbb{E}_i$ is the average across SNPs.

The expanded form can now also be written as shared-drift covariances:

$\operatorname{Cov}(X,Y)$ stands for the shared drift covariance between populations $X$ and $Y$.

1. General interpretation

The covariance terms can be grouped into two sets of cross-pairs:

Thus, without an outgroup assumption:

• When $f_4>0$, the combined shared-drift covariance of $(A,C)$ and $(B,D)$ is greater than that of $(A,D)$ and $(B,C)$.
• When $f_4<0$, the combined shared-drift covariance of $(A,D)$ and $(B,C)$ is greater than that of $(A,C)$ and $(B,D)$.
• When $f_4\approx0$, there is no detectable difference between the two combined pairings.

These are relative statements about the two sums of covariance terms; without further assumptions, the statistic cannot identify which individual population pair is responsible for the imbalance. For this reason, $f_4$-statistics are often used with one population specified as a suitable outgroup, as described below.

2. Interpretation with an outgroup

When $A$ is a suitable outgroup, it is expected to share approximately the same drift with $C$ and $D$, so:

These terms cancel in expectation, leaving:

Thus, the sign indicates the direction of excess allele sharing:

This also shows that an f4-statistic measures relative excess shared drift: the sign shows which pair shares more drift, and the estimate shows how large the difference is.

However, the estimate alone does not indicate whether the observed difference is statistically distinguishable from zero. This is assessed using the $Z$-score:

In this formula, the hat indicates a value estimated from the observed SNP data: $\widehat{f_4}$ is the $f_4$-estimate, and $\operatorname{SE}(\widehat{f_4})$ is its standard error. In AdmixPy, the standard error is estimated using a block jackknife. The sign of $Z$ is the same as the sign of $\widehat{f_4}$, while $\lvert Z\rvert$ measures how many standard errors the estimate lies from zero.

For example, $Z=-3$ means that the $f_4$-estimate is three standard errors below zero and corresponds to a two-sided $p$-value of approximately $0.0027$. More negative $Z$-scores provide stronger evidence that $B$ shares more drift with $C$ than with $D$.

There is no definitive significance threshold. A value of $\lvert Z\rvert\ge3$ is often used as a convention, lower absolute $Z$-scores could still be considered depending on the aim of the analysis. The $Z$-score is a useful summary metric because it represents the estimate relative to its standard error. It is also usually on a more convenient numerical scale than the small decimal values of the $f_4$-estimates themselves, which makes results easier to compare.

f4-Examples with AdmixPy

In this previous post, I gave already some general absolute drift share examples in “admixtools2”, usually against only a single axis which mostly focus on sign convention.

The population names in the examples are dataset-specific. In most .ind files, the population label is in the third column, so replace the names below with the labels used in your own file.

Below, I use AdmixPy. For installation instructions, see Introducing AdmixPy. If you already have AdmixPy installed, update it to the latest version first.

After activating your virtual environment, start a Python REPL by entering python in the terminal. Then import AdmixPy and set the dataset prefix:

import admixpy as ap

prefix = "v66_compatibility"

1. Which Northeast Asian source best represents the non-ANE ancestry in Tarim EMBA1?

The Early Bronze Age Tarim Basin population, known from the well-preserved Tarim mummies, derived most of its ancestry from Ancient North Eurasians but also had a Northeast Asian-related contribution. To test which reference groups share additional drift with Tarim EMBA1 relative to Afontova Gora, run:

ap.f4(
    prefix,
    "Chimp",
    [
        "Russia_Vologda_Mesolithic",
        "Turkey_Epipaleolithic",
        "Iran_BeltCave_Mesolithic",
        "China_TianyuanCave_UP",
        "China_AmurRiverBasin_N",
        "Russia_PrimorskyKrai_AmurRiver_N",
        "USA_WA_Kennewick_8800BP",
        "China_AmurRiverBasin_Mesolithic",
    ],
    "Tarim_EMBA1",
    "AfontovaGora_UP",
)

I use Chimp as the outgroup, although an African outgroup such as Mbuti would also work here. For pop2, I pass a Python list with the groups to test in position $B$. Tarim EMBA1 is in position $C$, and the ANE-related Afontova Gora group in position $D$. Swapping these last two populations would reverse the signs.

This returns:

    pop1                              pop2         pop3             pop4           est           se      z         p       n
0  Chimp         Russia_Vologda_Mesolithic  Tarim_EMBA1  AfontovaGora_UP   0.000195698  0.000486525    0.4     0.688  164434
1  Chimp             Turkey_Epipaleolithic  Tarim_EMBA1  AfontovaGora_UP   0.000531405  0.000673312   0.79      0.43  141828
2  Chimp          Iran_BeltCave_Mesolithic  Tarim_EMBA1  AfontovaGora_UP  -0.000203625  0.000895436  -0.23      0.82   63899
3  Chimp             China_TianyuanCave_UP  Tarim_EMBA1  AfontovaGora_UP   -0.00158086  0.000593296  -2.66     0.008  150232
4  Chimp            China_AmurRiverBasin_N  Tarim_EMBA1  AfontovaGora_UP   -0.00393877  0.000541723  -7.27  3.57e-13  153428
5  Chimp  Russia_PrimorskyKrai_AmurRiver_N  Tarim_EMBA1  AfontovaGora_UP   -0.00353571  0.000539207  -6.56  5.48e-11  164427
6  Chimp           USA_WA_Kennewick_8800BP  Tarim_EMBA1  AfontovaGora_UP   7.71892e-05  0.000917437   0.08     0.933   62992
7  Chimp   China_AmurRiverBasin_Mesolithic  Tarim_EMBA1  AfontovaGora_UP   -0.00385563  0.000615036  -6.27  3.63e-10  132437

As expected, the West Eurasian and Kennewick comparison populations yield $Z$-scores near zero. None shares significantly more drift with Tarim EMBA1 than with the Ancient North Eurasian Afontova Gora group.

The three Amur-related populations in position $B$ have significant negative $Z$-scores. Tianyuan also produces a negative result ($Z=-2.66$) that approaches the conventional significance threshold and points in the direction of excess East Asian-related affinity in Tarim EMBA1. As explained earlier, a negative result in this arrangement indicates excess affinity between $B$ (pop2) and $C$ (pop3).

Differences in SNP coverage (column n) can affect $f_4$-estimates and their $Z$-scores. The SNP counts above are not highly imbalanced, but I nevertheless reran the statistics for only the three significantly deviating groups with allsnps=False. This forces all three statistics to use the same intersecting set of SNPs:

ap.f4(
    prefix,
    "Chimp",
    [
        "China_AmurRiverBasin_N",
        "Russia_PrimorskyKrai_AmurRiver_N",
        "China_AmurRiverBasin_Mesolithic",
    ],
    "Tarim_EMBA1",
    "AfontovaGora_UP",
    allsnps=False,
)

This results in:

    pop1                              pop2        pop3              pop4          est           se      z         p       n
0  Chimp            China_AmurRiverBasin_N  Tarim_EMBA1  AfontovaGora_UP  -0.00390121  0.000559719  -6.97  3.17e-12  130128
1  Chimp  Russia_PrimorskyKrai_AmurRiver_N  Tarim_EMBA1  AfontovaGora_UP   -0.0037304  0.000592335   -6.3  3.02e-10  130128
2  Chimp   China_AmurRiverBasin_Mesolithic  Tarim_EMBA1  AfontovaGora_UP  -0.00387299  0.000616228  -6.28  3.28e-10  130128

Amur River Basin Neolithic still produces the largest absolute $Z$-score, although the differences among the three results are modest.

These findings could be used for a qpAdm model. I would test Amur River Basin Neolithic alongside Afontova Gora as source populations and place Amur River Basin Mesolithic among the right populations to anchor the Northeast Asian source.

This approach can also be extended into a proximity ranking. Instead of comparing Tarim EMBA1 with a single population in pop4, one can pass a list of candidate populations and use a less redundant, broader set of relevant references in pop2. The resulting statistics can then be summarised with pandas and numpy, both of which are already installed as AdmixPy dependencies. Candidates whose $f_4$-statistics show the smallest overall deviations from zero are the most similar to the target relative to that reference panel: the candidate and target have the most symmetric relationships to the populations in pop2.

2. Phylogenetic placement and unrooted pairings

An outgroup is not required. Under a simple tree without admixture, $f_4(A,B;C,D)$ is expected to be zero when $A$ and $B$ lie on one side of the internal split and $C$ and $D$ lie on the other.

For example:

ap.f4(prefix, "Andaman_100BP", "Onge", "Papuan", "Han")

            pop1  pop2    pop3 pop4          est           se     z      p       n
0  Andaman_100BP  Onge  Papuan  Han  0.000207602  0.000267998  0.77  0.439  697025

This near-zero result represents the unrooted split $(\text{Andaman},\text{Onge})\mid(\text{Papuan},\text{Han})$. Placing the closely related Andaman and Onge in positions $A$ and $C$ instead:

ap.f4(prefix, "Andaman_100BP", "Papuan", "Onge", "Han")

returns:

            pop1    pop2  pop3 pop4         est          se      z          p       n
0  Andaman_100BP  Papuan  Onge  Han  0.00971223  0.00033594  28.91  8.79e-184  697025

This large positive value shows that the combined covariance of the Andaman–Onge and Papuan–Han pairings is greater than that of the Andaman–Han and Papuan–Onge pairings. So it supports the same unrooted split as the first result.

The outgroup-based interpretation would be misleading here because Andaman is closely related to Onge and does not share approximately equal drift with Onge and Han. If Andaman were incorrectly treated as an outgroup, the second result might be read as evidence that Papuan shares more drift with Han than with Onge. Replacing Andaman with Chimp tests that interpretation:

ap.f4(prefix, "Chimp", "Papuan", "Onge", "Han")

    pop1    pop2  pop3 pop4          est           se     z      p       n
0  Chimp  Papuan  Onge  Han  0.000382461  0.000314833  1.21  0.224  679634

With a suitable outgroup, there is no significant evidence that Papuan shares meaningfully more drift with Han than with Onge. The large value in the previous arrangement was therefore driven by the Andaman–Onge relationship.

In qpGraph, these results would justify an internal node joining Andaman 100BP and Onge. Placing that branch, together with Han and Papuan, under a shared East Eurasian ancestral node would come from prior expectations rather than anything directly tested here. Tianyuan and Ust-Ishim might be useful reference groups for testing that deeper placement. Additional tests would be needed to place this structure relative to deeper Eurasian groups and to identify admixture edges.

3. Detecting deviations from a population tree with $f_4$-statistics

As mentioned above, under a simple bifurcating tree without admixture, one of the three possible arrangements of four populations should produce a non-significant $f_4$-statistic around zero. If none does, their relationships cannot be represented by a single unrooted split. The following three arrangements test Kotias Klde Mesolithic, Samara Yamnaya, Barcin Neolithic, and Vologda Mesolithic.

ap.f4(
    prefix,
    "Georgia_KotiasKlde_Mesolithic",
    "Turkey_Barcin_Neolithic-DG",
    "Russia_Samara_EBA_Yamnaya",
    "Russia_Vologda_Mesolithic",
)

ap.f4(
    prefix,
    "Georgia_KotiasKlde_Mesolithic",
    "Russia_Samara_EBA_Yamnaya",
    "Turkey_Barcin_Neolithic-DG",
    "Russia_Vologda_Mesolithic",
)

ap.f4(
    prefix,
    "Georgia_KotiasKlde_Mesolithic",
    "Russia_Vologda_Mesolithic",
    "Turkey_Barcin_Neolithic-DG",
    "Russia_Samara_EBA_Yamnaya",
)

To cover the three permutations, one population can remain fixed in position $A$, while each of the other three appears once in position $B$. The other two populations occupy positions $C$ and $D$.

These return:

test          est           se       z          p       n
1     0.000879912  0.000115676    7.61   2.81e-14  1241213
2      0.00364597  0.000167812   21.73  1.15e-104  1241213
3      0.00276606  0.000168548   16.41   1.59e-60  1241213

The three statistics are permutations of the same quartet and are not independent. Together, however, they show that no single unrooted split fits these four populations. Based on prior knowledge, the likely admixed population here is the Bronze Age Yamnaya group. These tests cannot determine which or how many populations are admixed, the direction of gene flow, or the mixture proportions. For that an admixture model with qpAdm or qpGraph is needed.

4. Outgroups and ancestry-pole contrasts in a population with sub-Saharan African admixture

To compare the modern BedouinB group with the Early Medieval Bedouin group from Tell Qarassa along several ancestry axes, I first use Chimp in position $A$, as in the first example:

ap.f4(
    prefix,
    "Chimp",
    [
        "Turkey_Central_Boncuklu_PPN",
        "Iran_GanjDareh_N",
        "Jordan_PPNB",
        "Dinka",
        "Yoruba",
        "Morocco_Iberomaurusian",
        "Natufian",
    ],
    "Syria_TellQarassa_EarlyMedieval",
    "BedouinB",
)

This returns:

    pop1                         pop2                             pop3      pop4           est           se      z         p       n
0  Chimp  Turkey_Central_Boncuklu_PPN  Syria_TellQarassa_EarlyMedieval  BedouinB    -0.0020768  0.000345682  -6.01  1.88e-09  671213
1  Chimp             Iran_GanjDareh_N  Syria_TellQarassa_EarlyMedieval  BedouinB   -0.00105039  0.000339625  -3.09     0.002  666707
2  Chimp                  Jordan_PPNB  Syria_TellQarassa_EarlyMedieval  BedouinB   -0.00197691  0.000357698  -5.53  3.26e-08  530640
3  Chimp                        Dinka  Syria_TellQarassa_EarlyMedieval  BedouinB  -8.93486e-05  0.000283272  -0.32     0.752  677979
4  Chimp                       Yoruba  Syria_TellQarassa_EarlyMedieval  BedouinB    0.00013197  0.000203751   0.65     0.517  678088
5  Chimp       Morocco_Iberomaurusian  Syria_TellQarassa_EarlyMedieval  BedouinB  -0.000985543  0.000336736  -2.93     0.003  624606
6  Chimp                     Natufian  Syria_TellQarassa_EarlyMedieval  BedouinB   -0.00256286  0.000385414  -6.65  2.94e-11  427631

The negative results for the Near Eastern references show that they share more drift with Tell Qarassa than with BedouinB. African admixture could explain the reduced proportion of Near Eastern ancestry in BedouinB, although other explanations are possible. With Chimp as the outgroup, BedouinB shares notably less drift with the tested Near Eastern references than Tell Qarassa does. Still BedouinB is not significantly closer to either Dinka or Yoruba. For a predominantly West Eurasian population, this combination suggests ancestry that shifts BedouinB away from the Near Eastern profile but is possibly not well represented by either African reference. Sub-Saharan African admixture is a plausible explanation, although the Chimp-based tests alone do not establish it.

A more sensitive contrast can be constructed by replacing Chimp with Ust-Ishim:

ap.f4(
    prefix,
    "Russia_UstIshim_IUP",
    [
        "Turkey_Central_Boncuklu_PPN",
        "Iran_GanjDareh_N",
        "Jordan_PPNB",
        "Dinka",
        "Yoruba",
        "Morocco_Iberomaurusian",
        "Natufian",
    ],
    "Syria_TellQarassa_EarlyMedieval",
    "BedouinB",
)

Now the returned table is:

    pop1                              pop2         pop3                    pop4           est           se      z         p       n
0  Russia_UstIshim_IUP  Turkey_Central_Boncuklu_PPN  Syria_TellQarassa_EarlyMedieval  BedouinB  -0.000452437  0.000229414  -1.97     0.049  1216469
1  Russia_UstIshim_IUP             Iran_GanjDareh_N  Syria_TellQarassa_EarlyMedieval  BedouinB   0.000233737  0.000230595   1.01     0.311  1210551
2  Russia_UstIshim_IUP                  Jordan_PPNB  Syria_TellQarassa_EarlyMedieval  BedouinB  -0.000211399  0.000258228  -0.82     0.413   917557
3  Russia_UstIshim_IUP                        Dinka  Syria_TellQarassa_EarlyMedieval  BedouinB    0.00138816  0.000400902   3.46  0.000535   705767
4  Russia_UstIshim_IUP                       Yoruba  Syria_TellQarassa_EarlyMedieval  BedouinB    0.00118368  0.000234436   5.05  4.44e-07  1232680
5  Russia_UstIshim_IUP       Morocco_Iberomaurusian  Syria_TellQarassa_EarlyMedieval  BedouinB   0.000333061  0.000242472   1.37      0.17  1055263
6  Russia_UstIshim_IUP                     Natufian  Syria_TellQarassa_EarlyMedieval  BedouinB  -0.000703807  0.000270607   -2.6     0.009   711444

Ust-Ishim is approximately equally related to the West Eurasian ancestry shared by Tell Qarassa and BedouinB. The additional African ancestry in BedouinB nevertheless reduces its overall affinity to Ust-Ishim. Placing Dinka or Yoruba on the opposite side of the comparison makes this African shift easier to detect. The significant positive results place BedouinB closer to the sub-Saharan African side of the contrast than Tell Qarassa, although they do not identify the exact African source or estimate its proportion. The Near Eastern references now also remain closer to zero. Boncuklu and especially Natufian still show a tendency toward greater affinity with Tell Qarassa.

This ancestry-pole setup is mainly useful when the two populations differ along an African–non-African axis. For comparisons among predominantly Eurasian populations, Chimp or an African population will usually be approximately symmetric to both, making the standard outgroup approach sufficient.

In another post, I will build on these tests to model the examples used here with qpAdm.