Testing Ancestry Models with qpAdm in AdmixPy: A Tarim Basin Example

qpAdm tests whether a target can be modeled as a mixture of proposed source populations. It uses $f_4$-statistics to compare how the target and sources relate to a set of reference populations, then estimates mixture proportions and evaluates how well the proposed model fits the data.

Data and prerequisites

To run qpAdm models, you need a genetic dataset. See Downloading and Converting AADR v66 for the dataset used here.

I’ll use AdmixPy to run the models. For installation instructions, see Introducing AdmixPy.

You can paste the code below into a .py file (and run it with python run.py or python3 run.py) or run it in the Python REPL. Jupyter Notebook offers a useful compromise between easily navigable code, interactive testing, and variables that remain available during the session. Setup instructions can be found on their website.

Example: Modeling the Tarim EMBA1 group

For qpAdm, you need specific hypotheses about ancestry that you want to test.

For this example, I’ll model the Tarim EMBA1 group. Earlier analysis suggested that this group is genetically similar to the Afontova Gora group, with additional Northeast Asian ancestry.

Target and source populations

To select groups for the analysis, use the exact population labels from your dataset’s .ind file. The labels here are specific to my dataset.

You can pass the population labels directly to qpadm(), but I’ll define them as variables here:

import admixpy as ap

data = "dataset_prefix"

target = "Tarim_EMBA1"
left = ["AfontovaGora_UP", "China_AmurRiverBasin_N"]

Replace "dataset_prefix" with your dataset’s path and prefix, without a file extension. The target is a single string in quotes, while the sources are a list of strings in square brackets. I’ve named the source variable left.

Right populations

As mentioned, right populations are the references that provide the $f_4$-statistics used to fit the model.

With a fixed reference in position $A$, the arrangement can be written as:

Here, $R_0$ is the fixed right population, $R_j$ runs over the other right populations, $T$ is the target, and $S_i$ runs over the sources.

The right groups should be differentially related to the sources and target, providing contrasts that help distinguish their ancestry. I prefer to use a distant outgroup as the fixed reference $R_0$.

You can think of these contrasts as anchoring ancestry axes. Without references that distinguish the proposed sources, their weights can be poorly resolved, with large standard errors. The previous qpAdm article shows an example of this.

Here, I’ll start with three groups:

refs = ["Chimp", "Russia_Malta_UP", "China_AmurRiverBasin_Mesolithic"]

Mal’ta UP is an Ancient North Eurasian group related to the Afontova Gora source. The Mesolithic Amur River Basin group comes from the same region as the Amur Neolithic source and predates it. I use these references to help distinguish the two sources. As the first element in refs, Chimp is automatically used as the fixed right population $R_0$.

Running qpAdm

Does Afontova Gora fit alongside Amur Neolithic?

The earlier $f_4$ tests found excess Amur-related affinity in Tarim EMBA1 relative to Afontova Gora. That motivates the first model: Afontova Gora plus Amur Neolithic, using Mal’ta and Amur Mesolithic as references for these two ancestry components.

With the variables defined, we can call qpadm():

result = ap.qpadm(data, target, left, refs)
print(result)

The arguments are passed in the order data, target, left, and right. Alternatively, you can name them explicitly, for example ap.qpadm(data, target=target, left=left, right=refs). Named arguments can appear in any order after data.

In a REPL or notebook, entering result displays the output without print(). You can also call ap.qpadm(data, target, left, refs) directly to display the result without assigning it to a variable.

This gives the following output:

QpAdmResult(target='Tarim_EMBA1')

weights:
                  left weight    se     z
       AfontovaGora_UP  0.821 0.016 49.96
China_AmurRiverBasin_N  0.179 0.016 10.91

rankdrop:
f4rank dof   chisq     p p_nested
     1   1    1.62 0.203        0
     0   4 3321.07     0      NaN

popdrop:
pat                dropped f4rank dof  chisq        p  feasible status
 00                             1   1   1.62    0.203      True   PASS
 01 China_AmurRiverBasin_N      0   2  89.57 3.54e-20      True   FAIL
 10        AfontovaGora_UP      0   2 2989.6        0      True   FAIL

In this model, Tarim EMBA1 is estimated to have $82.1\% \pm 1.6\%$ Afontova Gora-related Ancient North Eurasian ancestry and $17.9\% \pm 1.6\%$ Amur Neolithic-related Northeast Asian ancestry (± one standard error).

The feasible column indicates whether the estimated weights lie between 0 and 1. Here, PASS marks the feasible model that is not rejected at the $0.05$ threshold.

The model is not rejected at the commonly used $p = 0.05$ threshold, with $p = 0.203$. Dropping either source leads to rejection of the remaining one-source model.

This model has only one degree of freedom and is evaluated using the genetic contrasts provided by just three right populations. Adding references that capture additional ancestry relationships can reveal discrepancies that the initial set cannot detect as will be seen below.

The chi-square statistic ($\chi^2 = 1.62$) measures the discrepancy between observed and fitted $f_4$-statistics, relative to their uncertainty and covariance. Together with the degrees of freedom (dof), it determines the p-value.

Do Baikal and Irkutsk EBA sources fit better?

Zhang et al. (2021) modeled Tarim ancestry using Afontova Gora and Baikal_EBA as proxies. For comparison, I keep Afontova Gora and the initial three references unchanged and replace Amur Neolithic with the following Baikal/Irkutsk EBA groups:

The Irkutsk and Shamanka models are not rejected at the 0.05 threshold, although their p-values are relatively close to it. The Kurma model is rejected. All three have larger discrepancies and lower p-values than the Amur Neolithic model in this setup. This is a comparison using my reference set, not an exact replication of the published model. However, I do not see evidence here that the Northeast Asian contribution must have come specifically from such Baikal EBA groups.

Does the model also fit with Vologda as a reference?

To test the first model further, I add Vologda Mesolithic, a western Eurasian hunter-gatherer group with substantial ANE-related ancestry. It provides an additional contrast that can reveal differences not captured by the initial ANE and Northeast Asian references.

In the earlier $f_4$ analysis, Vologda showed no significant excess affinity to Tarim EMBA1 relative to Afontova Gora ($Z = 0.4$). Vologda can still reveal differences when either group is compared with Amur Mesolithic. The proposed mixture must also reproduce these relationships.

# or use refs.append("Russia_Vologda_Mesolithic") in the REPL
refs = ["Chimp", "Russia_Malta_UP", "China_AmurRiverBasin_Mesolithic", "Russia_Vologda_Mesolithic"]

Rerunning the model with this expanded reference set results in:

QpAdmResult(target='Tarim_EMBA1')

weights:
                  left weight    se     z
       AfontovaGora_UP  0.874 0.013 65.58
China_AmurRiverBasin_N  0.126 0.013  9.47

rankdrop:
f4rank dof   chisq        p p_nested
     1   2   22.76 1.14e-05        0
     0   6 5643.07        0      NaN

popdrop:
pat                dropped f4rank dof   chisq        p  feasible status
 00                             1   2   22.76 1.14e-05      True   FAIL
 01 China_AmurRiverBasin_N      0   3   95.54 1.41e-20      True   FAIL
 10        AfontovaGora_UP      0   3 5127.36        0      True   FAIL

The model is rejected at $p = 1.14 \times 10^{-5}$, even though its weights remain feasible. The estimates shift to 87.4% Afontova Gora and 12.6% Amur Neolithic. The initial pass therefore depended on a reference set that did not capture this mismatch.

This rejection does not by itself necessitate a Vologda-related ancestry contribution. I also tried adding Vologda Mesolithic as a third source and replacing it in the references with another Eastern Hunter-Gatherer (EHG) group. That model was feasible but still rejected with no improvement in the fit.

Does Tutkaul fit better as the ANE-related source?

Another possibility is that Afontova Gora is an imperfect proxy for the ANE-related ancestry. To test an alternative, I return to two sources and replace it with the ANE-rich Central Asian Tutkaul Mesolithic group:

left = ["Tutkaul_Mesolithic", "China_AmurRiverBasin_N"]

Rerunning qpadm() with the four-reference set above gives:

QpAdmResult(target='Tarim_EMBA1')

weights:
                  left weight    se     z
    Tutkaul_Mesolithic  0.961 0.017 56.91
China_AmurRiverBasin_N  0.039 0.017   2.3

rankdrop:
f4rank dof   chisq   p p_nested
     1   2     4.6 0.1        0
     0   6 5286.25   0      NaN

popdrop:
pat                dropped f4rank dof   chisq     p  feasible status
 00                             1   2     4.6   0.1      True   PASS
 01 China_AmurRiverBasin_N      0   3    9.44 0.024      True   FAIL
 10     Tutkaul_Mesolithic      0   3 5127.36     0      True   FAIL

With the same four references, the Tutkaul model is not rejected ($p = 0.1$). It estimates 96.1% Tutkaul-related and 3.9% Amur Neolithic-related ancestry.

Does Tutkaul also fit with an Iranian Plateau reference?

To check whether Tutkaul shares more drift with Iranian Plateau groups than Tarim EMBA1 does, I run:

ap.f4(
    data,
    "Chimp",
    ["Iran_BeltCave_Mesolithic", "Iran_GanjDareh_N"],
    "Tutkaul_Mesolithic",
    "Tarim_EMBA1",
)

    pop1                      pop2                pop3         pop4          est           se      z         p       n
0  Chimp  Iran_BeltCave_Mesolithic  Tutkaul_Mesolithic  Tarim_EMBA1  -0.00162587   0.00102051  -1.59     0.111   52835
1  Chimp          Iran_GanjDareh_N  Tutkaul_Mesolithic  Tarim_EMBA1  -0.00210138  0.000538013  -3.91  9.39e-05  133886

Negative values in this arrangement indicate greater affinity between the Iranian reference and Tutkaul than between that reference and Tarim EMBA1. The Ganj Dareh result supports this difference ($Z = -3.91$), Belt Cave points in the same direction but is not strongly significant ($Z = -1.59$).

In the earlier $f_4$ post, Belt Cave also showed no significant excess affinity to Tarim EMBA1 relative to Afontova Gora ($Z = -0.23$). An Iranian Plateau reference therefore provides a useful additional test of the Tutkaul model.

I add Belt Cave Mesolithic to the references, keeping the two sources unchanged:

refs = ["Chimp", "Russia_Malta_UP", "China_AmurRiverBasin_Mesolithic", "Russia_Vologda_Mesolithic", "Iran_BeltCave_Mesolithic"]

Rerunning qpadm() gives:

QpAdmResult(target='Tarim_EMBA1')

weights:
                  left weight    se     z
    Tutkaul_Mesolithic  0.959 0.017 57.02
China_AmurRiverBasin_N  0.041 0.017  2.41

rankdrop:
f4rank dof   chisq    p p_nested
     1   3   11.33 0.01        0
     0   8 5297.04    0      NaN

popdrop:
pat                dropped f4rank dof   chisq     p  feasible status
 00                             1   3   11.33  0.01      True   FAIL
 01 China_AmurRiverBasin_N      0   4   16.62 0.002      True   FAIL
 10     Tutkaul_Mesolithic      0   4 5134.52     0      True   FAIL

The model is now rejected at $p = 0.01$, while the weights remain close to the previous estimates. Tutkaul therefore fits the smaller reference set but fails when the Iranian Plateau contrast is included.

These results suggest that a closer proxy could have been a more eastern, ANE-rich Mesolithic group with additional Northeast Asian ancestry relative to Upper Paleolithic Ancient North Eurasians. Such a proxy does not seem to be available among currently sampled populations.

Model comparison

The initial reference set is Chimp, Russia_Malta_UP, and China_AmurRiverBasin_Mesolithic. The table summarizes the models with output shown above; all use Tarim EMBA1 as the target.

Rotating qpAdm can automate comparisons by moving candidate populations between the source and reference sets. I’ll cover this in another article.