Applying f4-statistics and admixture graphs: Theory and examples

f-statistics (Patterson et al., 2012; Reich, Thangaraj, Patterson, Price, & Singh, 2009) are a widely used toolkit for making inferences about phylogeny and admixture from population genetic data, particularly in humans. The statistics measure correlations in allele frequencies among sets of two, three or four populations. Observed values reflect degrees of shared ancestry and can serve as a means for testing hypotheses regarding population split orders and past gene flow events under historical models. As compared to some other common methods in population genetics, f-statistics are quite simple and flexible, but interpreting them is not always straightforward. Additionally, one of the primary applications of f-statistics is in building admixture graphs (i.e., phylogenetic trees augmented with admixture events) with more than four populations, which introduces a greater level of complexity. In this note, I hope to clarify some of these potential difficulties and provide a range of tips for practitioners. Some of the topics have been addressed previously but are covered here as well for the sake of completeness. 2 | f-STATISTIC S AND ADMIX TURE

. For example, adding fixed sites to the SNP set will shrink f-statistics towards zero. As a result, when comparing multiple f-statistics, it is important that each one should be computed on the same set of SNPs (or as similar as possible). In applications involving ancient DNA, where missing data is common, I typically make the assumption that the SNPs covered for each individual or population are a random subset with respect to allele frequency. By contrast, comparisons across different genotyping arrays are likely to be biased.

| Interpreting nonzero f 4 -statistics
If a set of four populations are unadmixed relative to each other, then some permutation of them will yield an f 4 -statistic of zero (in expectation), as in Figure 1a. Equivalently, if all three permutations of f 4 -statistics for a certain set of four populations are (significantly) nonzero, then at least one of the populations must be admixed; this is one of the most common signals of admixture used in the literature. In this paper, I will use the example of a quartet consisting of four present-day human populations: Mixe (from Mexico), Han Chinese, French and Baka (hunter-gatherers from Cameroon). The common ancestral population of all Native Americans is known to have been admixed with ~70% ancestry from an eastern Eurasian lineage and 30% from a western Eurasian lineage (Figure 2; Raghavan et al., 2014). Thus, in the context of this quartet, Mixe can be modelled as admixed with ancestry related to Han (~70%) and to French (~30%). I computed the three possible f 4 -statistics for the quartet and obtained significantly nonzero values, with the signs as expected based on the known history (Table 1). (These and all results in the paper are computed from previously published whole-genome sequence data (Fan et al., 2019;Mallick et al., 2016), on a set of ~1.1 million autosomal SNPs (Mathieson et al., 2015), using the implementation in admixtools (Patterson et al., 2012), including standard errors estimated by block jackknife.) In this case, there is prior knowledge available about the admixture in Mixe, but in general, without additional information, the existence of such a quartet does not identify which of the four populations is admixed. Here, for example, it could also be that Han is admixed with most of its ancestry related to Mixe but a small amount related to Baka, and likewise for the other two (see further discussion in the admixture graph sections below). In real-world applications, it can also be true that more than one population is admixed, making the interpretation more complicated. Sometimes, in fact, two admixture events together can cause an f 4 -statistic to be close to zero and thereby mask the signal of admixture (at first glance).
Another observation is that as depicted in Figure 1, f 4 -statistics are not only zero or nonzero but also carry quantitative information about amounts of shared drift between populations. One implication is that populations sharing more drift (i.e., yielding longer  (c) With population C admixed, the path from B to C can be decomposed into two components. Under the model shown, with a proportion of B-related ancestry and 1 − D-related ancestry, the former yields a path (lighter red) that has a weight of but does not intersect the path from A to D, while the latter yields a path (darker red) that has a weight of 1 − and intersects the path from A to D over the branch with length y.
intersecting paths in an admixture graph) will have greater-magnitude f 4 -statistics associated with them. For example, in the trees of C, D′) to conclude that D is a better proxy than D′ for the ancestry in C (the component with proportion 1 − ). However, this procedure is complicated by the fact that if the D-related source was in fact itself admixed, with ancestry related to X and Y, then the f 4 -statistic can sometimes be maximized by X or Y instead of by D, even though one would consider D to be a better proxy (Pickrell et al., 2014). It is also good to remember that if a certain signal is weak compared to the noise in the data-for example, if one were testing for admixture in C and the shared drift branch length y was short-then one may not have enough power to identify it.
Finally, f-statistics can be subject to certain kinds of biases and batch effects (to varying degrees, as with other methods) arising from SNP ascertainment, sample type and processing (ancient vs.
present-day, sequencing platform, etc.), and other aspects of the data, so it is important to keep such factors in mind when interpreting results. For ancient DNA data, challenges include C-to-T errors induced by postmortem deamination (Hofreiter, Jaenicke, Serre, Haeseler, & Pääbo, 2001), as well as short fragment lengths and (often) low coverage, which can exacerbate reference bias (Günther & Nettelblad, 2019). All of these effects can cause ancient individuals to appear artificially closely related to one another and to certain other populations (e.g., deep outgroups). In general, statistics f 4 (A, B; C, D) in which A and C share a data type and B and D share a different data type are most prone to this kind of artefact.

| Fitting an admixture graph with qpGraph
In addition to their stand-alone usage, f-statistics can serve as a means to fit admixture graphs from allele frequency data. (Other kinds of statistics can also be used to fit admixture graphs, but I will not discuss such methods in detail here; see Discussion.) In this con- The software I typically use to build admixture graphs is qpGraph (also referred to as admixturegraph; Patterson et al., 2012). In qp-Graph, the user manually specifies the topology of the model, and the program then solves for the optimal values of the parameters.
In theory, one might wish to search the entire space of all topologies and parameter values (for a given number of admixture events) to find the best-fitting model, but the size of the space (exponential in the number of populations) makes this impractical for larger graphs (Leppälä, Nielsen, & Mailund, 2017). The set of basis statistics used for fitting is the set (2) alluded to in the previous paragraph, with the first population listed in the input file as the 'base' population.
In its standard mode, qpGraph attempts to minimize the quantity vector of predicted f-statistics under the model, and Q is the (estimated) covariance matrix of the statistics. Assuming multivariate normal errors, the score gives the negative log-likelihood of the model; it measures the total amount by which the system of f-statistic equations (one for each basis statistic) fails to be satisfied, taking into account the empirical correlation among the statistics (see also the next section on fit quality). To help insure that Q −1 does not become unstable, one can use the 'diag' input parameter to add a small number ('diag: 0.0001' works well in my experience, but smaller values may be sufficient as well) to the diagonal entries of Q. The program can also be run using simple least-squares optimization without the Q matrix by specifying 'lsqmode: YES', but in this case highly correlated statistics will be treated as independent for the sake of the fitting, and the score will no longer represent a log-likelihood, both of which make the full objective function preferable. Other input parameters I typically set are 'outpop: NULL' (meaning no specified outgroup population in which SNPs are required to be polymorphic) and 'lambdascale: 1' (leaving the f-statistics in typical units rather than scaling into approx-

| Parameters and constraints
An important consideration is whether the system of equations used to infer the parameters of an admixture graph is over-or under-determined. As mentioned above, a model with n populations linearly independent constraints (i.e., equations). In the absence of admixture, there are 2n − 3 parameters, which is the number of branches in an unrooted binary tree with n leaf nodes (with the settings I have described, qpGraph results should not depend on where the root of a graph is specified). Converting a population from unadmixed to admixed adds two parameters: one for the mixture proportion and one for the split position of the new source of ancestry. Thus, with a admixture events, the total number of free parameters is 2n + 2a − 3. One point to note is that in the case of an admixed population with two unsampled sources (which is the typical scenario), the three branch lengths surrounding the admixture event (in Figure 3a, from the node 'East1' to 'East2', from 'West1' to 'West2' and from 'pAM1' to Mixe) cannot be determined individually but instead form a single compound parameter 2 x + (1 − ) 2 y + z (where is the mixture proportion, x and y are the branch lengths to the two corresponding sources, and z is the terminal branch length

| Fit quality
To my knowledge, no absolute measure of model fit has been de- Deviations between model predictions and the observed data can be caused either by an incorrectly specified topology or un-modelled admixture. In the first case, assuming that the program does not get stuck at a local optimum, it will try to move the populations as close as possible to their correct positions but will be constrained by the input topology. Thus, an incorrectly specified split order usually manifests as an inferred length-zero internal branch; when such branches (i.e., trifurcations) appear in the results, the order of splits should be adjusted and re-tried. (The default qpGraph visualization output rounds branch lengths to the nearest integer, so some nonzero length but very short branches may initially appear as zero.) As noted in the f-statistics section above, however, one may not have sufficient power to resolve short branches, so some sets of three lineages may be found to be statistically consistent with forming a trifurcation, with all three possible split orders having similar fit quality.
In the case of un-modelled admixture, the observed deviations The score of the final graph is also returned as an output from the program, so it can be used to compare the fit quality of different models with the same set of populations, preferring the one with the lower score. (If the equations being fit were independent, then one could apply a chi-squared test for the overall fit, but in practice they are heavily correlated. qpGraph returns a naive degrees of freedom count and p-value alongside the score, but they are not well calibrated.) As above, while this approach provides a useful heuristic, evaluating statistical significance is complicated, and I do not have a rigorous set of recommendations. One recent direction that seems promising is using the score to compare alternative models with the same populations and same number of admixture events. In that case, the score difference can be interpreted in an AIC/BIC framework, with the likelihood difference as a Bayes factor (Flegontov et al., 2019;Leppälä et al., 2017;Shinde et al., 2019). The same idea could also be applied in cases with unequal numbers of free parameters-for example, adding one admixture event and testing whether the score improvement is significant. However, defining the change in degrees of freedom is not straightforward in this situation: as noted above, a new admixture event creates two additional parameters in the model, but that does not account for whether the admixture comes from a prespecified source or from a source that is allowed to be located anywhere in the graph. Finally, the score can additionally be used to compute confidence intervals on parameters (by considering the likelihood as a function of a single branch length or mixture proportion value), although it is worth keeping in mind that the results are model-dependent.

| ADMIX TURE G R APH S: E X AMPLE S
One of the strengths of f-statistic-based admixture graphs is that they are computationally tractable enough that programs such as

| Four populations
The first examples I will present are four-population admixture graphs containing Mixe, Han, French and Baka. Given the observed nonzero f 4 -statistics in Table 1, there must be at least one admixture event present in order to fit the data. However, in light of the discussions above about determining which population is admixed and about parameters and constraints in admixture graphs, it would be expected that these models should be insufficiently constrained to determine which population is admixed. Indeed, they have constraints but 2 (4) + 2 (1) − 3 = 7 free parameters. Confirming this expectation, perfectly fitting models (i.e., sets of branch length and mixture proportion parameters such that the six basis f-statistics are predicted exactly, yielding S (G) = 0) can be obtained with Mixe specified as admixed (Figure 3a) as well as with any of the other three populations (incorrectly) specified as admixed instead (Figure 3b-d).
Interestingly, in some scenarios, the admixed population can

| Five populations
In general, in order to be able to solve for the parameters of an admixture graph including one admixture event, it is necessary to use at least five populations, providing = 10 constraints for the 2n + 2a − 3 = 9 free parameters. Concurrently, in contrast to the fourpopulation examples above, having five populations present allows one to determine which of the populations is admixed, as long as the topological relationships of the populations are all unique relative to the true mixing sources. More detail on this last point can be found elsewhere Pease & Hahn, 2015). A simple version of this statement is that, at least in the case of a single admixture event, one four-population subset will be unadmixed, whereas the other four subsets will include the admixed population. Similarly, in order to solve for a given mixture proportion in a larger graph, there must four populations present (aside from the admixed one in question) in distinct positions, yielding a nonredundant five-population subgraph; having three populations in distinct positions allows one to detect the signal of admixture but not to determine the proportion uniquely.
As an example, I added Ulchi (from the Amur River Basin of northeastern Asia) as a fifth population alongside the four from above.
Ulchi splits closer to the eastern Eurasian source population for Mixe than does Han, which provides the additional degree of constraint.
The five-population model is a good fit to the data, but not a perfect one (Z = 1.9 for the most significant residual; Figure 5a). By contrast, if Baka are modelled as admixed instead of Mixe, the fit is poor (Z = 4.7 ; Figure 5b). I also show an example where the topology is incorrectly specified, with Han closer than Ulchi to the eastern Eurasian source population for Mixe ( Figure 5c); this version fits poorly (Z = 5.7), and the branch connecting the split positions of Ulchi and Han collapses to length zero. If I add a second admixture event into the models in Figure 5a,b, this creates more free parameters (11) than constraints, and indeed, there are choices of the parameters that yield perfect fits, even with Mixe modelled as unadmixed (not shown).
Having five populations present (with a single admixture event) also provides the ability to infer uniquely optimal parameter values.
In the four-population example model, the initial estimate of eastern Eurasian ancestry in Mixe was 71%, but with the proportion manually set at 75%, the fit is still perfect (Figure 6a). Outside of a certain range of mixture proportions (dependent on the values of the branch lengths), the fit will become worse, but within a finite interval, the likelihood is entirely flat. In terms of f 4 -statistics, the observed nonzero value is being fit as equal to a branch length in the admixture graph times the mixture proportion (as in Figure 1c) Most of the material in this paper pertaining to admixture graphs has been presented from the perspective of the qpGraph software, but other methods are also available, using both different kinds of data and different fitting schemes. At the level of mathematical formulation, the results have assumed that models are fit based on a distance metric (specifically, f-statistics). As an alternative example, the treemix algorithm (Pickrell & Pritchard, 2012) is based on a maximum-likelihood framework in terms of allele frequency covariances, although the information captured is the same; see Peter (2016) for the equivalence and a thorough exploration of alternative interpretations of f-statistics in terms of population genetic models. There are also methods that use richer summaries of the data (e.g., the full joint allele frequency spectrum) to infer more complicated demographic models that are similar in form, or in some cases essentially identical, Within the class of f-statistic-based (or equivalent) admixture graph methods, there are different approaches to automation and the selection of which populations to model as admixed. qpGraph leaves the choice of how many admixture events to include (and which populations are admixed) up to the user; some guidelines pertaining to this choice have been discussed above. For smaller models, it can also be possible to search some or all of the full graph space (Shinde et al., 2019) to determine best-fitting topologies for a given number of admixture events (e.g., using the similar admixturegraph r implementation (Leppälä et al., 2017) and admixturebayes (Nielsen, 2018); other techniques are the subject of ongoing work).
mixmapper (Lipson et al., 2013) provides an intermediate level of automation by attempting to infer an unadmixed submodel and then fitting one or two admixed populations onto this scaffold. With a small set of populations, this can sometimes be a useful approach, but it can largely be recapitulated within qpGraph, and the software does not support large models with more admixture events.
At the most automated end of the spectrum is treemix (Pickrell & Pritchard, 2012), which only asks the user to supply the list of populations and the number of admixture events and then returns a single inferred model. The advantage of this strategy is that the program does all of the work of building the graph, which is especially useful if one has limited prior knowledge about the populations. The main drawback, in my view, is that the way the program builds the graph is by starting with an optimal mixture-free tree and then adding admixture events to account for deviations between the predictions of the tree model and the observed data. Depending on the true histories of the populations, this approach can be successful, but it can also increase the chances of falling into local optima imposed by the initial tree (especially if many populations are admixed; see Lipson et al., 2013). Additionally-as in other methods-the choice of how many admixture events to include, which can sometimes be difficult, is still left to the user.
In my experience, I have found f-statistics and admixture graphs to be very useful tools for learning about phylogeny and admixture.
I hope that this guide will help others to get the most out of these tools in a range of real-world applications.

ACK N OWLED G M ENTS
I would like to thank David Reich, Vagheesh Narasimhan, Nick Patterson, Robert Maier, Iosif Lazaridis and Pavel Flegontov for helpful discussions, and three anonymous reviewers for comments.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available through the European Nucleotide Archive (ENA), under accession numbers PRJEB9586 and ERP010710, and at the European Genome-phenome Archive (EGA), under accession number EGAS00001001959 (Fan et al., 2019;Mallick et al., 2016).