Author response: On the limits of fitting complex models of population history to f-statistics
Year: 2023
Type: peer-review
Abstract: Full text Figures and data Side by side Abstract Editor's evaluation Introduction Results Discussion Appendix 1 Appendix 2 Data availability References Decision letter Author response Article and author information Metrics Abstract Our understanding of population history in deep time has been assisted by fitting admixture graphs (AGs) to data: models that specify the ordering of population splits and mixtures, which along with the amount of genetic drift and the proportions of mixture, is the only information needed to predict the patterns of allele frequency correlation among populations. The space of possible AGs relating populations is vast, and thus most published studies have identified fitting AGs through a manual process driven by prior hypotheses, leaving the majority of alternative models unexplored. Here, we develop a method for systematically searching the space of all AGs that can incorporate non-genetic information in the form of topology constraints. We implement this findGraphs tool within a software package, ADMIXTOOLS 2, which is a reimplementation of the ADMIXTOOLS software with new features and large performance gains. We apply this methodology to identify alternative models to AGs that played key roles in eight publications and find that in nearly all cases many alternative models fit nominally or significantly better than the published one. Our results suggest that strong claims about population history from AGs should only be made when all well-fitting and temporally plausible models share common topological features. Our re-evaluation of published data also provides insight into the population histories of humans, dogs, and horses, identifying features that are stable across the models we explored, as well as scenarios of populations relationships that differ in important ways from models that have been highlighted in the literature. Editor's evaluation This is a rigorous and critical analysis of the performance of a popular suite of methods for inferring population history, accompanied by improvements. Should be of broad interest to anyone interested in human history. https://doi.org/10.7554/eLife.85492.sa0 Decision letter eLife's review process Introduction Admixture graph models provide a powerful intellectual framework for describing the relationships among populations that allows not only branching of populations from a common ancestor but also mixture events. An admixture graph (abbreviated below as AG), as fit in the widely used software packages ADMIXTOOLS (Patterson et al., 2012) and TreeMix (Pickrell and Pritchard, 2012; Molloy et al., 2021), is a directed acyclic bifurcating graph with two types of edges: those representing genetic drift, and those representing gene flow. Each admixture event is represented as a confluence of two gene flow edges. Nodes of such a graph represent unsampled intermediate populations, and terminal nodes (leaves) represent sampled present-day or ancient groups (see a mathematical definition in Soraggi and Wiuf, 2019). An attractive feature of AGs is that they can summarize important features of population history without requiring specification of all parameters such as population sizes, split times, mixture times, and distinguishing between sudden splits or drawn-out separations. All these parameters describe important features of demographic history and are fit by many methods for fitting demographic models (Gutenkunst et al., 2009; Gronau et al., 2011; Schiffels et al., 2016; Flegontov et al., 2019; Kamm et al., 2020, Rogers, 2019; Hubisz et al., 2020). However, the fact that it is possible to factor this difficult problem by first inferring important aspects of the topology (AGs fitted to allele frequency correlation statistics), and then fitting additional demographic parameters to data such as site frequency spectra, simplifies demographic inference (Patterson et al., 2012; Pickrell and Pritchard, 2012; Lipson et al., 2013; Leppälä et al., 2017; Lipson et al., 2020b; Molloy et al., 2021; Yan et al., 2021). AGs thus serve both as conceptual frameworks that allow us to think about the relationships of populations deep in time, and as mathematical models we can fit to genetic data. AGs are fitted to f-statistics (Reich et al., 2009; Patterson et al., 2012; Peter, 2016; Soraggi and Wiuf, 2019). For convenience, below we use a concise definition of f-statistics by Lipson, 2020a: ‘The most general definition is that of the f4-statistic f4(A, B; C, D), which measures the average correlation in allele frequency differences between (1) populations A and B and (2) populations C and D that is, (pA – pB)∗ (pC – pD), for allele frequencies p, typically averaged over many biallelic single-nucleotide polymorphisms. This f4-statistic is the same as the D-statistic up to a normalization factor.’ The other f-statistics (f2 and f3) can be defined as special cases of f4-statistics: f2(A, B) = f4(A, B; A, B) and f3(A; B, C) = f4(A, B; A, C). f4-Statistics can be written as linear combinations of f3- or f2-statistics, and f3-statistics can be written as linear combinations of f4- and f2-statistics. f2-, f3-, and f4-statistics have straightforward interpretations in terms of drift edges along the tree, see Figure 2 in Patterson et al., 2012 and Appendix 1—figure 1b. A challenge for fitting AG models is that they are often not uniquely constrained by the data, with many providing equally good fits to the f2-, f3-, and f4-statistics used to constrain them within the limits of statistical resolution. Previously published methods for finding fitting AGs (mainly qpGraph, Patterson et al., 2012 and TreeMix, Pickrell and Pritchard, 2012; Molloy et al., 2021) were not well equipped to handle the large range of equally well-fitting models for three reasons: (1) They did not reliably provide information on whether there is a uniquely fitting parsimonious model or alternatively whether there are many models that fit equally well to the limits of statistical resolution, (2) they did not provide formal goodness-of-fit tests, and related to this, (3) they did not provide tests for whether the difference between the fits of any two models is statistically significant. As a consequence, and as we demonstrate in what follows, many published AG models have been interpreted as providing more confidence than is merited about the extent to which genetic data allows us to disentangle ancestral relationships. To appreciate these problems, we first need to consider two main approaches that were utilized to study demographic history with AGs. The first approach is to identify AGs automatically, either without human intervention or with guidance. It is possible in theory to exhaustively test all possible graphs for a given set of populations and pre-specified number of admixture events, as implemented, for example, in the admixturegraph R package (Leppälä et al., 2017). An exhaustive approach can provide a complete view of the range of models that are consistent with the data for a specified level of parsimony (total number of admixture events allowed in the graph), which is not biased by the algorithm used to explore the space of possible AGs. However, this approach is limited to small graphs (typically up to six groups, two admixture events) due to the rapid increase in the number of possible AGs as the number of populations and admixture events grows. As we show in our discussion of case studies, the simple models explored with an exhaustive approach can lead to misleading conclusions about population history because not including additional populations can blind users to additional mixture events that occurred (and whose existence is revealed by examining data from additional populations). Furthermore, models with additional admixture events that are qualitatively different to the best-fitting parsimonious graph and that capture the true history, will sometimes be completely missed when constraining the number of gene flows. Alternatively, the programs TreeMix (Pickrell and Pritchard, 2012; Molloy et al., 2021), MixMapper (Lipson et al., 2013), miqoGraph (Yan et al., 2021), and AdmixtureBayes (Nielsen et al., 2023 preprint) all address the problem of how to rapidly explore the vast space of AGs relating a set of populations by applying algorithmic ideas or heuristics; all of these methods speed up model search by orders of magnitude. The second approach to fitting AGs is to manually build them up by grafting additional populations onto simpler smaller graphs that fit the data. This approach involves stepwise addition of populations in an order that is chosen based on the best judgment of the user, and for each newly added population involves adding admixture events or tweaks in the graph until a fit is obtained; the user then moves on to adding the next population (see Reich et al., 2009; Reich et al., 2011; Reich et al., 2012; Lazaridis et al., 2014; Seguin-Orlando et al., 2014; Fu et al., 2016; Skoglund et al., 2016; Yang et al., 2017; McColl et al., 2018; Moreno-Mayar et al., 2018; Tambets et al., 2018; van de Loosdrecht et al., 2018; Flegontov et al., 2019; Sikora et al., 2019; Wang et al., 2019; Lipson et al., 2020b; Shinde et al., 2019; Yang et al., 2020; Hajdinjak et al., 2021; Wang et al., 2021; Bergström et al., 2022 for examples). The program qpGraph in the ADMIXTOOLS package (Patterson et al., 2012) has been the most common computational method used for testing fits of individual AGs. Most AGs in the literature have been constructed manually in this way, often acknowledging the existence of alternative models by presenting plausible models side-by-side, and this approach has been the basis for many claims about population history (Lazaridis et al., 2014; Yang et al., 2017; Posth et al., 2018; Sikora et al., 2019; Shinde et al., 2019; Bergström et al., 2020; Lipson et al., 2020b; Hajdinjak et al., 2021; Wang et al., 2021; Bergström et al., 2022). A strength of this approach is that it takes advantage of human judgment and outside knowledge about what graphs best fit the history of the human or animal populations being analyzed. This external information is powerful as it can incorporate non-genetic evidence such as geographic plausibility and temporal ordering of populations or linguistic similarity, or other genetic data such as estimates of population split times, or shared Y chromosomes, or rejection of proposed scenarios based on joint analysis of much larger numbers of populations than can reasonably be analyzed within a single AG. Thus, while manual approaches explore many orders of magnitude fewer topologies than automatic approaches often do, they still may provide inferences about population history that are more useful than those provided by automatic approaches. These methods’ strength is also their weakness: by relying on intuition, following a manual approach has the potential to validate the biases users have as to what types of histories are most plausible (these may be the only types of histories that will be carefully explored). This can blind users to surprises: to profoundly different topologies that may correspond more closely to the true history, and we discuss examples of this in the Results section. In this study, we introduce a new method, findGraphs, that belongs to the first class of algorithms (those for automated AG topology inference). Algorithmic innovations and speedups in findGraphs enable us to explore a much larger proportion of plausible AG space than many other methods reported to date. The findGraphs method combines the advantages of automated and manual topology exploration by allowing users to encode various sources of information as constraints on the space of AGs, which is then explored automatically. However, the main innovations in findGraphs are not computational, but instead conceptual. Instead of finding one or a few AGs fitting the data well, we use findGraphs for exploring AG spaces and assessing if any reliable information on population history can be extracted from a given AG space (defined by a population set and parsimony constraints) in the first place. Results Regardless of the approach used to search through the space of possibly fitting AGs, a challenge in the effort to find a uniquely well-fitting AG (or group of topologically similar AGs) is that it has been difficult to quantify the absolute goodness of fit of a model to date. We have not been entirely successful with this and are not aware of other work that has been successful. It is also difficult to assess the relative fits of multiple models, especially if they differ in complexity. Performance gains relative to the original implementation of qpGraph allow us to address this problem by obtaining bootstrap confidence intervals and p-values for estimated parameters of single models, as well as for the difference in fit quality of two models (see Appendix 1, Sections 1.B.3 and 2.E). In combination with the approach to automating the search of well-fitting AGs, this leads to a situation where we are able to find and test a large number of models, many of which fit equally well despite often having very different topological features. Published approaches to comparing the fits of AG models based on Akaike information criterion (AIC) or Bayesian information criterion (BIC), see Flegontov et al., 2019; Shinde et al., 2019 have the problem that it is often not clear what the effective number of degrees of freedom is in the two models being compared since in the case of AGs it depends not only on the number of graph edges, but also on graph topology. The methods for automated graph topology inference and model comparison relying on bootstrap resampling are implemented in ADMIXTOOLS 2, a comprehensive platform for learning about population history from f-statistics. It is built to provide a stand-alone workspace for research in this area and is implemented as an R package. For all computations, ADMIXTOOLS 2 exhibits large speedups relative to previously published platforms for f-statistic analysis (e.g., popstats and ADMIXTOOLS version 6.0 which we call ‘Classic ADMIXTOOLS’ in what follows to distinguish it from updated ADMIXTOOLS version 7.0.2 which implements some of the speedup ideas also implemented in ADMIXTOOLS 2). This is achieved by deploying a series of algorithmic improvements, most notably storage of precomputed f-statistics in random access memory, which avoids having to rely on reading in extremely large genotype matrices to perform most computations. In addition to the new algorithmic ideas allowing efficient searching through the space of AGs and comparing the fits of two AGs, ADMIXTOOLS 2 also provides a solution to the question of which parameters of an AG are identifiable in the limit of infinite data. Methodological details are presented in Appendix 1, and below we focus on documenting problems of AG inference on simulated data and revisiting AGs from the literature to understand the extent to which methodological challenges with AG fitting biased previous studies. Topological diversity of well-fitting models and effects of parsimony constraints on simulated data First, we explored the performance of the findGraphs method for automated topology inference on simulated AGs of random topology, focusing on the following questions: (1) among findGraphs results, how common are AGs fitting nominally or significantly better than the true one but different topologically; (2) what is the degree of topological diversity among these models fitting the data better than the true one? For this purpose, we simulated AGs of four complexity classes using msprime v.1.1.1: eight or nine non-outgroup populations, and four or five admixture events. Only simulations where pairwise FST for groups were in the range characteristic for anatomically modern and archaic humans were selected for further analysis, resulting in 20 random topologies per complexity class, each including a distant outgroup that facilitates automated exploration of the topology space. We ran findGraphs on each simulated dataset starting from random graphs and pre-specifying the true number of admixture events (n), or n − 1, or n + 1 events. For each of these graph complexity levels, we performed 100 independent findGraphs runs and recorded 5 AGs from each run having the best log-likelihood (LL) scores. Topologically redundant AGs were discarded, and for the remaining AGs we calculated worst f4-statistic residuals (WR) and tested if the newly found models fit significantly better than the true model, using the bootstrap model comparison method developed in this study (see Appendix 1, Sections 1.B.3 and 2.E). In Figure 1a–c, we show the following statistics for each simulated AG, summarized across simulated complexity classes and parsimony levels allowed at the stage of topology exploration: fraction of topologies found with findGraphs that fit better than the true AG (according to LL score), or that fit significantly better than the true AG, or those with plausible absolute fits (WR < 3 SE). It is clear that for the great majority of simulated datasets, even a shallow exploration of the topology space with findGraph (100 independent runs) uncovers AGs that fit nominally better than the true topology (Figure 1a and d) and are topologically diverse (see Figure 1e for examples). When allowing for n admixture events, at least one AG fitting significantly better than the true one was found for 60% of simulated datasets. When n + 1 admixture events were allowed, this grew to 100% (all 80 datasets). It should be noted that some admixture events are indistinguishable with f-statistics; for instance, successive gene flows between two lineages, with no other edges branching off between the gene flows. If such gene flows were included in the random topologies we simulated, AGs with n events were overly complex for representing the true history. Thus, if we are dealing with random histories, choosing an optimal complexity class for topology search is not straightforward. Figure 1 Download asset Open asset Computer simulations show that when the true admixture graph (AG) topology is complex, findGraphs frequently finds AGs fitting the data better than the true AG. (a) Fractions of distinct AGs found with findGraphs that fit the data nominally better than the true AG (according to log-likelihood [LL] scores). The simulated datasets are grouped by complexity class (eight or nine leaves, four or five admixture events) and by the number of admixture events allowed at the topology search stage (n − 1 on the left, n in the middle, and n + 1 on the right, where n is the true number of simulated admixture events). Each dot represents a simulated random history, and 20 such histories were simulated for each complexity class. (b) Fractions of distinct AGs found with findGraphs that fit the data significantly better than the true AG of the bootstrap model comparison method Fractions of distinct AGs found with findGraphs that fit the data well in absolute terms (WR < 3 SE). AGs found for a simulated history (eight groups and four admixture events) in the LL and Only best-fitting graphs with < 3 are The fit of the true topology is in and topologies that fit the data significantly better than the true one are in The true topology was not by our findGraphs The true model from and two alternative models found with findGraphs, both fitting significantly better than the true one on the bootstrap and very different This is presented as an of very topological diversity among well-fitting parameters that were to be (see Appendix 1, are in These results on simulated data about the extent to which fitting AG topologies provide reliable information about population history. for histories including eight or nine groups, an and four or five admixture events, data, and groups as as and anatomically modern complexity class that is simpler than many models fitted to genetic data in published fitting the data as well as or better than the true one are and their topological diversity is in most cases that it inference of topology by analysis of multiple As we demonstrate on set of simulated data in Appendix 1 the of finding a model that fits better than the true one with graph and that is with both findGraphs and TreeMix 1—figure 2). We this problem to be even more when are dealing with complex However, often rely on external constraints on AG topologies as temporal plausibility of a topology, results from and linguistic that were not used for the results of our topology on simulated data. Thus, it is possible in that published AG models are more than our results on simulated data, and we explore this in in the next and in Appendix published AGs We AGs from eight publications (Lazaridis et al., 2014; Shinde et al., 2019; Sikora et al., 2019; Bergström et al., 2020; Lipson et al., 2020b; Hajdinjak et al., et al., 2021; Wang et al., 2021) with the of comparing published models to models identified by our algorithm for inferring optimal AGs In all but one study, qpGraph or automated reimplementation Leppälä et al., was used for fitting topologies to genetic data, while et al., on the automated method et al., 2021). The main question we were interested in is whether we can find alternative models which (1) fit as well or better than the published (2) differ in important ways from the published and (3) be based on other evidence such as temporal The studies were selected to the criterion that an AG model in the study is used as evidence for at least one about population history in the main text of the In other the AG method was used in the original studies to new conclusions about population history, and not to show that there is a model that that not results of other genetic an approach that is a use of AGs and has been in some studies (e.g., Seguin-Orlando et al., 2014; et al., 2019; Wang et al., 2019). are many published studies that have been included in our re-evaluation as they our key criterion (e.g., Yang et al., 2017; McColl et al., 2018; Posth et al., 2018; Flegontov et al., 2019; et al., et al., 2021; Bergström et al., Lipson et al., et al., 2022). However, critical re-evaluation of each published graph is an and the of studies we is diverse to identify some general 1 Published graphs in the of found We compared graphs from eight publications to alternative graphs on the same or very similar data (see 1 for in the original worst alternative topologies better fitting better fitting fitting fitting et al., et al., et al., et al., et al., et al., et al., et al., of the first author and of the Figure in the original Figure number in the original where the AG is The number of populations in each Admixture The number of admixture events in each The number of single-nucleotide with no data at the group used for fitting the AGs. For all case studies, we tested the original data population and the published graph and model fits very similar to the published However, for the of efficient topology we in some cases for population or graph complexity as noted in the in 1, and in the The worst f-statistic of the published graph fitted to the set in the in alternative topologies The number of distinct newly found topologies from the published one. better fitting The of distinct alternative topologies that fit significantly better than the published graph to the bootstrap model comparison test If the number of distinct topologies was very a of models to of models along the log-likelihood was compared to the published one and the in this and following were calculated on this better fitting The of distinct topologies that fit better than the published graph to the bootstrap model comparison test fitting The of distinct topologies that fit than the published graph to the bootstrap model comparison test fitting The of distinct topologies that fit significantly than the published graph to the bootstrap model comparison test The population was see 1 and the gene flows were from the published model for see 1 and the we a of these Discussion of individual case studies follows and for details see the in Appendix For of published graphs we we were able to find at least but graphs of the same complexity of groups and admixture with an LL that was nominally better than that of the published graph (see results for selected graphs in 1 and results for all graphs in The graphs were from the publications as there were multiple graphs presented in some of the publications et al., 2019; Sikora et al., 2019; et al., 2021), or we selected in the model process et al., 2020; Lazaridis et al., 2014; Lipson et al., 2020b; Wang et al., 2021), or we an outgroup not used in the original study et al., 2021; Sikora et al., or we tested additional graph complexity classes admixture events (Lipson et al., 2020b; Sikora et al., 2019). These alternative graphs often fit not significantly better than the published one into across single-nucleotide In the following at least one model that fits significantly better than the published one to our bootstrap model comparison method was the Bergström et and Lazaridis et the et al., graph with 3 admixture the Hajdinjak et graphs with or without adding a the Lipson et al., intermediate graphs with groups and admixture events and with groups and admixture the Wang et and the Sikora et graphs for and for with or admixture events In nearly all cases for the Lazaridis et Shinde et graph with eight populations and three admixture events, and the et al., graph with four admixture we also identified many additional graphs that fit the data not significantly than the published In example, some of these graphs have topologies that are qualitatively different in important ways from those of the published such as which populations are or of gene or the order of split events, if not constrained a are not the same between alternative fit
Source: eLife
Language: English
Subfield: Genetics
Domain: Life Sciences
Open Access status: gold
Estimated APC paid: $2,500
OpenAlex ID: W4382724322