qpAdm analysis tutorial: from a raw DNA file to a model with a p-value

This is the tutorial I wish had existed when I first tried to run qpAdm on a consumer DNA file. It goes from the raw export a testing company gives you to a model with a p-value you can defend, and it stops at every point where the method can quietly go wrong. Nothing here requires a purchase: the file check and the browser workbench are free, and the R code runs on any laptop. If you would rather have the whole thing done and audited for you, that is what a qpAdm analysis from us is, and the last section explains how the two paths meet. But read the tutorial first, because even a customer reads a report better after having watched a model fail.

What qpAdm actually tests#

qpAdm is a method from the ADMIXTOOLS package, introduced in the supplementary material of Haak and colleagues (2015) and now the standard admixture-modelling tool in ancient-DNA research. It takes a target (here, your genome), a proposed set of source populations (the "left" set), and a set of outgroups (the "right" set), and asks one question: can the target be written as a mixture of the sources, given how each of them relates to the outgroups?

The mechanism is f4-statistics. For every source and every outgroup the method measures how much genetic drift the pair shares, then checks whether the target's pattern of shared drift with the outgroups is a weighted average of the sources' patterns. If it is, the model is admissible and the weights are the mixture proportions. If it is not, the model is rejected. That rejection is the whole point: a p-value that can say no is what separates a formal test from a percentage generator. The fuller treatment is in Understanding qpAdm.

Harney and colleagues (2021) tested the method on simulated histories and found it reliable when sources and outgroups are chosen with care, and easy to fool when they are not. Two of their findings shape every step below. First, an outgroup that shares recent drift with a source breaks the test in both directions, rejecting good models and admitting bad ones. Second, if you rotate sources and outgroups until something passes, something will pass, and the model that passes is often wrong. A tutorial that skipped those two facts would be teaching you to produce confident nonsense.

The reproducible workflow#

Seven steps, in an order that does not change:

1. Check the raw file and learn its coverage.
2. Merge it into a reference panel (AADR v66).
3. Choose the left sources and the right outgroups.
4. Run the model, in the browser or in R.
5. Read the p-value, the weights, the standard errors, the Z-scores and the nested-model table.
6. Handle the rejections.
7. Write the model record.

Each step has a check you can perform before moving on. If you skip a check, the later steps will still produce numbers; they will just not mean anything.

Step 1: the raw file and its coverage#

A raw-data export from 23andMe, AncestryDNA, MyHeritage, FamilyTreeDNA or LivingDNA is a text table: one row per marker, with a chromosome, a position and a genotype. Depending on the vendor and the chip version it holds somewhere between 600,000 and 700,000 markers on a modern array, fewer on older ones. Whole-genome sequencing gives a VCF instead, which holds far more positions but needs converting to the reference panel's markers before any of this works.

The number that matters is not the marker count in the file but the number of markers that overlap the reference panel after the merge. That overlap is your coverage, and it bounds everything downstream. Standard errors in qpAdm are estimated by block jackknife across the genome; the fewer informative positions, the wider the jackknife distribution, the larger the standard error. No amount of skill in step 3 can shrink an error that step 1 has already fixed.

Before doing anything else, run the file through the free Raw DNA File Check. It parses the export, reports the marker count, the build (GRCh37 or GRCh38) and the expected overlap with AADR, and tells you whether the file can support a model with standard errors under 0.10 at all. The file is analysed and discarded. If the check says the coverage is thin, the honest response is to know that now, not after two hours in R.

Two things the check will catch that people miss: a file that has been re-saved through a spreadsheet program and lost its position column formatting, and a "raw" file that is in fact a vendor's imputed output with far more rows than the chip ever read. Both look fine to the eye and neither merges cleanly.

Step 2: merging into AADR v66#

The Allen Ancient DNA Resource (AADR) is the curated compendium of published ancient and modern genomes maintained by the Reich laboratory (Mallick and colleagues, 2024). Version 66 holds roughly 23,265 individuals across roughly 6,015 population labels, in EIGENSTRAT format: a .geno file of genotypes, a .snp file of marker positions and an .ind file of individual labels. The whole thing is explained in The AADR explained.

A merge is the operation that puts your genotypes and the panel's at the same positions in the same file. Concretely:

• your file is converted from the vendor text layout to EIGENSTRAT, with the chromosome and position columns lifted to the panel's build where needed;
• indels and strand-ambiguous SNPs (A/T and C/G) are dropped, because they cannot be aligned unambiguously between two datasets;
• the remaining positions are intersected with the panel's .snp list, and only the intersection survives;
• your sample is appended as one new individual in the .ind file, with a population label of its own.

The intersection is where coverage becomes a hard number. A 650,000-marker consumer file intersected with the 1240K panel typically keeps somewhere in the low hundreds of thousands of positions; older chips keep fewer. That is what the file check estimates in step 1. If you are doing the merge yourself, Poseidon's trident is the tool that makes this reproducible, and a VCF needs an extra conversion step first; the details of that conversion are in Upload a whole-genome VCF.

The check for this step: open the merged .ind file and confirm your sample is present with the label you gave it, then count the lines of the merged .snp file. That count is the number you will report as "merged SNPs" in step 7.

Step 3: choosing left sources and right outgroups#

This is the step that decides whether the result means anything. The rules are short and every one of them has a failure mode behind it.

The left set (sources)#

Sources are the populations you propose the target descends from. Start with two or three. Keep them distinct: two closely related sources inflate each other's standard errors because the data cannot tell them apart. Keep them in one era: a model that mixes a Mesolithic source with a Roman-period one asks the method to treat a later, already-mixed population as an ancestor of an earlier structure, and it will usually reject that, correctly. Respect chronology: a source that postdates the target cannot be its ancestor, however good the fit.

For a European target, a defensible first left set is the classic three-way split used across the literature since Haak (2015): Anatolia_N (early Anatolian farmers), WHG (western hunter-gatherers) and Yamnaya_Samara (Bronze Age steppe pastoralists). For targets further east, Iran_N, CHG (Caucasus hunter-gatherers) and EHG (eastern hunter-gatherers) enter the picture; for the Levant, Levant_N.

The right set (outgroups)#

The right set gives the method its power to tell sources apart. Three rules:

• Distant. Every outgroup must be distant from every source. An outgroup that shares recent drift with one source (say, a Neolithic European population beside Anatolia_N) makes the f4-statistics that involve that source lie.
• Informative. The outgroups must, between them, be differentially related to the sources. If every outgroup is equally distant from every source, the sources are indistinguishable and the standard errors explode.
• Enough of them. The literature usually uses ten to fifteen. Six weak outgroups let almost anything pass; a p-value of 0.9 against six says less than 0.2 against fourteen.

A right set that satisfies all three for a European target, drawn from the AADR labels, is:

Mbuti.DG, Ust_Ishim.DG, Kostenki14, MA1, Han.DG, Papuan.DG, Karitiana.DG, Onge.DG,
Iran_N, Levant_N, CHG, EHG

The first eight are deep outgroups: an African population, a 45,000-year-old Siberian, an Upper Palaeolithic European, a 24,000-year-old Siberian, and four present-day populations from East Asia, Oceania, the Americas and the Andaman Islands. The last four are closer, and they are what makes the set informative: Iran_N and CHG are differentially related to Anatolia_N versus Yamnaya_Samara, and EHG is differentially related to Yamnaya_Samara versus WHG. That is exactly the leverage the model needs.

Rotate#

The sanity check on any left/right split is the rotating scheme: move a population from the right set into the left set and see whether it earns a place, or move a source to the right and see whether the model still passes without it. Rotation is a diagnostic, not a search. Used to test a model you already believe, it tells you how robust it is. Used to find a model, it is the practice Harney (2021) warns about: enough rotations will produce a passing model by chance.

The check for this step: write down, before running anything, why each outgroup is in the right set and which pair of sources it helps to separate. If you cannot answer for one of them, it does not belong there.

Step 4: run it#

In the browser, free, no R#

The free AdmixTools 2 Lab runs the real ADMIXTOOLS 2 package on our server over the public AADR Human Origins panel. You pick a target from the panel, add sources to the left set and outgroups to the right set, and press run; the output is the same table qpadm() prints in R, including the nested-model block. It needs a free, email-verified account and runs over the public panel only, so the target is a panel population rather than your own file. That is the right place to learn the method before spending an evening on a local install, and it is where the worked example below can be reproduced by anyone.

The R equivalent with ADMIXTOOLS 2#

ADMIXTOOLS 2 (Maier and colleagues, 2023) is an R package. Install it, point it at the merged EIGENSTRAT prefix from step 2, and the workflow is three calls.

install.packages("remotes")
remotes::install_github("uqrmaie1/admixtools")
library(admixtools)

First, extract the f2-statistics. This is the slow part: it reads the whole genotype file once and writes block-wise f2 values for every pair of populations you name to a directory, so that every later model runs in seconds.

prefix  <- "merged/my_kit_aadr_v66"    # .geno / .snp / .ind
f2_dir  <- "f2/my_kit"
 
pops <- c(
  "MyKit",
  "Anatolia_N", "WHG", "Yamnaya_Samara",
  "Mbuti.DG", "Ust_Ishim.DG", "Kostenki14", "MA1", "Han.DG",
  "Papuan.DG", "Karitiana.DG", "Onge.DG", "Iran_N", "Levant_N", "CHG", "EHG"
)
 
extract_f2(
  prefix, f2_dir,
  pops        = pops,
  maxmiss     = 0,        # only SNPs present in every population
  blgsize     = 0.05,     # 5 cM jackknife blocks
  overwrite   = TRUE
)

maxmiss = 0 keeps only positions typed in every population, which is the conservative default and the reason the "min SNPs per f4" figure in the output is often much lower than the merged SNP count. Raising it admits missing data, which ADMIXTOOLS 2 handles with a correction, at the cost of some comparability with the classic implementation.

Second, load the precomputed statistics and run the model.

f2 <- f2_from_precomp(f2_dir)
 
left  <- c("Anatolia_N", "WHG", "Yamnaya_Samara")
right <- c("Mbuti.DG", "Ust_Ishim.DG", "Kostenki14", "MA1", "Han.DG",
           "Papuan.DG", "Karitiana.DG", "Onge.DG", "Iran_N", "Levant_N",
           "CHG", "EHG")
 
fit <- qpadm(f2, left = left, right = right, target = "MyKit")

Third, read the output. fit is a list of tables:

fit$weights     # one row per source: weight, se, z
fit$rankdrop    # the rank test: f4rank, dof, chisq, p
fit$popdrop     # the nested-model table: every subset of sources
fit$f4          # the individual f4-statistics behind the fit

fit$weights is the table most people stop at, and step 5 explains why that is a mistake. fit$rankdrop carries the p-value for the full-rank model on its first row. fit$popdrop is the nested-model table: for every subset of the sources it reports the refitted weights, the p-value and whether the weights are feasible (all between 0 and 1). fit$f4 is what you open when a model is rejected and you want to know which outgroup did it.

Step 5: reading the result#

Read in this order, every time: p-value, then standard errors, then Z-scores, then the weights, then the nested-model table. The full reading guide is How to read qpAdm results; the short version follows.

The p-value asks whether the target's pattern of shared drift with the outgroups is compatible with the proposed mixture. Above 0.05 the model is not refuted. Below it, the model is rejected and the weights are not worth reading. It is not the probability the model is true, it is not comparable across different right sets, and a higher value is not a stronger result once the threshold is cleared.

The standard error on each weight comes from the block jackknife. Read every weight as a range: 0.37 with SE 0.034 is "roughly 30 to 44% at two standard errors". Coverage sets it; the analyst cannot.

The Z-score is the weight divided by its standard error. It tells you whether a source is distinguishable from contributing nothing. A source at Z = 1.5 is not measurably present, whatever its headline weight.

The weights, read last, are the mixture proportions for this model against this right set. Change either and they move.

The nested-model table is the check most tutorials skip. For every simpler model (each source removed in turn) it shows whether the simpler version also passes. If the model without your weakest source is admissible with feasible weights, that source was never earning its place, and the simpler model is the one to report.

The bar we publish against, at every tier of our own service, is p > 0.05, every |Z| > 3 and every SE < 0.10. Those thresholds are a choice, and a stricter one than the bare p > 0.05 most papers use, but they have the virtue of being stated before the run rather than after.

Step 6: the rejection cases, with a worked example#

What follows is an illustrative example, constructed to be internally consistent with the bar above. It is not any customer's result and the numbers are not taken from a real run. The target is a present-day genome from the Balkans, merged as in step 2, with the right set from step 3.

The accepted three-way model#

target: Balkan_Kit    left: Anatolia_N, WHG, Yamnaya_Samara
p-value: 0.21    f4 rank: 2    dof: 9

source            weight    SE      Z
Anatolia_N        0.52      0.031   16.8
WHG               0.11      0.028    3.9
Yamnaya_Samara    0.37      0.034   10.9

Reading it in order: p = 0.21 clears 0.05, so the model is admissible. Every SE is under 0.10. Every |Z| exceeds 3, and the smallest, WHG at 3.9, is the one to watch. The weights can now be read: about half early farmer, a little over a third steppe, a tenth western hunter-gatherer, which is a familiar shape for the region. The nested-model table is what confirms WHG's place:

nested model                   p-value   feasible
Anatolia_N + Yamnaya_Samara    0.004     yes
Anatolia_N + WHG               <0.001    no
WHG + Yamnaya_Samara           <0.001    no

The two-way model without WHG is rejected at p = 0.004. So WHG is not decoration: removing it breaks the fit. The three-way model is the one to report.

The rejected two-way alternative#

Suppose you had started with the simpler proposal, Anatolia_N + Yamnaya_Samara, on the reasoning that hunter-gatherer ancestry in the Balkans is small enough to ignore:

target: Balkan_Kit    left: Anatolia_N, Yamnaya_Samara
p-value: 0.004    f4 rank: 1    dof: 10

source            weight    SE      Z
Anatolia_N        0.58      0.029   20.0
Yamnaya_Samara    0.42      0.029   14.5

Every standard error is small and both Z-scores are enormous. A reader who skipped the p-value would report a clean 58/42 split. The p-value says the split is incompatible with the data. The weights are not worth reading, and the correct next move is to look at fit$f4 to see which outgroups the residual is concentrated on. In this constructed case they would be the ones that separate WHG from the other two sources, Kostenki14 and EHG, pointing directly at the missing source.

The other ways a model fails#

A bad outgroup. If Levant_N is replaced with a Neolithic European population that shares drift with Anatolia_N, the three-way model can be rejected outright, or accepted with the Anatolia_N weight pushed to an implausible value. The fix is not a different left set; it is the right set.

Two sources too close. Asking the model to split steppe ancestry between Yamnaya_Samara and a closely related Corded Ware source produces two weights with standard errors of 0.15 or more and Z-scores under 2 on any consumer file. The data cannot tell them apart. Report the single steppe source.

A weight below zero or above one. The nested table's "feasible" column flags these. An infeasible model with a good p-value is not a good model; it usually means a source is standing in for something absent from the left set.

Coverage. On a thin file every standard error sits above 0.10 whatever you do. There is no model-side fix. Say so.

Step 7: reporting the model#

A qpAdm result is only a finding if someone else can reproduce it. That means the record has to carry more than the weights. The minimum:

Every one of those fields, and what each number means, is read line by line in The model record explained. Our own reports publish all of them in the Reading view and download as a plain-text file at every tier. If a product prints percentages with a qpAdm label and none of the fields above, you are not looking at a qpAdm result you can check.

If you would rather not do this yourself#

Everything above is what a qpAdm analysis from Ancestrify does on your behalf, with two differences from an evening in R. The merge runs on our infrastructure against AADR v66, and the model search is done by a person who then has to clear the same bar for every source in every era: p > 0.05, every |Z| > 3, every SE < 0.10. The search covers every proxy label for every candidate source, several hundred hand-built runs, whichever of the four tiers (29.99 to 59.99 EUR, one-time) you choose; the tiers change how the search is presented, not the bar. The report publishes the full record above, the complete right set and the analyst's written explanation of why the genome resolved the way it did.

Afterwards, the Model Lab (10 EUR, one unlock) puts your own merged sample in the workbench: you choose the sources and outgroups, qpadm() runs on your genome, up to 100 runs per rolling 24 hours, and the exact EIGENSTRAT bundle downloads for the R workflow above. The walkthrough is Run your own qpAdm models. Whole-genome VCFs are accepted through a 10 EUR upload add-on. And if you want to see what a finished report looks like before deciding, the demo mirrors a real one. The starting point either way is /qpadm.

References#

• Haak, W. et al. (2015). Massive migration from the steppe was a source for Indo-European languages in Europe. Nature, 522, 207 to 211. Supplementary Information 10 introduces qpAdm.
• Harney, E., Patterson, N., Reich, D. and Wakeley, J. (2021). Assessing the performance of qpAdm: a statistical tool for studying population admixture. Genetics, 217(4), iyaa045.
• Maier, R. et al. (2023). On the limits of fitting complex models of population history to f-statistics. eLife, 12, e85492.
• Mallick, S. et al. (2024). The Allen Ancient DNA Resource (AADR): a curated compendium of ancient human genomes. Scientific Data, 11, 182.
• Patterson, N. et al. (2012). Ancient admixture in human history. Genetics, 192(3), 1065 to 1093.
• Lazaridis, I. et al. (2016). Genomic insights into the origin of farming in the ancient Near East. Nature, 536, 419 to 424.