Interpreting f4-Statistics with AdmixPy

f4-statistics can be used to test asymmetries in allele sharing between populations. They measure the covariance between allele-frequency differences across two pairs of populations.

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$. A larger covariance indicates that the two populations share more allele-frequency changes and more shared genetic drift, while a smaller covariance indicates less shared drift.

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 pairings.

These are relative statements about the two sums of covariance terms; without further assumptions, the statistic cannot identify which 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 each other 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. 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.

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.

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. Replace them with the exact labels in the third column of your .ind 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 a outgroup here, although an African outgroup such as Mbuti would also work for this example. 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, as would swapping the outgroup from position $A$ to $B$.

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 reference 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 has 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 for the significant statistics above are not highly imbalanced, but I nevertheless reran the tests with allsnps=False for those three groups. 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 has 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 Mbuti tests that interpretation:

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

    pop1    pop2  pop3 pop4           est          se      z      p       n
0  Mbuti  Papuan  Onge  Han  -0.000253514  0.00027654  -0.92  0.359  706613

With a suitable outgroup, there is no significant evidence that Papuan shares more drift with Han than with Onge. The large value in the previous arrangement was therefore driven by the Andaman–Onge relationship. This direct arrangement tests whether Han is closer to Onge or Papuan:

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

returns:

    pop1 pop2  pop3    pop4          est          se      z        p       n
0  Mbuti  Han  Onge  Papuan  -0.00216992  0.00030816  -7.04  1.9e-12  706613

With Mbuti as the outgroup, the significantly negative value indicates that Han shares more drift with Onge than with Papuan ($Z=-7.04$). This does not contradict the earlier unrooted split: Papuan and Han can lie on the same side of that split without forming a clade in the rooted tree. A topology in which Papuan branches first, then Han, with Andaman and Onge as sisters, would result in the same unrooted split when Mbuti is excluded.

In qpGraph, these results support testing an internal node joining Andaman 100BP and Onge, with Han joining their lineage before Papuan. Placing these branches under a shared East Eurasian ancestral node would still require additional testing rather than following directly from these results. Tianyuan and Ust-Ishim might be useful reference groups for testing that deeper placement.

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 result in 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 should be placed once in position $B$. The other two populations occupy positions $C$ and $D$.

These return:

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

The three statistics are permutations of the same quartet. They indicate that no single unrooted split fits these four populations, meaning that at least one population is admixed. Based on prior knowledge, the likely admixed population here is the Bronze Age Yamnaya group.

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 Mbuti in position $A$, as in the first example:

ap.f4(
    prefix,
    "Mbuti",
    [
        "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  Mbuti  Turkey_Central_Boncuklu_PPN  Syria_TellQarassa_EarlyMedieval  BedouinB   -0.00163561  0.000199886  -8.18  2.78e-16  1216503
1  Mbuti             Iran_GanjDareh_N  Syria_TellQarassa_EarlyMedieval  BedouinB  -0.000982149  0.000189909  -5.17  2.32e-07  1210579
2  Mbuti                  Jordan_PPNB  Syria_TellQarassa_EarlyMedieval  BedouinB   -0.00144038  0.000208195  -6.92  4.57e-12   917580
3  Mbuti                        Dinka  Syria_TellQarassa_EarlyMedieval  BedouinB   -0.00031521  0.000223465  -1.41     0.158   705778
4  Mbuti                       Yoruba  Syria_TellQarassa_EarlyMedieval  BedouinB  -4.26945e-05  7.75272e-05  -0.55     0.582  1232719
5  Mbuti       Morocco_Iberomaurusian  Syria_TellQarassa_EarlyMedieval  BedouinB  -0.000821848  0.000187372  -4.39  1.15e-05  1055293
6  Mbuti                     Natufian  Syria_TellQarassa_EarlyMedieval  BedouinB   -0.00175053  0.000209975  -8.34  7.63e-17   711464

The negative results for the Near Eastern references indicate that they share more drift with Tell Qarassa than with BedouinB. African admixture could explain the reduced Near Eastern affinity of BedouinB. However, 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. Sub-Saharan African admixture is a plausible explanation, although the Mbuti-based tests alone do not seem to indicate that.

A more sensitive contrast can be constructed by replacing Mbuti 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. The Near Eastern references now also remain closer to zero. Boncuklu and especially Natufian still show a tendency toward greater affinity with Tell Qarassa. Because of the imbalance in SNP counts across the SSA groups, it could be worth rerunning them separately with allsnps=Falsehere.

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

In another post, I’ll use qpAdm to model the examples shown here.