1. Computational and Systems Biology
  2. Evolutionary Biology
Download icon

Fast and flexible estimation of effective migration surfaces

  1. Joseph Marcus  Is a corresponding author
  2. Wooseok Ha  Is a corresponding author
  3. Rina Foygel Barber  Is a corresponding author
  4. John Novembre  Is a corresponding author
  1. Department of Human Genetics, University of Chicago, United States
  2. Department of Statistics, University of California, Berkeley, United States
  3. Department of Statistics, University of Chicago, United States
  4. Department of Ecology and Evolution, University of Chicago, United States
Tools and Resources
  • Cited 0
  • Views 384
  • Annotations
Cite this article as: eLife 2021;10:e61927 doi: 10.7554/eLife.61927

Abstract

Spatial population genetic data often exhibits ‘isolation-by-distance,’ where genetic similarity tends to decrease as individuals become more geographically distant. The rate at which genetic similarity decays with distance is often spatially heterogeneous due to variable population processes like genetic drift, gene flow, and natural selection. Petkova et al., 2016 developed a statistical method called Estimating Effective Migration Surfaces (EEMS) for visualizing spatially heterogeneous isolation-by-distance on a geographic map. While EEMS is a powerful tool for depicting spatial population structure, it can suffer from slow runtimes. Here, we develop a related method called Fast Estimation of Effective Migration Surfaces (FEEMS). FEEMS uses a Gaussian Markov Random Field model in a penalized likelihood framework that allows for efficient optimization and output of effective migration surfaces. Further, the efficient optimization facilitates the inference of migration parameters per edge in the graph, rather than per node (as in EEMS). With simulations, we show conditions under which FEEMS can accurately recover effective migration surfaces with complex gene-flow histories, including those with anisotropy. We apply FEEMS to population genetic data from North American gray wolves and show it performs favorably in comparison to EEMS, with solutions obtained orders of magnitude faster. Overall, FEEMS expands the ability of users to quickly visualize and interpret spatial structure in their data.

Introduction

The relationship between geography and genetics has had enduring importance in evolutionary biology (see Felsenstein, 1982). One fundamental consideration is that individuals who live near one another tend to be more genetically similar than those who live far apart (Wright, 1943; Wright, 1946; Malécot, 1948; Kimura, 1953; Kimura and Weiss, 1964). This phenomenon is often referred to as ‘isolation-by-distance’ (IBD) and has been shown to be a pervasive feature in spatial population genetic data across many species (Slatkin, 1985; Dobzhansky and Wright, 1943; Meirmans, 2012). Statistical methods that use both measures of genetic variation and geographic coordinates to understand patterns of IBD have been widely applied (Bradburd and Ralph, 2019; Battey et al., 2020). One major challenge in these approaches is that the relationship between geography and genetics can be complex. Particularly, geographic features can influence migration in localized regions leading to spatially heterogeneous patterns of IBD (Bradburd and Ralph, 2019).

Multiple approaches have been introduced to model spatially non-homogeneous IBD in population genetic data (McRae, 2006; Duforet‐Frebourg and Blum, 2014; Hanks and Hooten, 2013; Petkova et al., 2016; Bradburd et al., 2018; Al-Asadi et al., 2019; Safner et al., 2011; Ringbauer et al., 2018). Particularly relevant to our proposed approach is the work of Petkova et al., 2016 and Hanks and Hooten, 2013. Both approaches model genetic distance using the ‘resistance distance’ on a weighted graph. This distance metric is inspired by concepts of effective resistance in circuit theory models, or alternatively understood as the commute time of a random walk on a weighted graph or as a Gaussian graphical model (specifically a conditional auto-regressive process) (Chandra et al., 1996; Hanks and Hooten, 2013; Rue and Held, 2005). Additionally, the resistance distance approach is a computationally convenient and accurate approximation to spatial coalescent models (McRae, 2006), although it has limitations in asymmetric migration settings (Lundgren and Ralph, 2019).

Hanks and Hooten, 2013 introduced a Bayesian model that uses measured ecological covariates, such as elevation, to predict genetic distances across sub-populations. Specifically, they use a graph-based model for genotypes observed at different spatial locations. Expected genetic distances across sub-populations in their model are given by resistance distances computed from the edge weights. They parameterize the edge weights of the graph to be a function of known biogeographic covariates, linking local geographic features to genetic variation across the landscape.

Concurrently, the Estimating Effective Migration Surfaces (EEMS) method was developed to help interpret and visualize non-homogeneous gene-flow on a geographic map (Petkova, 2013; Petkova et al., 2016). EEMS uses resistance distances to approximate the between-sub-population component of pairwise coalescent times in a ‘stepping-stone’ model of migration and genetic drift (Kimura, 1953; Kimura and Weiss, 1964). EEMS models the within-sub-population component of pairwise coalescent times, with a node-specific parameter. Instead of using known biogeographic covariates to connect geographic features to genetic variation as in Hanks and Hooten, 2013, EEMS infers a set of edge weights (and diversity parameters) that explain the genetic distance data. The inference is based on a hierarchical Bayesian model and a Voronoi-tessellation-based prior to encourage piece-wise constant spatial smoothness in the fitted edge weights.

EEMS uses Markov Chain Monte Carlo (MCMC) and outputs a visualization of the posterior mean for effective migration and a measure of genetic diversity for every spatial position of the focal habitat. Regions with relatively low effective migration can be interpreted to have reduced gene-flow over time, whereas regions with relatively high migration can be interpreted as having elevated gene-flow. EEMS has been applied to multiple systems to describe spatial genetic structure, but despite EEMS’s advances in formulating a tractable solution to investigate spatial heterogeneity in IBD, the MCMC algorithm it uses can be slow to converge, in some cases leading to days of computation time for large datasets (Peter et al., 2020).

The inference problems faced by EEMS and Hanks and Hooten are related to a growing area referred to as ‘graph learning’ (Dong et al., 2019; Mateos et al., 2019). In graph learning, a noisy signal is measured as a scalar value at a set of nodes from the graph, and the aim is then to infer non-negative edge weights that reflect how spatially ‘smooth’ the signal is with respect to the graph topology (Kalofolias, 2016). In population genetic settings, this scalar could be an allele frequency measured at locations in a discrete spatial habitat with effective migration rates between sub-populations. Like the approach taken by Hanks and Hooten, 2013, one widely used representation of smooth graph signals is to associate the smoothness property with a Gaussian graphical model where the precision matrix has the form of a graph Laplacian (Dong et al., 2016; Egilmez et al., 2016). The probabilistic model defined on the graph signal then naturally gives rise to a likelihood for the observed samples, and thus much of the literature in this area focuses on developing specialized algorithms to efficiently solve optimization problems that allow reconstruction of the underlying latent graph. For more information about graph learning and signal processing in general see the survey papers of Dong et al., 2019 and Mateos et al., 2019.

To position the present work in comparison to the ‘graph learning’ literature, our contributions are twofold. First, in population genetics, it is impossible to collect individual genotypes across all the geographic locations and, as a result, we often work with many, often the majority, of nodes having missing data. As far as we are aware, none of the work in graph signal processing considers this scenario and thus their algorithms are not directly applicable to our setting. In addition, if the number of the observed nodes is much smaller than the number of nodes of a graph, one can project the large matrices associated with the graph to the space of observed nodes, therefore allowing for fast and efficient computation. Second, highly missing nodes in the observed signals can result in significant degradation of the quality of the reconstructed graph unless it is regularized properly. Motivated by the Voronoi-tessellation-based prior adopted in EEMS (Petkova et al., 2016), we propose regularization that encourages spatial smoothness in the edge weights.

Building on advances in graph learning, we introduce a method, Fast Estimation of Effective Migration Surfaces (FEEMS), that uses optimization to obtain penalized-likelihood-based estimates of effective migration parameters. In contrast to EEMS which uses a node-specific parameterization of effective migration, we optimize over edge-specific parameters allowing for more flexible migration processes to be fit, such as spatial anisotropy, in which the migration process is not invariant to rotation of the coordinate system (e.g. migration is more extensive along a particular axis). Although we developed this model as a Gaussian Markov Random Field, the resulting likelihood has key similarities to the EEMS model, in that it is a Wishart-distribution that is a function of a genetic distance matrix. Expected genetic distances in both models can be interpreted as ‘resistance distances’ (McRae, 2006).

To fit the model, rather than using MCMC, we develop a fast quasi-Newton optimization algorithm (Nocedal and Wright, 2006) and a cross-validation approach for choosing the penalty parameter used in the penalized likelihood. We demonstrate the method using coalescent simulations and an application to a dataset of gray wolves from North America. The output is comparable to the results of EEMS but is provided in orders of magnitude less time. With this improvement in speed, FEEMS opens up the ability to perform fast exploratory data analysis of spatial population structure.

Results

Overview of FEEMS

Figure 1 shows a visual schematic of the FEEMS method. The input data are genotypes and spatial locations (e.g. latitudes and longitudes) for a set of individuals sampled across a geographic region. We construct a dense spatial grid embedded in geographic space where nodes represent sub-populations, and we assign individuals to nodes based on spatial proximity (see Appendix 1—figure 1 for a visualization of the grid construction and node assignment procedure). The density of the grid is user defined and must be explored to appropriately balance model mis-specification and computational burden. As the density of the lattice increases, the model is similar to discrete approximations used for continuous spatial processes, but the increased density comes at the cost of computational complexity.

Schematic of the FEEMS model: The full panel shows a schematic of going from the raw data (spatial coordinates and genotypes) through optimization of the edge weights, representing effective migration, to convergence of FEEMS to a local optima.

(A) Map of sample coordinates (black points) from a dataset of gray wolves from North America (Schweizer et al., 2016). The input to FEEMS are latitude and longitude coordinates as well as genotype data for each sample. (B) The spatial graph edge weights after random initialization uniformly over the graph to begin the optimization algorithm. (C) The edge weights after 10 iterations of running FEEMS, when the algorithm has not converged yet. (D) The final output of FEEMS after the algorithm has fully converged. The output is annotated with important features of the visualization.

Details on the FEEMS model are described in the Materials and methods section, however at a high level, we assume exchangeability of individuals within each sub-population and estimate allele frequencies, f^j(k), for each sub-population, indexed by k, and single nucleotide polymorphism (SNP), indexed by j, under a simple Binomial sampling model. We also use the recorded sample sizes at each node to model the precision of the estimated allele frequency. With the estimated allele frequencies in hand, we model the data at each SNP using an approximate Gaussian model whose covariance is, up to constant factors, shared across all SNPs—in other words, after rescaling by SNP-specific variation factors, we assume that the set of observed frequencies at each SNP is an independent realization of the same spatial process. The latent frequency variables, fj(k), are modeled as a Gaussian Markov Random Field (GMRF) with a sparse precision matrix determined by the graph Laplacian and a set of residual variances that vary across SNPs. The pseudo-inverse of the graph Laplacian in a GMRF is inherently connected to the notion of resistance distance in an electrical circuit (Hanks and Hooten, 2013) that is often used in population genetics to model the genetic differentiation between sub-populations (McRae, 2006). The graph’s weighted edges, denoted by wij between nodes i and j, represent gene-flow between the sub-populations (Friedman et al., 2008; Hanks and Hooten, 2013; Petkova et al., 2016). The Gaussian approximation has the advantage that we can analytically marginalize out the latent frequency variables. The resulting likelihood of the observed frequencies shares a number of similarities to that of EEMS (see Materials and methods).

To prevent over-fitting we use penalized maximum likelihood to estimate the edge weights of the graph. Our overall goal is thus to solve the following optimization problem:

w^=argminlwu (w)+ϕλ(w),

where 𝒘 is a vector that stores all the unique elements of the weighted adjacency matrix, 𝒍 and 𝒖 are element-wise non-negative lower and upper bounds for 𝒘, (𝒘) is the negative log-likelihood function that comes from the GMRF model described above, and ϕλ(𝒘) is a penalty that controls how constant or smooth the output migration surface will be and is controlled by the hyperparameter λ>0. Writing 𝒱 to denote the set of nodes in the graph and (i)𝒱 to denote the subset of nodes that have edges connected to node i, our penalty is given by

ϕλ(𝒘)=λ2i𝒱k,(i)(log(ewik/w^0-1)-log(ewi/w^0-1))2.

This function serves to penalize large differences between the weights wik and wi on edges that are adjacent, that is, penalizing differences for any pair of edges that share a common node. The tuning parameter λ controls the overall strength of the penalization placed on the output of the migration surface—if λ is large, the fitted surface will favor a homogeneous set of inferred migration weights on the graph, while if λ is low, more flexible graphs can be fitted to recover richer local structure, but this suffers from the potential for over-fitting. The tuning parameter λ is selected by evaluating the model’s performance at predicting allele frequencies at held out locations using leave-one-out cross-validation (see Materials and methods ‘Leave-one-out cross-validation to select tuning parameters’).

The scale parameter w^0 is chosen first fitting a ‘constant w’ model, which is a spatially homogeneous isolation-by-distance model constrained to have a single w value for all edges. In the ϕλ penalty, for adjacent edges (i,k) and (i,), if wik and wi are large (relative to w^0) then the corresponding term of the penalty is approximately proportional to (wik-wi)2, penalizing differences among neighboring edges on a linear scale; if instead wik and wi are small relative to w^0, then the penalty is approximately proportional to (log(wik)-log(wi))2, penalizing differences on a logarithmic scale. In fact, it is also possible to consider treating this scale parameter as a second tuning parameter—we can define a penalty function ϕλ,α(𝒘)=λ2i𝒱k,(i)(log(eαwik-1)-log(eαwi-1))2, and explore the solution across different values of both λ and α. However, we find that empirically choosing α=1/w^0 offers good performance as well as an intuitive interpretation (i.e. scaling edge weights wik with reference to the constant-w model), and allows us to avoid the computational burden of searching a two-dimensional tuning parameter space.

We use sparse linear algebra routines to efficiently compute the objective function and gradient of our parameters, allowing for the use of widely applied quasi-Newton optimization algorithms (Nocedal and Wright, 2006) implemented in standard numerical computing libraries like scipy (Virtanen et al., 2020) (RRID:SCR_008058). See the Materials and methods section for a detailed description of the statistical models and algorithms used.

Evaluating FEEMS on ‘out of model’ coalescent simulations

While our statistical model is not directly based on a population genetic process, it is useful to see how it performs on simulated data under the coalescent stepping stone model (Figure 2, also see Appendix 1—figure 2 for additional scenarios). In these simulations we know, by construction, the model we fit (FEEMS) is different from the true model we simulate data under (the coalescent), allowing us to assess the robustness of the fit to a controlled form of model mis-specification.

FEEMS fit to coalescent simulations: We run FEEMS on coalescent simulations, varying the migration history (columns) and sampling design (rows).

In each simulation, we used leave-one-out cross-validation (at the level of sampled nodes) to select the smoothness parameter λ. The first column (A–C) shows the ground-truth and fit of FEEMS to coalescent simulations with a homogeneous migration history that is, a single migration parameter for all edge weights. Note that the ground-truth simulation figures (A,D,G) display coalescent migration rates, not fitted effective migration rates output by FEEMS. The second column (D–F) shows the ground truth and fit of FEEMS to simulations with a heterogeneous migration history that is, reduced gene-flow, with 10-fold lower migration, in the center of the habitat. The third column (G–I) shows the ground truth and fit of FEEMS to an anisotropic migration history with edge weights facing east-west having five fold higher migration than north-south. The second row (B,E,H) shows a sampling design with no missing observations on the graph. The final row (C,F,I) shows a sampling design with 80% of nodes missing at random.

The first migration scenario (Figure 2A–C) is a spatially homogeneous model where all the migration rates are set to be a constant value on the graph, this is equivalent to simulating data under an homogeneous isolation-by-distance model. In the second migration scenario (Figure 2D–E), we simulate a non-homogeneous process by representing a geographic barrier to migration, lowering the migration rates by a factor of 10 in the center of the habitat relative to the left and right regions of the graph. Finally, in the third migration scenario (Figure 2G–I), we simulate a pattern which corresponds to anisotropic migration with edges that point east/west being assigned to a fivefold higher migration rate than edges pointing north/south. For each migration scenario, we simulate two sampling designs. In the first ‘dense-sampling’ design (Figure 2B,E,I) we sample individuals for every node of the graph. Next, in the ‘sparse-sampling’ design (Figure 2C,F,J) we sample individuals for only a randomly selected 20% of the nodes.

For each coalescent simulation, we used leave-one-out cross-validation (at the level of sampled nodes) to select the smoothness parameter λ. In the homogeneous migration simulations, the best value for the smoothness parameter, as determined by the grid value with the lowest leave-one-out cross-validation error, is λcv=100 in both sampling scenarios with complete and missing data. In the heterogeneous migration simulations λcv=0.298 with no missing data and λcv=37.927 with missing data. Finally, in the anisotropic simulations with no missing data λcv=0.298 and with missing data λcv=0.042. We note the magnitude of the selected λ depends on the scale of the loss function so comparisons across different datasets are not generally interpretable.

With regard to the visualizations of effective migration, FEEMS performs best when all the nodes are sampled on the graph, that is, when there is no missing data (Figure 2B,E,H). Interestingly, in the simulated scenarios with many missing nodes, FEEMS can still partly recover the migration history, including the presence of anisotropic migration (Figure 2I). A sampling scheme with a central gap leads to a slightly narrower barrier in the heterogeneous migration scenario (Appendix 1—figure 2I) and for the anisotropic scenario, a degree of over-smoothness in the northern and southern regions of the center of the graph (Appendix 1—figure 2N). For the missing at random sampling design, FEEMS is able to recover the relative edge weights surprisingly well for all scenarios, with the inference being the most challenging when there is anisotropic migration. The potential for FEEMS to recover anisotropic migration is novel relative to EEMS, which was parameterized for fitting non-stationary isotropic migration histories and produces banding patterns perpendicular to the axis of migration when applied to data from anisotropic coalescent simulations (Petkova et al., 2016, Supplementary Figure 2; see also Appendix 1 ‘Edge versus node parameterization’ for a related discussion). Overall, even with sparsely sampled graphs, FEEMS is able to produce visualizations that qualitatively capture the migration history in coalescent simulations.

Application of FEEMS to genotype data from North American gray wolves

To assess the performance of FEEMS on real data, we used a previously published dataset of 111 gray wolves sampled across North America typed at 17,729 SNPs (Schweizer et al., 2016; Appendix 1—figure 5). This dataset has a number of advantageous features that make it a useful test case for evaluating FEEMS: (1) The broad sampling range across North America includes a number of relevant geographic features that, a priori, could conceivably lead to restricted gene-flow averaged throughout the population history. These geographic features include mountain ranges, lakes, and islands. (2) The scale of the data is consistent with many studies for non-model systems whose spatial population structure is of interest. For instance, the relatively sparse sampling leads to a challenging statistical problem where there is the potential for many unobserved nodes (sub-populations), depending the density of the grid chosen.

Before applying FEEMS, we confirmed a signature of spatial structure in the data through regressing genetic distances on geographic distances and top genetic PCs against geographic coordinates (Appendix 1—figures 6, 7, 8, 9). We also ran multiple replicates of ADMIXTURE for K=2 to K=8, selecting for each K the highest likelihood run among replicates to visualize (Appendix 1—figure 10). As expected in a spatial genetic dataset, nearby samples have similar admixture proportions and continuous gradients of changing ancestries are spread throughout the map (Bradburd et al., 2018). Whether such gradients in admixture coefficients are due to isolation by distance or specific geographic features that enhance or diminish the levels of genetic differentiation is an interpretive challenge. Explicitly modeling the spatial locations and genetic distance jointly using a method like EEMS or FEEMS is exactly designed to explore these types of questions in the data (Petkova, 2013; Petkova et al., 2016).

We first show FEEMS results for four different values of the smoothness parameter, λ from large λ=100 to small λ=0.0008 (Figure 3). One interpretation of our regularization penalty is that it encourages fitting models of homogeneous and isotropic migration. When λ is very large (Figure 3A), we see FEEMS fits a model where all of the edge weights on the graph nearly equal the mean value, hence all the edge weights are colored white in the relative log-scale. In this case, FEEMS is fitting a relatively homogeneous migration model where all the estimated edge weights get assigned nearly the same value on the graph. As we sequentially lower the penalty parameter, (Figure 3B,C,D) the fitted graph begins to appear more complex and heterogeneous as expected (discussed further below). Figure 3E shows the cross-validation error for a pre-defined grid of λ values (also see Appendix 1—figure 6 for visualizations of the fitted versus genetic distance on the full dataset).

The fit of FEEMS to the North American gray wolf dataset for different choices of the smoothing regularization parameter λ: (A) λ=100, (B) λ=2.06, (C) λ=0.04, and (D) λ=0.0008.

As expected, when λ decreases from large to small (A–D), the fitted graph becomes less smooth and eventually over-fits to the data, revealing a patchy surface in (D), whereas earlier in the regularization path FEEMS fits a homogeneous surface with all edge weights having nearly the same fitted value, as in (A). (E) shows the mean square error between predicted and held-out allele frequencies output by running leave-one-out cross-validation to select the smoothness parameter λ. The cross-validation error is minimized over a pre-selected grid at an intermediate value of λ=2.06 as shown in (B).

The cross-validation approach finds the optimal value of λ to be 2.06. This solution visually appears to have a moderate level of regularization and aligns with several known landscape features (Figure 4). Spatial features in the FEEMS visualization qualitatively matches the structure plot output from ADMIXTURE using K=6 (Appendix 1—figure 10). We add labels on the figure to highlight a number of pertinent features: (A) St. Lawrence Island, (B) the coastal islands and mountain ranges in British Columbia, (C) the boundary of Boreal Forest and Tundra eco-regions in the Shield Taiga, (D) Queen Elizabeth Islands, (E) Hudson Bay, and (F) Baffin Island. Many of these features were described in Schweizer et al., 2016 by interpretation of ADMIXTURE, PCA, and FST statistics. FEEMS is able to succinctly provide an interpretable view of these data in a single visualization. Indeed many of these geographic features plausibly impact gray wolf dispersal and population history (Schweizer et al., 2016).

FEEMS applied to a population genetic dataset of North American gray wolves: We show the fit of FEEMS applied to a previously published dataset of North American gray wolves.

Leave-one-out cross-validation (at the level of sampled nodes) was used to select the smoothness parameter λ=2.06. We show the fitted parameters in log-scale with lower effective migration shown in orange and higher effective migration shown in blue. The bold text letters highlights a number of known geographic features that could have plausibly influenced wolf migration over time: (A) St. Lawrence Island, (B) Coastal mountain ranges in British Columbia, (C) The boundary of Boreal Forest and Tundra eco-regions in the Shield Taiga, (D) Queen Elizabeth Islands, (E) Hudson Bay, and (F) Baffin Island. We also display two insets to help interpret the results and compare them to EEMS. In the top left inset we show a map of sample coordinates colored by an ecotype label provided by Schweizer et al., 2016. These labels were devised using a combination of genetic and ecological information for 94 ‘un-admixed’ gray wolf samples, and the remaining samples were labeled ‘Other’. We can see these ecotype labels align well with the visualization output provided by FEEMS. In the right inset, we display a visualization of the posterior mean effective migration rates from EEMS.

Comparison to EEMS

We also ran EEMS on the same gray wolf dataset. We used default parameters provided by EEMS but set the number of burn-in iterations to 20×106, MCMC iterations to 50×106, and thinning intervals to 2000. We were unable to run EEMS in a reasonable run time (3 days) for the dense spatial grid of 1207 nodes so we ran EEMS and FEEMS on a sparser graph with 307 nodes.

We find that FEEMS is multiple orders of magnitude faster than EEMS, even when running multiple runs of FEEMS for different regularization settings on both the sparse and dense graphs (Table 1). We note that constructing the graph and fitting the model with very low regularization parameters are the most computationally demanding steps in running FEEMS.

Table 1
Runtimes for FEEMS and EEMS on the North American gray wolf dataset.

We show a table of runtimes for FEEMS and EEMS for two different grid densities, a sparse grid with 307 nodes and a dense grid with 1207 nodes. The second row shows the FEEMS run-times for applying leave-one-out cross-validation to select λ. The third row shows the run-times when applying FEEMS at the best λ value selected using cross-validation. FEEMS is orders of magnitude faster than EEMS, even when using cross-validation to select λ. Runtimes are based on computation using Intel Xeon E5-2680v4 2.4 GHz CPUs with 5 Gb RAM reserved using the University of Chicago Midway2 cluster.

MethodSparse grid (run-time)Dense grid (run-time)
EEMS27.43 hrN/A
FEEMS (Cross-validation)10 min 32 s1.03 hr
FEEMS (Best λ)1.23 s4.08 s

We find that many of the same geographic features that have reduced or enhanced gene-flow are concordant between the two methods. The EEMS visualization, qualitatively, best matches solutions of FEEMS with lower λ values (Figure 4, Appendix 1—figure 11); however, based on the ADMIXTURE results, visual inspection in relation to known geographical features and inspection of the observed vs fitted dissimilarity values (Appendix 1—figures 14, 22), we find these solutions to be less satisfying compared to the FEEMS solution found with λ chosen by leave-one-out cross-validation. We note that in many of the EEMS runs the MCMC appears to not have converged (based on visual inspection of trace plots) even after a large number of iterations.

Discussion

FEEMS is a fast approach that provides an interpretable view of spatial population structure in real datasets and simulations. We want to emphasize that beyond being a fast optimization approach for inferring population structure, our parameterization of the likelihood opens up a number of exciting new directions for improving spatial population genetic inference. Notably, one major difference between EEMS and FEEMS is that in FEEMS each edge weight is assigned its own parameter to be estimated, whereas in EEMS each node is assigned a parameter and each edge is constrained to be the average effective migration between the nodes it connects (see Materials and methods and Appendix 1 ‘Edge versus node parameterization’ for details). The node-based parameterization in EEMS makes it difficult to incorporate anisotropy and asymmeteric migration (Lundgren and Ralph, 2019). As we have shown here, FEEMS’s simple and novel parameterization already has potential to fit anisotropic migration (as shown in coalescent simulations) and may be extendable to other more complex migration processes (such as long-range migration, see below).

One general challenge, which is not unique to this method, is selecting the tuning parameters controlling the strength of regularization (λ in our case). A natural approach is to use cross-validation, which estimates the out-of-sample fit of FEEMS for a particular choice of λ. We used leave-one-out cross-validation, leaving one sampled population out at a time, and find such an approach works well based on the coalescent simulations and application to the North American wolf data. That said, we sometimes found high variability in the selected solution when we used cross-validation with fewer folds (e.g. five-fold versus leave-one-out, results not shown). We expect this happens when the number of sampled populations is small relative to the complexity of the gene flow landscape, and we recommend using leave-one-out cross-validation in general. We also find it useful to fit FEEMS to a sequential grid of regularization parameters and to look at what features are consistent or vary across multiple fits. Informally, one can gain an indication of the strongest features in the data by looking at the order they appear in the regularization path that is, what features overcome the strong penalization of smoothness in the data and that are highly supported by the likelihood. For example, early in the regularization path, we see regions of reduced gene-flow occurring in the west coast of Canada that presumably correspond to Coastal mountain ranges and islands in British Columbia (Figure 3B) and this reduced gene-flow appears throughout more flexible fits with lower λ.

An important caveat is that the objective function we optimize is non-convex so any visualization output by FEEMS should be considered a local optimum and, keeping in mind that with different initialization one could get different results. That said, for the datasets investigated, we found the output visualizations were not sensitive to initialization, and thus our default setting is constant initialization fitted under an homogeneous isolation by distance model (See Materials and methods).

When comparing to EEMS, we found FEEMS to be much faster (Table 1). While this is encouraging, care must be taken because the goals and outputs of FEEMS and EEMS have a number of differences. FEEMS fits a sequential grid of solutions for different regularization parameters, whereas EEMS infers a posterior distribution and outputs the posterior mean as a point estimate. FEEMS is not a Bayesian method and unlike EEMS, which explores the entire landscape of the posterior distribution, FEEMS returns a particular point estimate: a local minimum point of the optimization landscape. Setting the prior hyper-parameters in EEMS act somewhat like a choice of the tuning parameter λ, except that EEMS uses hierarchical priors that in principle allow for exploration of multiple scales of spatial structure in a single run, but requires potentially long computation times for adequate MCMC convergence.

Like EEMS, FEEMS is based on an assumed underlying spatial graph of populations exchanging gene flow with neighboring populations. While the inferred migration rates explain the data under an assumed model, it is important for users and readers of FEEMS results to keep in mind the range and density of the chosen grid when interpreting results. We note that using a denser grid has the two potential advantages of providing improved approximation for continuously distributed species, as well as a more flexible model space to fit the data.

Depending on the scale of the analysis and the life history of the species, the process of assuming and assigning a single geographic location for each individual is a potential limitation of the modeling framework used here. For instance, the North American wolves studied here are understood to be generally territorial with individual ranges that are on the scale of 103 km2 (Burch et al., 2005), which is small relative to the greater than 106 km2 scale of our analysis. Thus, modeling individual wolves with single locations may not generally be problematic. However, at the boundary of the Boreal forest and Tundra, there are wolves with larger annual ranges and seasonal migrations that track caribou herds roughly north-south over distances of 1000 km (Musiani et al., 2007), and the wolves in the study were sampled in the winter (Musiani et al., 2007; Schweizer et al., 2016). If the samples were instead obtained in the summer, the position of the inferred low migration feature near the boundary of the Boreal Forest (marked 'C' in Figure 4) would presumably shift northward. The general cautionary lesson is that one must be careful when interpreting these maps to consider the life history of dispersal for the organism under study during the interpretation of results. Extending the methodology to incorporate knowledge of uncertainty in position or known dispersal may be an interesting direction for future work.

One natural extension to FEEMS, pertinent to a number of biological systems, is incorporating long-range migration (Pickrell and Pritchard, 2012; Bradburd et al., 2016). In this work, we have used a triangular lattice embedded in geographic space and enforced smoothness in nearby edge weights through penalizing their squared differences (see Materials and methods). We could imagine changing the structure of the graph by adding edges to allow for long-range connections; however, our current regularization scheme would not be appropriate for this setting. Instead, we could imagine adding an additional penalty to the objective, which would only allow a few long range connections to be tolerated. This could be considered to be a combination of two existing approaches for graph-based inference, graphical lasso (GLASSO) and graph Laplacian smoothing, combining the smoothness assumption for nearby connections and the sparsity assumption for long-range connections (Friedman et al., 2008; Wang et al., 2016). Another potential methodological avenue to incorporate long-range migration is to use a ‘greedy’ approach. We could imagine adding long-range edges one a time, guided by re-fitting the spatial model and taking a data-driven approach to select particular long-range edges to include. The proposed greedy approach could be considered to be a spatial graph analog of TreeMix (Pickrell and Pritchard, 2012).

Another interesting extension would be to incorporate asymmetric migration into the framework of resistance distance and Gaussian Markov Random Field based models. FEEMS, like EEMS, used a likelihood that is based on resistance distances, which are limited in their ability to model asymmetric migration (Lundgren and Ralph, 2019). Recently, Hanks, 2015 developed a promising new framework for deriving the stationary distribution of a continuous time stochastic process with asymmetric migration on a spatial graph. Interestingly, the expected distance of this process has similarities to the resistance distance-based models, in that it depends on the pseudo-inverse of a function of the graph Laplacian. Hanks, 2015 used MCMC to estimate the effect of known covariates on the edge weights of the spatial graph. Future work could adapt this framework into the penalized optimization approach we have considered here, where adjacent edge weights are encouraged to be smooth.

Finally, when interpreted as mechanistic rather than statistical models, both EEMS and FEEMS implicitly assume time-stationarity, so the estimated migration parameters should be considered to be ‘effective’ in the sense of being averaged over time in a reality where migration rates are dynamic and changing (Pickrell and Reich, 2014). The MAPS method is one recent advance that utilizes long stretches of shared haplotypes between pairs of individuals to perform Bayesian inference of time varying migration rates and population sizes (Al-Asadi et al., 2019). With the growing ability to extract high quality DNA from ancient samples, another exciting future direction would be to apply FEEMS to ancient DNA datasets over different time transects in the same focal geographic region to elucidate changing migration histories (Mathieson et al., 2018). There are a number of technical challenges in ancient DNA data that make this a difficult problem, particularly high levels of missing and low-coverage data. Our modeling approach could be potentially more robust, in that it takes allele frequencies as input, which may be estimable from dozens of ancient samples at the same spatial location, in spite of high degrees of missingness (Korneliussen et al., 2014).

In closing, we look back to a review titled ‘How Can We Infer Geography and History from Gene Frequencies?’ published in 1982 (Felsenstein, 1982). In this review, Felsenstein laid out fundamental open problems in statistical inference in population genetic data, a few of which we restate as they are particularly motivating for our work:

  • For any given covariance matrix, is there a corresponding migration matrix which would be expected to lead to it? If so, how can we find it?

  • How can we characterize the set of possible migration matrices which are compatible with a given set of observed covariances?

  • How can we confine our attention to migration patterns which are consistent with the known geometric co-ordinates of the populations?

  • How can we make valid statistical estimates of parameters of stepping stone models?

The methods developed here aim to help address these longstanding problems in statistical population genetics and to provide a foundation for future work to elucidate the role of geography and dispersal in ecological and evolutionary processes.

Materials and methods

Model description

Request a detailed protocol

See Appendix 1 ‘Mathematical notation’ for a detailed description of the notation used to describe the model. To visualize and model spatial patterns in a given population genetic dataset, FEEMS uses an undirected graph, 𝒢=(𝒱,) with |𝒱|=d, where nodes represent sub-populations and edge weights (wk)(k,) represent the level of gene-flow between sub-populations k and . For computational convenience, we assume 𝒢 is a highly sparse graph, specifically a triangular grid that is embedded in geographic space around the sample coordinates. We observe a genotype matrix, 𝒀n×p, with n rows representing individuals and p columns representing SNPs. We imagine diploid individuals are sampled on the nodes of 𝒢 so that yij(k){0,1,2} records the count of some arbitrarily predefined allele in individual i, SNP j, on node k𝒱. We assume a commonly used simple Binomial sampling model for the genotypes:

(1) yij(k)| fj(k)Binomial(2,fj(k)),

where conditional on fj(k) for all j,k, the yij(k)’s are independent. We then estimate an allele frequency at each node and SNP by maximum likelihood:

f^j(k)=i=1nkyij(k)2nk,

where nk is the number of individuals sampled at node k. We estimate allele frequencies at o of the observed nodes out of d total nodes on the graph. From Equation (1), the estimated frequency in a particular sub-population, conditional on the latent allele frequency, will approximately follow a Gaussian distribution:

f^j(k)| fj(k)𝒩(fj(k),fj(k)(1fj(k))2nk).

Using vector notation, we represent the joint model of estimated allele frequencies as:

(2) f^j| fj𝒩o(Afj,diag(df,n)),

where 𝒇^j is a o×1 vector of estimated allele frequencies at observed nodes, 𝒇j is a d×1 vector of latent allele frequencies at all the nodes (both observed and unobserved), and 𝑨 is a o×d node assignment matrix where 𝑨k=1 if the kth estimated allele frequency comes from sub-population and 𝑨k=0 otherwise; and diag(𝒅𝒇,𝒏) denotes a o×o diagonal matrix whose diagonal elements corresponds to the appropriate variance term at observed nodes.

To summarize, we estimate allele frequencies from a subset of nodes on the graph and define latent allele frequencies for all the nodes of the graph. The assignment matrix 𝑨 maps these latent allele frequencies to our observations. Our summary statistics (the data) are thus (𝑭^,𝒏) where 𝑭^ is a o×p matrix of estimated allele frequencies and 𝒏 is a o×1 vector of sample sizes for every observed node. We assume the latent allele frequencies come from a Gaussian Markov Random Field:

(3) 𝒇j𝒩d(μj𝟏,μj(1-μj)𝑳),

where 𝑳 is the graph Laplacian, † represents the pseudo-inverse operator, and μj represents the average allele frequency across all of the sub-populations. Note that the multiplication by the SNP-specific factor μj(1-μj) ensures that the variance of the latent allele frequencies vanishes as the average allele frequency approaches to 0 or 1. One interpretation of this model is that the expected squared Euclidean distance between latent allele frequencies on the graph, after being re-scaled by μj(1-μj), is exactly the resistance distance of an electrical circuit (McRae, 2006; Hanks and Hooten, 2013):

E[(fj(k)fj())2]μj(1μj)=rk,where rk=(oko)L(oko)=Lkk2Lk+L,

where 𝒐i is a one-hot vector (i.e. storing a 1 in element i and zeros elsewhere). It is known that the resistance distance rk is equivalent to the expected commute time between nodes k and of a random walker on the weighted graph 𝒢 (Chandra et al., 1996). Additionally, the model (Equation 3) forms a Markov random field, and thus any latent allele frequency fj(k) is conditionally independent of all other allele frequencies given its neighbors which are encoded by nonzero elements of 𝑳 (Lauritzen, 1996; Koller and Friedman, 2009). Since we use a triangular grid embedded in geographic space to define the graph 𝒢, the pattern of nonzero elements is prefixed by the structure of the sparse traingular grid.

Using the law of total variance formula, we can derive from (Equations 2, 3) an analytic form for the marginal likelihood. Before proceeding, however, we further approximate the model by assuming 12fj(k)(1-fj(k))σ2μj(1-μj) for all j and k (see Appendix 1 ‘Estimating the edge weights under the exact likelihood model’ for the data model without this approximation). This assumption is mainly for computational purposes and may be a coarse approximation in general. On the other hand, the assumption is not too strong if we exclude SNPs with extremely rare allele frequencies, and more importantly, we find it leads to a good empirical performance, both statistically and computationally. With this approximation, the residual variance parameter σ2 is still unknown and needs to be estimated.

Under (Equation 2, 3), the law of total variance formula leads to specific formulas for the mean and variance structure as given in (Equation 4). With those results, we arrive at the following approximate marginal likelihood:

(4) 𝒇^j𝒩o(μj𝟏,μj(1-μj)(𝑨𝑳𝑨+σ2diag(𝒏-1))),

where diag(𝒏-1) is a o×o diagonal matrix computed from the sample sizes at observed nodes. We note the marginal distribution of f^j is not necessarily a Gaussian distribution; however, we use a Gaussian approximation to facilitate computation.

To remove the SNP means we transform the estimated frequencies by a contrast matrix, 𝑪(o-1)×o, that is orthogonal to the one-vector:

(5) 𝑪𝒇^jμj(1-μj)𝒩o-1(,𝑪𝑨𝑳𝑨𝑪+σ2𝑪diag(𝒏-1)𝑪).

Let 𝚺^=1p𝑭^s𝑭^s be the o×o sample covariance matrix of estimated allele frequencies after re-scaling, that is, 𝑭^s is a matrix formed by rescaling the columns of 𝑭^ by μ^j(1-μ^j), where μ^j is an estimate of the average allele frequency (see above). We can then express the model in terms of the transformed sample covariance matrix:

(6) p𝑪𝚺^𝑪𝒲o-1(𝑪𝑨𝑳𝑨𝑪+σ2𝑪diag(𝒏-1)𝑪,p),

where 𝒲p denotes a Wishart distribution with p degrees of freedom. Note we can equivalently use the sample squared Euclidean distance (often refereed to as a genetic distance) as a summary statistic: letting 𝑫^ be the genetic distance matrix with 𝑫^k=j=1p(f^j(k)-f^j())2/pμ^j(1-μ^j), we have

𝑫^=𝟏diag(𝚺^)+diag(𝚺^)𝟏-2𝚺^,

and so

𝑪𝑫^𝑪=-2𝑪𝚺^𝑪,

using the fact that the contrast matrix 𝑪 is orthogonal to the one-vector. Thus, we can use the same spatial covariance model implied by the allele frequencies once we project the distances on to the space of contrasts:

-p2𝑪𝑫^𝑪𝒲o-1(𝑪𝑨𝑳𝑨𝑪+σ2𝑪diag(𝒏-1)𝑪,p).

Overall, the negative log-likelihood function implied by our spatial model is the following (ignoring constant terms):

(7) (w,σ2;CΣ^C)=ptr((CALAC+σ2Cdiag(n1)C)1CΣ^C)plogdet(CALAC+σ2Cdiag(n1)C)1,

where 𝒘m is a vectorized form of the non-zero lower-triangular entries of the weighted adjacency matrix 𝑾 (recall that the graph Laplacian is completely defined by the edge weights, 𝑳=diag(𝑾𝟏)-𝑾, so there is an implicit dependency here). Since the graph is a triangular lattice, we only need to consider the non-zero entries to save computational time, that is, not all sub-populations are connected to each other.

We note our model (Equation 6) assumes that the p SNPs are independent. This assumption is unlikely to hold when datasets are analyzed with SNPs that statistically covary (linkage disequilibrium). However, we note that the degree of freedom parameter does not affect the point estimate produced by FEEMS because it is treated as a constant term in the log-likelihood function.

One key difference between EEMS (Petkova et al., 2016) and FEEMS is how the edge weights are parameterized. In EEMS, each node is given an effective migration parameter mk for node k𝒱 and the edge weight is parameterized as the average between the nodes it connects, that is, wk=(mk+m)/2 for (k,). FEEMS, on the other hand, assigns a parameter to every nonzero edge-weight. The former has fewer parameters, with the specific consequence that it only allows isotropy and imposes an additional degree of similarity among edge weights; instead, in the latter, the edge weights are free to vary apart from the regularization imposed by the penalty. See Appendix 1 ‘Edge versus node parameterization’ and Appendix 1—figures 15, 17 for more details.

Penalty description

Request a detailed protocol

As mentioned previously, we would like to encourage that nearby edge weights on the graph have similar values to each other. This can be performed by penalizing differences between all edges connected to the same node, that is, spatially adjacent edges:

ϕλ,α(w)=λ2i𝒱k,(i)(log(eαwik1)log(eαwi1))2,

where, as before, (i) denotes the set of edges that is connected to node i. (As mentioned earlier, in practice we choose α=1/w^0, where w^0 is the solution for the ‘constant-w’ model, but we use the free parameter α here for full generality.) The function xlog(ex-1) (on positive values x(0,)) is approximately equal to x, for x much larger than 1, and is approximately equal to log(x), for x much smaller than 1. This means that our penalty function effectively penalizes differences on the log scale for edges (i,k) and (i,) with very small weights, but penalizes differences on the original non-log scale for edges with large weights. Using a logarithmic-scale penalty for edges with low weights (rather than simply penalizing (wik-wi)2) leads to smooth graphs for small edge values, and thus allow for an additional degree of flexibility across orders of magnitude of edge weights. The penalty parameter, λ, controls the overall contribution of the penalty to the objective function. It is convenient to write the penalty in matrix-vector form which we will use throughout:

(8) ϕλ,α(w)=λ2Δlog(eαw1)22,

where 𝚫 is a signed graph incidence matrix derived from a unweighted graph denoting if pairs of edges are connected to the same node. Specifically, in this expression, we treat 𝒘 as a vector of length || (i.e. the number of edges), and apply the function wlog(eαw-1) entrywise to this vector. For each pair adjacent edges (i,k) and (i,) in the graph, there is a corresponding row of 𝚫 with the value +1 in the entry corresponding to edge (i,k), a −1 in the entry corresponding to edge (i,), and 0’s elsewhere.

One might wonder whether it is possible to use the 1 norm in the penalty form Equation (8) in place of the 2 norm. While it is known that the 1 norm might increase local adaptivity and better capture the sharp changes of the underlying structure of the latent allele frequencies (e.g. Wang et al., 2016), in our case, we found an inferior performance when using the 1 norm over the 2 norm—in particular, our primary application of interest is the regime of highly missing nodes, that is, od, in which case the global smoothing seems somewhat necessary to encourage stable recovery of the edge weights at regions with sparsely observed nodes (see Appendix 1 ‘Smooth penalty with 1norm’). In addition, adding the penalty ϕλ,α(𝒘) allows us to implement faster algorithms to solve the optimization problem due to the differentiability of the 2 norm, and as a result, it leads to better overall computational savings and a simpler implementation.

Optimization

Request a detailed protocol

Putting Equation (7) and Equation (8) together, we infer the migration edge weights 𝒘^ by minimizing the following penalized negative log-likelihood function:

(9) w^=argminlwu (w,σ2;CΣ^C)+ϕλ,α(w)=argminlwu [ptr((CALAC+σ2Cdiag(n1)C)1CΣ^C)plogdet(CALAC+σ2Cdiag(n1)C)1+λ2Δlog(eαw1)22],

where 𝒍,𝒖+m represent respectively the entrywise lower- and upper bounds on 𝒘, that is, we constrain the lower- and upper bound of the edge weights to 𝒍 and 𝒖 throughout the optimization. When no prior information is available on the range of the edge weights, we often set 𝒍= and 𝒖=+.

One advantage of the formulation of Equation 9 is the use of the vector form parameterization 𝒘+m of the symmetric weighted adjacency matrix 𝑾+d×d. In our triangular graph 𝒢=(𝒱,), the number of non-zero lower-triangular entries is m=𝒪(d)d2, so working directly on the space of vector parameterization saves computational cost. In addition, this avoids the symmetry constraint imposed on the adjacency matrix 𝑾, hence making optimization easier (Kalofolias, 2016).

We solve the optimization problem using a constrained quasi-Newton optimization algorithm, specifically L-BFGS implemented in scipy (Byrd et al., 1995; Virtanen et al., 2020) (RRID:SCR_008058). Since our objective Equation 9 is non-convex, the L-BFGS algorithm is guaranteed to converge only to a local minimum. Even so, we empirically observe that local minima starting from different initial points are qualitatively similar to each other across many datasets. The L-BFGS algorithm requires gradient and objective values as inputs. Note the naive computation of the objective Equation 9 is computationally prohibitive because inverting the graph Laplacian has complexity 𝒪(d3). We take advantage of the sparsity of the graph and specific structure of the problem to efficiently compute gradient and objective values. In theory, our implementation has computational complexity of 𝒪(do+o3) per iteration which, in the setting of od, is substantially smaller than 𝒪(d3). It is possible to achieve 𝒪(do+o3) per-iteration complexity by using a solver that is specially designed for a sparse Laplacian system. In our work, we use sparse Cholesky factorization which may slightly slow down the per-iteration complexity (See Appendix Material for the details of the gradient and objective computation).

Estimating the residual variance and edge weights under the null model

Request a detailed protocol

For estimating the residual variance parameter σ2, we first estimate it via maximum likelihood assuming homogeneous isolation by distance. This corresponds to the scenario where every edge-weight in the graph is given the exact same unknown parameter value w0. Under this model, we only have two unknown parameters w0 and the residual variance σ2. We estimate these two parameters by jointly optimizing the marginal likelihood using a Nelder-Mead algorithm implemented in scipy (Virtanen et al., 2020) (RRID:SCR_008058). This requires only likelihood computations which are efficient due to the sparse nature of the graph. This optimization routine outputs an estimate of the residual variance σ^2 and the null edge weight w^0, which can be used to construct 𝑾(w^0) and in turn 𝑳(w^0).

One strategy we found effective is to fit the model of homogeneous isolation by distance and then fix the estimated residual variance σ^2 throughout later fits of the more flexible penalized models—See Appendix 1 ‘Jointly estimating the residual variance and edge weights’. Additionally, we find that initializing the edge weights to w^0 to be a useful and intuitive strategy to set the initial values for the entries of 𝒘 to the correct scale.

Leave-one-out cross-validation to select tuning parameters

Request a detailed protocol

FEEMS estimates one set of graph edge weights for each setting of the tuning parameters λ and α which control the smoothness of the fitted edge weights. Figure 3 shows that the estimated migration surfaces vary substantially depending on the particular choices of the tuning parameters, and indeed, due to the large fraction of unobserved nodes, it can highly over-fit the observed data unless regularized accordingly. To address the issue of selecting the tuning parameters, we propose using leave-one-out cross-validation to assess each fitted model’s generalization ability at held out locations.

To simplify the notation, we write the model Equation 4 for the estimated allele frequencies in SNP j as

𝒇^j𝒩o(𝝁j,𝚺j),

where

(10) 𝝁j=μj𝟏 and 𝚺j=μj(1-μj)(𝑨𝑳𝑨+σ2diag(𝒏-1)).

For each fold, we hold out one node from the set of observed nodes in the graph and use the rest of the nodes to fit FEEMS across a sequential grid of regularization parameters. Note that our objective function is non-convex, so the algorithm converges to different local minima for different regularization parameters, even with the same initial value w^0. To stabilize the cross-validation procedure, we recommend using a warm start strategy in which one solves the problem for the largest value of regularization parameters first and use this solution to initialize the algorithm at the next largest value of regularization parameters, and so on. Empirically, we find that using warm starts gives far more reliable model selection than with cold starts, where the problems over the sequence of parameters are solved independently with same initial value w^0. We suspect that the poor performance of leave-one-out cross-validation without warm starts is attributed to spatial dependency of allele frequencies and the large fraction of unobserved nodes. Without loss of generality, we assume that the last node has been held out. Re-writing the distribution of the observed frequencies according to the split of observed nodes,

f^j=(f^jtrf^jval)𝒩o((μjtrμjval),(ΣjtrΣjcovΣjcovΣjval)),

the conditional mean of the observed frequency f^jval on the held out node, given the rest, is given by

f^jval,pred=𝔼[f^jval|𝒇^jtr]=μjval+𝚺jcov𝚺jtr(𝒇^jtr-𝝁jtr)-1.

Using this formula, we can predict allele frequencies at held out locations using the fitted graph 𝑳^=𝑳^(λ,α) for each setting of tuning parameters λ and α. Note that in Equation (10), the parameters μj and σ are also unknown, and we use an estimate of the average allele frequency μ^j and the estimated residual variance σ^ from the ‘constant-w’ model (they are not dependent on λ and α). Then we select the tuning parameters λ and α that output the minimum prediction error averaged over all SNPs 1pjf^jval,predf^jval22, averaged over all the held out nodes (with o observed nodes in total). As mentioned earlier, in practice we choose α=1/w^0 and hence we can use the leave-one-out cross-validation to search for λ only, which allows us to avoid the computational cost of searching over the two-dimensional parameter space.

Comparison between FEEMS and EEMS models

Request a detailed protocol

At a high level, we can summarize the differences between FEEMS and EEMS as follows: (1) the likelihood functions of FEEMS and EEMS are slightly different as a function of the graph Laplacian 𝑳; (2) the migration rates are parameterized in terms of edge weights or in terms of node weights; and (3) EEMS is based on Bayesian inference and thus chooses a prior and studies the posterior distribution, while FEEMS is an optimization-based approach and thus chooses a penalty function and minimizes the penalized log-likelihood (in particular, the EEMS prior and the FEEMS penalty are both aiming for locally constant type migration surfaces). The last two points were already discussed in the above sections, so here we focus on the difference of the likelihoods between the two methods.

FEEMS develops the spatial model for the genetic differentiation through Gaussian Markov Random Field, but the resulting likelihood has similarities to EEMS (Petkova et al., 2016) which considers the pairwise coalescent times. Using our notation, we can write the EEMS model as

(11) -ν2𝑪𝑫^*𝑪𝒲o-1(σ*𝑪(12𝑨𝑳𝑨+diag(𝑨𝒒))𝑪,ν),

where ν[o-1,p] is the effective degree of freedom, σ*>0 is the scale nuisance parameter, and 𝒒 is a d×1 vector of the within-sub-population coalescent rates. 𝑫^* represents the genetic distance matrix without re-scaling, where the (k,)-th element is given by 𝑫^k=j=1p(f^j(k)-f^j())2/p. That is, unlike FEEMS, EEMS does not consider the SNP-specific re-scaling factor μj(1-μj) to account for the vanishing variance of the observed allele frequencies as the average allele frequency approaches to 0 or 1.

In Equation (11), the effective degree of freedom ν is introduced to account for the dependency across SNPs in close proximity. Because EEMS uses a hierarchical Bayesian model to infer the effective migration rates, ν is being estimated alongside other model parameters. On the other hand, FEEMS uses an optimization-based approach and the degrees of freedom has no influence on the point estimate of the migration rates. Besides the effective degree of freedom and the SNP-specific re-scaling by μj(1-μj), the EEMS and FEEMS likelihoods are equivalent up to constant factors, as long as only one individual is observed per node and the residual variance σ2 is allowed to vary across nodes—See Appendix 1 ‘Jointly estimating the residual variance and edge weights’ for details. The constant factors, such as σ*, can be effectively absorbed into the unknown model parameters 𝑳 and 𝒒 and therefore they do not affect the estimation of effective migration rates, up to constant factors.

Data description and quality control

Request a detailed protocol

We analyzed a population genetic dataset of North American gray wolves previously published in Schweizer et al., 2016. For this, we downloaded plink (RRID:SCR_001757) formatted files and spatial coordinates from https://doi.org/10.5061/dryad.c9b25. We removed all SNPs with minor allele frequency less than 5% and with missingness greater then 10%, resulting in a final set of 111 individuals and 17,729 SNPs.

Population structure analyses

Request a detailed protocol

We fit the Pritchard, Donnelly, and Stephens model (PSD) and ran principal components analysis on the genotype matrix of North American gray wolves (Price et al., 2006; Pritchard et al., 2000). For the PSD model, we used the ADMIXTURE software (RRID:SCR_001263) on the un-normalized genotypes, running five replicates per choice of K, from K=2 to K=8 (Alexander et al., 2009). For each K, we choose the one that achieved the highest likelihood to visualize. For PCA, we centered and scaled the genotype matrix and then ran sklearn (RRID:SCR_019053) implementation of PCA, truncated to compute 50 eigenvectors.

Grid construction

Request a detailed protocol

To create a dense triangular lattice around the sample locations, we first define an outer boundary polygon. As a default, we construct the lattice by creating a convex hull around the sample points and manually trimming the polygon to adhere to the geography of the study organism and balancing the sample point range with the extent of local geography using the following website https://www.keene.edu/campus/maps/tool/. We often do not exclude internal ‘holes’ in the habitat (e.g. water features for terrestrial animals), and let the model instead fit effective migration rates for those features to the extent they lead to elevated differentiation. We also emphasize the importance of defining the lattice for FEEMS as well as EEMS and suggest this should be carefully curated with prior biological knowledge about the system.

To ensure edges cover an equal area over the entire region, we downloaded and intersected a uniform grid defined on the spherical shape of earth (Sahr et al., 2003). These defined grids are pre-computed at a number of different resolutions, allowing a user to test FEEMS at different grid densities which is an important feature to explore.

Code availability

Request a detailed protocol

The code to reproduce the results of this paper and more can be found at https://github.com/jhmarcus/feems-analysis (Marcus and Ha, 2021a, copy archived at swh:1:rev:f2d7330f25f8a11124db09000918ae38ae00d4a7, Marcus and Ha, 2021b). A python (RRID:SCR_008394) package implementing the method can be found at https://github.com/Novembrelab/feems.

Appendix 1

Mathematical notation

We denote matrices using bold capital letters 𝑨. Bold lowercase letters are vectors 𝒂, and non-bold lowercase letters are scalars a. We denote by 𝑨-1 and 𝑨 the inverse and (Moore-Penrose) pseudo-inverse of 𝑨, respectively. We use 𝒚Np(𝝁,𝚺) to express that the random vector 𝒚 is modeled as a p-dimensional multivariate Gaussian distribution with fixed parameters 𝝁 and 𝚺 and use the conditional notation 𝒚|𝝁Np(𝝁,𝚺) if 𝝁 is random.

A graph is a pair 𝒢=(𝒱,), where 𝒱 denotes a set of nodes or vertices and 𝒱×𝒱 denotes a set of edges. Throughout we assume the graph 𝒢 is undirected, weighted, and contains no self loops, that is, (k,)(,k) and (k,k) and each edge (k,) is given a weight wk=wk>0. We write 𝑾 to indicate the symmetric weighted adjacency matrix, that is,

𝑾k={wk,if (k,),0,otherwise.

𝒘m is a vectorized form of the non-zero lower-triangular entries of 𝑾 where m=||/2 is the number of non-zero lower triangular elements. We denote by 𝑳=diag(𝑾𝟏)-𝑾 the graph Laplacian.

Gradient computation

In practice, we make a change of variable from 𝒘+m to 𝒛=log(𝒘)m and the algorithm is applied to the transformed objective function:

(exp(𝒛),σ2;𝑪𝚺^𝑪)+ϕλ,α(exp(𝒛))=~(𝒛,σ2;𝑪𝚺^𝑪)+ϕ~λ,α(𝒛).

After the change of variable, the objective value remains the same, whereas it follows from the chain rule that (~(𝒛)+ϕ~λ,α(𝒛))=((𝒘)+ϕλ,α(𝒘))𝒘 where indicates the Hadamard product or elementwise product—for notational convenience, we drop the dependency of on the quantities σ2 and 𝑪𝚺^𝑪. Furthermore, the computation of ϕλ,α(𝒘) is relatively straightforward, so in the rest of this section, we discuss only the computation of the gradient of the negative log-likelihood function with respect to 𝒘, that is, (𝒘).

Recall, by definition, the graph Laplacian 𝑳 implicitly depends on the variable 𝒘 through 𝑳=diag(𝑾𝟏)-𝑾. Throughout we assume the first o rows and columns of 𝑳 correspond to the observed nodes. With this assumption, our node assignment matrix has block structure 𝑨=[𝐈o×o| 0o×(d-o)]. To simplify some of the equations appearing later, we introduce the notation: we define

(12) 𝑳full:=𝑳+𝟏𝟏d,𝚺:=𝑨𝑳full-1𝑨+σ2diag(𝒏-1),

and

𝑴:=𝑪((𝑪𝚺𝑪)-1(𝑪𝚺^𝑪)(𝑪𝚺𝑪)-1-(𝑪𝚺𝑪)-1)𝑪.

Applying the chain rule and matrix derivatives, we can calculate:

(𝒘)=(𝒘)vec(𝑳)vec(𝑳)𝒘,

where vec is the vectorization operator and /vec(𝑳) and vec(𝑳)/𝒘 are 1×d2 vector and d2×d matrix, respectively, given by

(13) (𝒘)vec(𝑳)=pvec(𝑳full-1𝑨𝑴𝑨𝑳full-1,),vec(𝑳)𝒘=𝑺-𝑻.

Here, 𝑺 and 𝑻 are linear operators that satisfy 𝑺𝒘=diag(𝑾𝟏) and 𝑻𝒘=𝑾. Note 𝑺 and 𝑻 both have 𝒪(d) many nonzero entries, so we can perform sparse matrix multiplication to efficiently compute the matrix-vector multiplication /vec(𝑳)(𝑺-𝑻). On the other hand, the computation of /vec(𝑳) is more challenging as it requires inverting the full d×d matrix 𝑳full. Next, we develop a procedure that efficiently computes /vec(𝑳). We proceed by dividing the task into multiple steps.

1. Computing 𝚺-1

Recalling the block structure 𝑨=[𝐈o×o| 0o×(d-o)] of the node assignment matrix, we can write 𝚺 as:

Σ=(Lfull1)o×o+σ2diag(n1),

where (𝑳full-1)o×o denotes the o×o upper-left block of 𝑳full-1. Following Petkova et al., 2016, the inverse 𝚺-1 has the form

(14) 𝚺-1=𝑿+σ-2diag(𝒏),

for some matrix 𝑿o×o. Equating 𝚺𝚺-1=𝐈, it follows that

(15) [(Lfull1)o×o+σ2diag(n1)](X+σ2diag(n))=I[(Lfull1)o×o+σ2diag(n1)]X=σ2(Lfull1)o×odiag(n).

Therefore, 𝚺-1 can be obtained by solving the o×o linear system Equation (15) and plugging the solution into Equation (14). The challenge here is to compute (𝑳full-1)o×o without matrix inversion of the full-dimensional 𝑳full.

2. Computing (𝑳full-1)o×o

Let 𝑳full,o×o be the o×o block matrix corresponding to the observed nodes of 𝑳full, and similarly let 𝑳full,(d-o)×(d-o) and 𝑳full,o×(d-o)=𝑳full,(d-o)×o be the corresponding block matrices of 𝑳full, respectively. The inverse of (𝑳full-1)o×o is then given by the Schur complement of 𝑳full,(d-o)×(d-o) in 𝑳:

(16) [(𝑳full-1)o×o]-1=𝑳full,o×o-𝑳full,o×(d-o)(𝑳full,(d-o)×(d-o))-1𝑳full,(d-o)×o.

See also Hanks and Hooten, 2013, Petkova et al., 2016. Since every term in Equation (16) has sparse + rank-1 structure, the matrix multiplications can be performed fast. In addition, for the term (𝑳full,(d-o)×(d-o))-1, we can use the Sherman-Morrison formula so that the inverse is given explicitly by

(Lfull,(do)×(do))1=(L(do)×(do)+11d)1=L(do)×(do)11d+1L(do)×(do)11L(do)×(do)111L(do)×(do)1.

Hence, in order to compute (𝑳full,(d-o)×(d-o))-1𝑳full,(d-o)×o, we need to solve two systems of linear equations:

L(do)×(do)U=Lfull,(do)×o and L(do)×(do)u=1.

Note that the matrix 𝑳(d-o)×(d-o) is sparse, so both systems can be solved efficiently by performing sparse Cholesky factorization on 𝑳(d-o)×(d-o) (Hanks and Hooten, 2013). Alternatively, one can implement fast Laplacian solvers (Vishnoi, 2013) that solve the Laplacian system in time nearly linear in the dimension 𝒪(d). After we obtain [(𝑳full-1)o×o]-1 via sparse + rank-1 matrix multiplication and sparse Cholesky factorization, we can invert the o×o matrix to get (𝑳full-1)o×o.

3. Computing (𝑳full-1)d×o

We write

(𝑳full-1)d×o=[(𝑳full-1)o×o(𝑳full-1)(d-o)×o].

Using the inversion of the matrix in a block form, the (d-o)×o block component is given by

(17) (Lfull1)(do)×o=(Lfull,(do)×(do))1Lfull,(do)×o(A)(Lfull1)o×0(B)

Since each of the two terms (A) and (B) has been already computed in the previous step, there is no need to recompute them. In total, it requires a (do)×o matrix and o×o matrix multiplication.

4. Computing the full gradient

Going back to the expression of (𝒘) in Equation (13), and noting the block structure of the assignment matrix 𝑨, we have:

(w)vec(L)=pvec((Lfull1)d×oM(Lfull1)d×o).

Define Π1=1(1Σ11)11Σ1 which acts as a sort of projection to the space of constant vectors with respect to the inner product 𝒙,𝒚=𝒙𝚺-1𝒚. Using the identity 𝐈-𝚷𝟏=𝚺𝑪(𝑪𝚺𝑪)-1𝑪 (McCullagh, 2009), then we can write 𝑴 in terms of 𝚷𝟏:

(18) M=Σ1(IΠ1)Σ^Σ1(IΠ1)Σ1(IΠ1).

Since Π𝟏 is a rank-1 matrix, this expression of 𝑴 allows easier computation. Finally we can put together Equation (14), Equation (15), Equation (17), and Equation (18), to compute the gradient of the negative log-likelihood function with respect to the graph Laplacian.

Objective computation

The graph Laplacian 𝑳 is orthogonal to the one vector 1, so using the notation introduced in Equation (12), we can express our objective function as

(w)+ϕλ,α(w)=ptr((CΣC)1CΣ^C)plogdet(CΣC)1+λ2Δlog(eαw1)22.

With the identity 𝐈-𝚷𝟏=𝚺𝑪(𝑪𝚺𝑪)-1𝑪, the trace term is:

tr((CΣC)1CΣ^C)=tr(C(CΣC)1CΣ^)=tr(Σ1(IΠ1)Σ^).

The matrix inside the trace has been constructed in the gradient computation, see Equation (18). In terms of the determinant, we use the same approach considered in Petkova et al., 2016—in particular, concatenating 𝑪 and 𝟏, the matrix [𝑪| 1] is orthogonal, so it can be shown that

det(Σ)=det(11)det(CΣC)det(CC)det(1Σ11).

Rearranging terms and using the fact det(𝑼-1)=det(𝑼)-1 for any matrix 𝑼, we obtain:

det(CΣC)1=det(11)det(Σ1)det(CC)det(1Σ11)=o1Σ11det(Σ1).

We have computed 𝚺-1 in Equation (14), so each of the terms above can be computed without any additional matrix multiplications. Finally, the signed graph incidence matrix 𝚫 defined on the edges of the graph is, by construction, highly sparse with 𝒪(d) many nonzero entries. Hence we implement sparse matrix multiplication to evaluate the penalty function ϕλ,α(𝒘) while avoiding the full-dimensional matrix-vector product.

Estimating the edge weights under the exact likelihood model

When we developed the FEEMS model, we used the approximation 12fj(k)(1-fj(k))σ2μj(1-μj) for all SNPs j and all nodes k (see Equation 4) and estimated the residual variance σ2 under the homogeneous isolation by distance model. The primary reason of using this approximation was primarily computational. While the approximation is not too strong if SNPs with rare allele frequencies are excluded, it is also critical that the estimation quality of the migration rates is not affected. In this subsection we introduce the inferring procedure of the migration rates under the exact likellihood model and compare it with FEEMS.

Note that without approximation, we can calculate the exact analytical form for the marginal likelihood of the estimated frequency as follows (after removing the SNP means):

(19) Cf^jμj(1μj)𝒩o1(0,CALAC+Cdiag(n1)Adiag({1Lkk2}k=1d)AC),

where {ak}k=1d represents the vector 𝒂=(a1,,ad). Compared to the model Equation (5), this expression does not introduce the unknown residual variance parameter σ2 and instead each node has its own residual parameter given by (1-Lkk)/2. Because the residual parameters must be positive, this means that we have to search for the graphs that ensure Lkk1 for all nodes k. With that said, we can consider the following constrained optimization problem:

(20) w^=argminlwu{exact(w;CΣ^C)+ϕλ,α(w):Lkk1 for all k𝒱},

where 𝖾𝗑𝖺𝖼𝗍 is the negative log-likelihood function based on Equation (19) and ϕλ,α is the smooth penalty function defined earlier. The main difficulty of solving Equation (20) is that enforcing the constraint Lkk1 for all nodes k𝒱, requires full computation of the pseudo-inverse of a d×d matrix 𝑳 which is computationally demanding. We instead relax the constraint and consider the following form as a proxy for optimization Equation (20):

(21) w^=argminlwu {exact(w;CΣ^C)+ϕλ,α(w):Lkk1 for all observed nodes k}.

Note that the constraint Lkk1 is now placed at the observed nodes only, which can lead to computational savings if od. The problem Equation (21) can be solved efficiently using any gradient-based algorithms where we can calculate the gradient of 𝖾𝗑𝖺𝖼𝗍 with respect to 𝑳 as

𝖾𝗑𝖺𝖼𝗍(𝒘)vec(𝑳)=pvec(𝑳full-1𝑨𝑴𝑨𝑳full-1,)-pdiag(𝑴)diag((2𝒏)-1)𝑵,

where 𝑴 is a o×o matrix defined in Equation (18), and 𝑵 is a o×d2 matrix whose rows correspond to the observed subsets of the rows of the d2×d2 matrix 𝑳full-1𝑳full-1.

Appendix 1—figure 12 shows the result when the penalized maximum likelihood Equation (21) is applied to the North American wolf dataset with a setting of λ=2.06 (the same value of λ as given in Figure 4) and α=1/w^0, where w^0 is the solution for the ‘constant-w’ model. We can see that the resulting estimated migration surfaces are qualitatively similar to that shown in Figure 4. We also observed similar results between FEEMS and the penalized maximum likelihood Equation (21) across multiple datasets. On the other hand, we found that at the fitted surface the residual variances 1Lkk are not always positive because the constraints are enforced only at the observed nodes. This is problematic because it can cause the model to be ill-defined at the unobserved nodes and make the algorithm numerically unstable. Note that FEEMS avoids this issue by decoupling the residual variance parameter σ2 from the graph-related parameters w. The resulting model Equation (6) also has more resemblance to spatial coalescent model used in EEMS (Petkova et al., 2016), and we thus recommend using FEEMS as a primary method for inferring migration rates.

Jointly estimating the residual variance and edge weights

One simple strategy we have used throughout the paper was to fit σ2 first under a model of homogeneous isolation by distance and prefix the estimated residual variance to the resulting σ^2 for later fits of the effective migration rates. Alternatively, we can consider estimating the unknown residual variance simultaneously with the edge weights, instead of prefixing it from the estimation of the null model—the hope here is to simultaneously correct the model misspecification and allow for improving model fit to the data. To develop the framework for simultaneous estimation of the residual variance and edge weights, let us consider a model that generalizes both Equation (6) and Equation (19), that is,

(22) pCΣ^C𝒲o1(CALAC+Cdiag(n1)Adiag(σ2)AC,p),

where 𝝈2 is a d×1 vector of node specific residual variance parameters, that is, each deme has its own residual parameter σk. If the parameters σk’s are assumed to be the same across nodes, this reduces to the FEEMS model Equation (6) while setting σk=(1Lkk)/2 gives the model Equation (19). Then we solve the following optimization problem

w^,σ^2=argminlwu,σ2>0joint(w,σ2;CΣ^C)+ϕλ,α(w),

where 𝗃𝗈𝗂𝗇𝗍 is the negative log-likelihood function based on Equation (22). Note that the residual variances and edge weights are both searched in the optimization for finding the optimal solutions. To solve the problem, we can use the quasi-newton algorithm for optimizing the objective function.

Appendix 1—figure 13 shows the fitted graphs with different strategies of estimating the residual variances. Appendix 1—figure 13A shows the result when the model has a single residual variance σ2, and Appendix 1—figure 13B shows the result when the residual variances are allowed to vary across nodes. In both cases, estimating the residual variances jointly with the edge weights yields similar and comparable outputs to the default setting of prefixing it from the null model (Figure 4), except that we can further observe reduced effective migration around Queen Elizabeth Islands as shown in Appendix 1—figure 13B. In EEMS, in order to estimate the genetic diversity parameters for every spatial location, which play a similar role as the residual variances in FEEMS, a Voronoi-tessellation prior is placed to encourage sharing of information across adjacent nodes and prevent over-fitting. Similarly, we can place the spatial smooth penalty on the residual variances (i.e. ϕλ,α defined on the variable 𝝈2), but it introduces additional hyperparameters to tune, without substantially improving the model’s fit to the data. In this work, we choose to fit the single residual variance σ2 under the null model and prefix it as a simple but effective strategy with apparent good empirical performance.

Edge versus node parameterization

One of the novel features of FEEMS is its ability to directly fit the edge weights of the graph that best suit the data. This direct edge parameterization may increase the risk of model’s overfitting, but also allows for more flexible estimation of migration histories. Furthermore, as seen in Figure 2 and Appendix 1—figure 2, it has potential to recover anisotropic migration processes. This is in contrast to EEMS wherein every spatial node is assigned an effective migration parameter mk and a migration rate on each edge joining nodes k and is given by the average effective migration wk=(mk+m)/2. Not surprisingly, by assigning each edge to be the average of connected nodes, a form of implicit spatial regularization is imposed because multiple edges connected to the same node would average that node’s parameter value. In some cases, this has the desirable property of imposing an additional degree of similarity across edge weights, but at the same time it also restricts the model’s capacity to capture a richer set of structure present in the data (e.g. Petkova et al., 2016, Supplementary Figure 2). To be concrete, Appendix 1—figure 15 displays two different fits of FEEMS based on edge parameterization (Appendix 1—figure 15A) and node parameterization (Appendix 1—figure 15B), run on a previously published dataset of human genetic variation from Africa (see Peter et al., 2020 for details on the description of the dataset). Running FEEMS with a node-based parameterization is straightforward in our framework—all we have to do is to reparameterize the edge weights by the average effective migration and solve the corresponding optimization problem (Optimization) with respect to 𝒎. It is evident from the results that FEEMS with edge parameterization exhibits subtle correlations that exist between the annotated demes in the figure, whereas node parameterization fails to recover them. We also compare the model fit of FEEMS to the observed genetic distance (Appendix 1—figure 16) and find that edge-based parameterization provides a better fit to the African dataset. Appendix 1—figure 17 further demonstrates that in the coalescent simulations with anisotropic migration, the node parameterization is unable to recover the ground truth of the underlying migration rates even when the nodes are fully observed.

Smooth penalty with 1 norm

FEEMS’s primary optimization objective (see Equation 9) is:

Minimizelwu(w,σ2;CΣ^C)+ϕλ,α(w),

where the spatial smoothness penalty is given by an 2-based penalty function: ϕλ,α(w)=λ2Δlog(eαw1)22. It is well known that an 1-based penalty can lead to a better local adaptive fitting and structural recovery than 2-based penaltyies (Wang et al., 2016), but at the cost of handling non-smooth objective functions that are often computationally more challenging. In a spatial genetic dataset, one major challenge is to deal with the relatively sparse sampling design where there are many unobserved nodes on the graph. In this statistically challenging scenario, we found that an 2-based penalty allows for more accurate and reliable estimation of the geographic features.

Specifically, writing ϕλ,α1(w)=λΔlog(eαw1)1, we considered the alternate following composite objective function:

(23) (w,σ2;CΣ^C)+ϕλ,α1(w).

To solve Equation (23), we apply linearized alternating direction method of multipliers (ADMM) (Boyd, 2010), a variant of the standard ADMM algorithm, that iteratively optimizes the augmented Lagrangian over the primal and dual variables. The derivation of the algorithm is a standard calculation so we omit the detailed description of the algorithm. As opposed to the common belief about the effectiveness of the 1 norm for structural recovery, the recovered graph of FEEMS using 1-based smooth penalty shows less accurate reconstruction of the migration patterns, especially when the sampling design has many locations with missing data on the graph (Appendix 1—figure 18A, Appendix 1—figure 19H). We can see that the 1-based penalty function is not able to accurately estimate edge weights at regions with little data, partially due to its local adaptation, in contrast to the 2-based method that considers regularization more globally. This suggests that in order to use the 1 penalty ϕλ,α1(𝒘) in the presence of many missing nodes, one may need an additional degree of regularization that encourages global smoothness of the graph’s edge weights, such as a combination of ϕλ,α1(w) and ϕλ,α(𝒘) (in the same spirit as elastic net [Zou and Hastie, 2005]), or ϕλ,α1(𝒘) on top of node-based parameterization (see Appendix 1—figure 18B).

Coalescent simulations with weak migration

In Figure 2, we evaluated FEEMS by applying it to ‘out-of-model’ coalescent simulations. In these simulations, we generated genotype data under a coalescent model with structured meta-populations organized on a spatial triangular lattice. In a relatively ‘strong’ heterogeneous migration scenario (Figure 2D,E,F), we set the coalescent migration rate to be an order of magnitude lower (10-fold) in the center of the spatial grid than on the left and right regions, emulating a depression in gene-flow caused, for example, by a mountain range or body of water. The variation in migration rates should create a spatially varying covariance structure in the genetic variation data. To get a sense of the level of genetic divergence implied by this simulation setting, we visualized Wright’s fixation index (FST, Patterson’s estimator [Patterson et al., 2012]) plotted against the geographic distance between nodes (Appendix 1—figure 20). We see in the strong heterogeneous migration simulation there is a clear signal of two clusters of data points (Appendix 1—figure 20B). These clusters correspond to pairwise FST comparisons of two nodes on the same side of the central depression in gene flow, where gene-flow roughly follows a homogeneous ‘isolation-by-distance’ like pattern, or two nodes across the central depression where gene-flow is reduced, hence increasing the expected FST between such nodes.

While simulating this strong reduction of gene-flow provides an illustrative and clear example where FEEMS has a lot of signal for accurate inference, we wanted to understand the qualitative performance of FEEMS in an less idealized scenario with weaker signal. To this end, we performed coalescent simulations with only a 25% reduction of gene-flow in the center of the habitat (Appendix 1—figure 21). In Appendix 1—figure 21A, when all the nodes are observed on the spatial graph, FEEMS is still able to detect this subtle reduction of gene-flow. While FEEMS is able to detect this signal, there remain particularly erroneous estimates among the lower than average edge weights, implying the fit could benefit from additional smoothing by increasing the level regularization on the smoothness penalty. In contrast to the strong heterogeneous migration simulations, we see that the pairwise FST in this weak migration scenario does not obviously show a ‘clustering’ like effect in the data (Appendix 1—figure 20A). The average FST between all pairs of demes is approximately three times lower (mean FST=.1175 for the weak heterogeneity simulation versus mean FST = 0.3411 for the strong heterogeneity simulation). When the nodes are sparsely observed on the graph in this weak migration simulation, we see that the FEEMS output is overly smooth (Appendix 1—figure 21B). In the absence of data and thus a weak signal for spatial variation in migration, a smooth visualization is arguably a sensible outcome given the regularization acts like a prior distribution favoring spatial homogeneity in levels of effective migration.

In practice, weak population structure can be more accurately dissected when increasing the number of informative SNPs included in the analysis (Novembre and Peter, 2016). In conjunction with running FEEMS, we recommend for users to create exploratory visualizations such as variograms and PCA bi-plots to assess the level of population structure in their data, and to consider the number of SNPs used in the analysis.

Appendix 1—figure 1
Visualization of grid construction and node assignment: (A) Map of sample coordinates (black points) from a dataset of gray wolves from North America.

The input to FEEMS are latitude and longitude coordinates as well as genotype data for each sample. (B) Map of sample coordinates with an example dense spatial grid. The nodes of the grid represent sub-populations and the edges represent local gene-flow between adjacent sub-populations. (C) Individuals are assigned to nearby nodes (sub-populations) and summary statistics (e.g. allele frequencies) are computed for each observed location.

Appendix 1—figure 2
Application of FEEMS to an extended set of coalescent simulations: We display an extended set of coalescent simulations with multiple migration scenarios and sampling designs.

The sample sizes across the grid are represented by the size of the gray dots at each node. The migration rates are obtained by solving FEEMS objective function Equation 9 where the the smoothness parameter λ was selected using leave-one-out cross-validation. (A, F, K) display the ground truth of the underlying migration rates. (B, G, L) shows simulations where there is no missing data on the graph. (C, H, M) shows simulations with sparse observations and nodes missing at random. (D, I, N) shows simulations of biased sampling where there are no samples from the center of the simulated habitat. (E, J, O) shows simulations of biased sampling where there are only samples on the right side of the habitat.

Appendix 1—figure 3
Application of FEEMS to a heterogeneous migration scenario with a ‘missing at random’ sampling design: We run FEEMS on coalescent simulation with a non-homogeneous process while varying hyperparameters λ (rows) and α (columns).

We randomly sample individuals for 20% of nodes. When λ grows, the fitted graph becomes overall smoother, whereas α effectively controls the degree of similarity among low migration rates.

Appendix 1—figure 4
Application of FEEMS to an anisotropic migration scenario with a ‘missing at random’ sampling design: We run FEEMS on coalescent simulation with an anisotropic process while varying hyperparameters λ (rows) and α (columns).

We randomly sample individuals for 20% of nodes. When λ grows, the fitted graph becomes overall smoother, whereas α effectively controls the degree of similarity among low migration rates.

Appendix 1—figure 5
SNP and individual quality control.

(A) displays a visualization of the sample site frequency spectrum. Specifically, we display a histogram of minor allele frequencies across all SNPs. We see a relatively uniform histogram which reflects the ascertainment of common SNPs on the array that was designed to genotype gray wolf samples. (B) visualization of allele frequencies plotted against genotype frequencies. Each point represents a different SNP and the colors represent the three possible genotype values. The black dashed lines display the expectation as predicted from a simple binomial sampling model i.e., Hardy-Weinberg equilibrium. (C) displays a histogram of the missingness fraction per SNP. We observe the missingness tends to be relatively low for each SNP. (D) displays a histogram of the missingness fraction per sample. Generally, the missingness tends to be low for each sample.

Appendix 1—figure 6
Comparing predictions of observed genetic distances: We display different predictions of observed genetic distances using geographic distance or the fitted genetic distance output by FEEMS.

(A) The x-axis displays the geographic distance between two individuals, as measured by the great circle distance (haversine distance). The y-axis displays the squared Euclidean distance between two individuals averaged over all SNPs. (B–D) The x-axis displays the fitted genetic distance as predicted by the FEEMS model and y-axis displays the squared Euclidean distance between two sub-populations averaged over all SNPs. For (B–D), we display the fit of λ getting subsequently smaller, λ=100,2.06,0.04 (the same values of λ used in Figure 3A,B,C), and as expected the fit appears better because we tolerate more complex surfaces and we are not evaluating the fit on out-of-sample data.

Appendix 1—figure 7
Summary of top axes of genotypic variation: We display a visual summary of Principal Components Analysis (PCA) applied to the normalized genotype matrix from the North American gray wolf dataset.

(A–D) displays PC bi-plots of the top seven PCs plotted against each other. The colors represent predefined ecotypes defined in Schweizer et al., 2016. We can see that the top PCs delineate these predefined ecotypes. (E) shows a ‘scree’ plot with the proportion of variance explained for each of the top 50 PCs. As expected by genetic data (Patterson et al., 2006), the eigenvalues of the genotype matrix tend to be spread over many PCs.

Appendix 1—figure 8
Relationship between top axes of genetic variation and latitude: In each sub-panel, we plot the PC value against latitude for each sample in gray the wolf dataset.

We see many of the top PCs are significantly correlated with latitude as tested by linear regression.

Appendix 1—figure 9
Relationship between top axes of genetic variation and longitude: In each sub-panel, we plot the PC value against longitude for each sample in the gray wolf dataset.

We see many of the top PCs are significantly correlated with longitude as tested by linear regression.

Appendix 1—figure 10
Summary of ADMIXTURE results: (A–G) Visualization of ADMIXTURE results for K=2 to K=8.

We display admixture fractions for each sample as colored slices of the pie chart on the map. For each K, we ran five replicate runs of ADMIXTURE and in this visualization, we display the solution that achieves the highest likelihood amongst the replicates. The ADMIXTURE results qualitatively reveal a spatial signal in the data as admixture fractions tend to be spatially clustered.

Appendix 1—figure 11
Application of EEMS to the North American gray wolf dataset: We display a visualization of EEMS applied to the North American gray wolf dataset.

The more orange colors represent lower than average effective migration on the log-scale and the more blue colors represent higher than average effective migration on the log-scale. The results of EEMS are qualitatively similar to FEEMS when lower regularization penalties are applied.

Appendix 1—figure 12
Application of FEEMS on the North American gray wolf dataset with an exact likelihood model: We display the fit of FEEMS based on the formulation Equation (21) to the North American gray wolf dataset.

This fit corresponds to a setting of tuning parameters at λ=2.06 and α=1/w^0. Additionally, we set the lower bound of the edge weights to 𝒍=10-6, to ensure that the diagonal elements of 𝑳 does not become too small—this has an implicit effect on Lkk, preventing it from blowing up at unobserved nodes. The more orange colors represent lower than average effective migration on the log-scale and the more blue colors represent higher than average effective migration on the log-scale. Visually, the result is comparable to that of the FEEMS fit (Figure 4) based on the approximate formulation (Optimization).

Appendix 1—figure 13
Application of FEEMS on the North American gray wolf dataset with joint estimation of the residual variances and graph’s edge weights: We show visualizations of fits of FEEMS to the North American gray wolf dataset when the residual variances and edge weights of the graph are jointly estimated.

Both fits correspond to a setting of tuning parameters at λ=2.06 and α=1/w^0. (A) displays the estimated effective migration surfaces where every deme shares a single residual parameter σ. (B) displays the estimated effective migration surfaces where each node has its own residual parameter σk. Both approaches yield similar results to the procedure that prefixes σ from the homogeneous isolation by distance model (Figure 4). The node-specific residual parameters may allow for more flexible graphs to be fitted, and we can further observe reduced effective migration around C (Queen Elizabeth Islands) in (B).

Appendix 1—figure 14
Relationship between fitted versus observed genetic dissimilarities on the North American gray wolf dataset: We display scatter plots of fitted genetic distance versus observed genetic distance from FEEMS fits on the gray wolf dataset.

(A) Corresponds to the result shown in Figure 4. (B) Corresponds to the result shown in Appendix 1—figure 13B. The x-axis displays the fitted genetic distance as predicted by the FEEMS model and y-axis displays the squared Euclidean distance between two sub-populations averaged over all SNPs. The simple linear regression fit is shown in orange dashed lines and R2 is given.

Appendix 1—figure 15
Application of FEEMS to a dataset of human genetic variation from Africa with different parameterization: We display visualizations of FEEMS to a dataset of human genetic variation from Africa with different parameterization of the graph’s edge weights.

See Peter et al., 2020 for the description of the dataset. (A) displays the recovered graph under the edge parameterization. (B) displays the recovered graph under the node parameterization. Both parameterization have their own regularization parameters λ and α, but these parameters are not on the same scale. We set λ=0.2 and α=0.02 for the node parameterization which is seen to yield similar results to those in Peter et al., 2020. For the edge parameterization, we set λ=0.5 and α=0.05 so that the resulting graph reveals similar geographic structure to the node parameterization. We also set the lower bound 𝒍=10-6. From the plots, it is worth noting two important distinctions: (1) We see the migration surfaces shown in (B) recover sharper edge features while the migration surfaces in (A) are overall smoother. This is attributed to the fact that node parameterization has its own additional regularization effect on the edge weights, and in order to achieve similar degree of regularization strength for the edge parameterization, it needs a higher regularization parameters, which results in more blurring edges than the node parameterization. (2) When measuring correlation of the estimated allele frequencies among nodes, we find that Deme B is the node with the second highest correlation to Deme A, whereas Deme C (and nearby demes) is not as much correlated to Deme A compared to Deme B. Panel (A) reflects this feature by exhibiting a corridor between Deme A and Deme B and reduced gene-flow beneath that corridor. This reduced gene-flow disappears in (B), even if the regularization parameters are varied over a range of values. Additionally, Deme D is most highly correlated to Deme E, F, and G, and this is implicated by a long-range corridor connecting those demes appearing in Panel (A) while not shown in (B). These results suggest that the form of the node parameterization is perhaps too strong and in this case limits the model’s ability to capture desirable geographic features that are subtle to detect.

Appendix 1—figure 16
Relationship between fitted versus observed genetic dissimilarities on a dataset of human genetic variation from Africa: We display scatter plots of fitted genetic distance versus observed genetic distance from FEEMS fits on the African dataset.

(A) Corresponds to the result shown in Appendix 1—figure 15A. (B) Corresponds to the result shown in Appendix 1—figure 15B. The x-axis displays the fitted genetic distance as predicted by the FEEMS model and y-axis displays the squared Euclidean distance between two sub-populations averaged over all SNPs. The simple linear regression fit is shown in orange dashed lines and R2 is given.

Appendix 1—figure 17
Application of FEEMS based on node parameterization to an extended set of coalescent simulations: We display an extended set of coalescent simulations with the same migration scenarios and sampling designs as Appendix 1—figure 2.

The sample sizes across the grid are represented by the size of the gray dots at each node. The migration rates are obtained by solving the FEEMS objective function (Optimization) with node parameterization where the regularization parameters are specified at λ=5 and α=0.01. (A, F, K) display the ground truth of the underlying migration rates. (B, G, L) shows simulations where there is no missing data on the graph. (C, H, M) shows simulations with sparse observations and nodes missing at random. (D, I, N) shows simulations of biased sampling where there are no samples from the center of the simulated habitat. (E, J, O) shows simulations of biased sampling where there are only samples on the right side of the habitat.

Appendix 1—figure 18
Application of 1-norm-based FEEMS to a dataset of human genetic variation from Africa: We display visualizations of FEEMS to a dataset of human genetic variation from Africa with the 1-based penalty function.

See Peter et al., 2020 for the description of the dataset. (A) displays the recovered graph under the edge parameterization with 1 norm based penalty where the regularization parameters are specified at λ=0.05 and α=5. (B) displays the recovered graph under the node parameterization with 1 norm-based penalty where the regularization parameters are specified at λ=0.05 and α=1. To minimize the objective Equation (23), linearized ADMM is applied with 20,000 number of iterations. The lower bound is set to be 𝒍=10-6 for both parameterizations. Note that due to the high degrees of missingness, the estimated effective migration surfaces using solely 1-based penalty exhibit many likely artifacts (e.g. high migration edges forming long paths, seen in panel A) unless an additional regularization is added to encourage global smoothness of the edge weights, such as a combination of 1 norm penalty function and node parameterization as shown in (B).

Appendix 1—figure 19
Application of 1-norm-based FEEMS to an extended set of coalescent simulations: We display an extended set of coalescent simulations with the same migration scenarios and sampling designs as Appendix 1—figure 2.

The sample sizes across the grid are represented by the size of the gray dots at each node. The migration rates are obtained by solving 1 norm based FEEMS objective Equation (23) where the regularization parameters are specified at λ=5 and α=0.01. (A, F, K) display the ground truth of the underlying migration rates. (B, G, L) shows simulations where there is no missing data on the graph. (C, H, M) shows simulations with sparse observations and nodes missing at random. (D, I, N) shows simulations of biased sampling where there are no samples from the center of the simulated habitat. (E, J, O) shows simulations of biased sampling where there are only samples on the right side of the habitat.

Appendix 1—figure 20
Comparing pairwise Fst between strong and weak heterogeneous migration coalescent simulations: We visualize of the relationship between geographic distance and estimated pairwise Fst (genetic distance) between nodes on the spatial grid for a weak heterogeneous migration simulation (A) and a strong heterogeneous migration (B).

As expected the average Fst is lower for the weak migration setting, and we observe a clear clustering like effect in the data for the strong heterogeneous migration simulation. This strong clustering effect can be attributed to pairwise comparisons of nodes across the region of reduced gene flow. Distances between nodes were set to one in the simulation, and so the units of geographic distance here are in units of the inter-node distance.

Appendix 1—figure 21
Applications of FEEMS to weak migration coalescent simulations: We visualize the FEEMS fit to coalescent simulations with weak heterogeneous migration.

The coalescent migration rate in the center of the habitat is set to be 25% lower than the left or right regions. Note that the color-scale limits are set to 10-.15 and 10.15, respectively. The top panel shows the fit when all the nodes of the spatial graph are observed, whereas the bottom panel shows the fit when a sparse subset of nodes are observed. We see that FEEMS can still detect a signal of heterogeneity by displaying reduced gene-flow in the center of the habitat in the top panel. When we observe only a few nodes, in this weak migration setting, the FEEMS visualization looks overly smooth.

Appendix 1—figure 22
EEMS fitted versus observed genetic dissimilarities from the North American gray wolf dataset: We visualize fitted versus observed genetic dissimilarities corresponding to the EEMS visualization run on the North American gray wolf dataset in Figure 4.

EEMS was run on a sparse grid with 307 nodes due to long run-times on the dense grid. Generally there is good concordance between the fitted and observed dissimilarities except for a small set of points whose fitted genetic dissimilarity over-estimates the observed dissimilarity, implying a relatively poorly fit for these points. Note, we do not see these poorly fit points in visualizations of the fitted versus observed distances when using FEEMS (see Appendix 1—figure 6).

Data availability

Genotyping data can be found at https://doi.org/10.5061/dryad.c9b25 and stored in the FEEMS python package at https://github.com/Novembrelab/feems (copy archived at https://archive.softwareheritage.org/swh:1:rev:2df82f92ba690f5fd98aee6612b155d973ffb12d).

The following previously published data sets were used
    1. Schweizer RM
    2. vonHoldt JC
    3. R
    4. BM
    5. Harrigan
    6. Knowles
    7. Musiani M
    8. Coltman D
    9. Novembre J
    10. Wayne RK
    (2016) Dryad Digital Repository
    Genetic subdivision and candidate genes under selection in North American grey wolves.
    https://doi.org/10.5061/dryad.c9b25

References

  1. Conference
    1. Kalofolias V
    (2016)
    How to learn a graph from smooth signals
    Proceedings of the 19th International Conference on Artificial Intelligence and Statistics. pp. 920–929.
  2. Conference
    1. Kimura M
    (1953)
    Stepping stone model of population
    Annual Report of the National Institute of Genetics Japan. pp. 62–63.
  3. Book
    1. Koller D
    2. Friedman N
    (2009)
    Probabilistic Graphical Models: Principles and Techniques
    MIT Press.
  4. Book
    1. Lauritzen SL
    (1996)
    Graphical Models
    Clarendon Press.
  5. Book
    1. Malécot G
    (1948)
    Les mathématiques de L’hérédité
    Paris, France: Masson et cie.
    1. McCullagh P
    (2009)
    Marginal likelihood for distance matrices
    Statistica Sinica 19:631–649.
  6. Book
    1. Nocedal J
    2. Wright S
    (2006)
    Numerical Optimization
    Springer Science & Business Media.
  7. Book
    1. Petkova D
    (2013)
    Inferring Effective Migration From Geographically Indexed Genetic Data
    Chicago, United States: The University of Chicago Press.
  8. Book
    1. Rue H
    2. Held L
    (2005)
    Gaussian Markov Random Fields: Theory and Applications
    CRC press.
    1. Vishnoi NK
    (2013) Lx = b
    Foundations and Trends in Theoretical Computer Science 8:1–141.
    https://doi.org/10.1561/0400000054
    1. Wang Y-X
    2. Sharpnack J
    3. Smola AJ
    4. Tibshirani RJ
    (2016)
    Trend filtering on graphs
    The Journal of Machine Learning Research 17:3651–3691.

Decision letter

  1. George H Perry
    Senior and Reviewing Editor; Pennsylvania State University, United States
  2. Isabel Alves
    Reviewer; University of Nantes, France
  3. Wesley Tansey
    Reviewer

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

The authors of the manuscript present a new implementation of the previously developed statistical method called "Estimating Effective Migration Surfaces", which displays on geographical map regions of low or high effective migration under a broad model of isolation by distance. In this new implementation migration surfaces are estimated under a penalized-likelihood approach coupled with optimization instead of MCMC leading to faster running times. The new implementation facilitates faster running times to make its usage computationally possible for a wider range of research groups and likely be applied to an even larger number of species/populations.

Decision letter after peer review:

Thank you for submitting your article "Fast and Flexible Estimation of Effective Migration Surfaces" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by George Perry as the Senior and Reviewing Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Isabel Alves (Reviewer #1); Wesley Tansey (Reviewer #2).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

As the editors have judged that your manuscript is of interest, but as described below that additional analyses are required before it is published, we would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). First, because many researchers have temporarily lost access to the labs, we will give authors as much time as they need to submit revised manuscripts. We are also offering, if you choose, to post the manuscript to bioRxiv (if it is not already there) along with this decision letter and a formal designation that the manuscript is "in revision at eLife". Please let us know if you would like to pursue this option. (If your work is more suitable for medRxiv, you will need to post the preprint yourself, as the mechanisms for us to do so are still in development.)

Summary:

The authors of the manuscript present a new implementation of the previously developed statistical method called "Estimating Effective Migration Surfaces", which displays on a geographical map regions of low or high effective migration under a broad model of isolation by distance. In this new implementation migration surfaces are estimated under a penalized-likelihood approach coupled with optimization instead of MCMC leading to faster running times. The new implementation appears very promising as faster running times will make its usage computationally possible for a wider range of research groups and likely be applied to an even larger number of species/populations. Overall, we value the approach for its pragmatism but felt that it falls short at the very end by failing to provide any quantitative, objective way to choose the hyperparameters, which needs to be addressed as per the essential revisions detailed below.

Essential revisions:

1. The authors need to provide principled (or at least reproducible) ways to select the hyperparameters. Specific reviewer comments include:

"This is a major benefit of the L1 penalty. You can use BIC as the model selection criterion in the L1 case since the degrees of freedom are well-described. In the squared L2 case, it's not really possible. The authors discuss the issue and note LOOCV did not produce stable results, but they do not provide any data or examples. A more thorough investigation of hyperparameter settings is needed along with a recommendation that does not rely on biologists' subjective preference of the results on each dataset."

"FEEMS outcomes are very sensitive to user-based settings such as grid density and tuning parameters, as well as to aspects of the real data (eg. sampling design) that may result in an arbitrary choice of the outcome and lead to over-interpretations. I know the authors recommend to explore several combinations of regularization parameters and then compare FEEMS results with clustering/differentiation patterns based on approaches like ADMIXTURE or FST distances in order to support the results, nevertheless it is still difficult to grasp what's the best strategy to assess if a fitted graph is over-fitting the observed data or instead is pointing out to a real area of, let's say, low effective migration rate. Sentences like: "…, while setting up the tuning parameter ɑ to a value that we found that worked for multiple data applications.…" (lines: 225-226) or "it is helpful to look more closely at particular solutions that find balance between spatial homogeneity and complexity…" (lines: 243-245) are confusing and make difficult the choice of the final regularization parameters."

"The grid design is another arbitrary aspect of the method whose influence on the identification of regions of low or high migration isn't clear. Imagining one has the computational power to construct a very dense grid, is it worth doing it once there is observed data in 1% of the nodes? Is there a good relationship between density and number of sampled points? Does it affect the outcome?"

"I think the points above would be clearer if the authors would provide for instance, step-by-step guidelines to help future users of FEEMS and referring to specific examples in the manuscript in order to more clearly justify their parameter choice (eg ɑ = 50 line 226)."

2. Please clarify the modeling decision and its comparison to the L1 approach. For instance, isn't smoothing over nodes vs edges really just the same thing with different penalties? The authors form a lifted graph, where edges are now nodes and they penalize differences between neighboring edges. In the L1 penalty case with a squared error loss, the fused lasso / total variation penalty on neighboring edges is equivalent to linear trend filtering on the nodes. The choice then of linear trend filtering could have been replaced with a higher order trend filtering step to achieve the smoothness that the authors seem to say is lacking in the L1 model.

3. Are the data points we observe actually sampled at random? Is some sort of latent confounding likely? For example, maybe wolves migrate based on the season and the scientists collecting the data only look in one spot in one particular time of year?

4. Please clarify the simulation results in Supp Fig19, panels I and J. Without any data points in the orange regions for panel I, the model somehow infers that there is a band of different edge weights. How? In the 1d case, it's as if someone showed you: [5, 5, 5, missing, missing, missing, 5 ,5 5], and you come back and told me [5, 5, 5, -3, -3, -3, 5, 5, 5]. Is this possible?

5. At present, there is no real quantitative assessment of how good the FEEMS solutions are relative to the EEMS solution. This should be provided.

6. It would be useful to provide another example of a heterogeneous migration scenario where the reduction in migration is less than one order of magnitude in order to give an idea to the user of how the method performs in a less heterogeneous scenario (ie the lower bound).

https://doi.org/10.7554/eLife.61927.sa1

Author response

Essential revisions:

1. The authors need to provide principled (or at least reproducible) ways to select the hyperparameters. Specific reviewer comments include:

"This is a major benefit of the L1 penalty. You can use BIC as the model selection criterion in the L1 case since the degrees of freedom are well-described. In the squared L2 case, it's not really possible. The authors discuss the issue and note LOOCV did not produce stable results, but they do not provide any data or examples. A more thorough investigation of hyperparameter settings is needed along with a recommendation that does not rely on biologists' subjective preference of the results on each dataset."

"FEEMS outcomes are very sensitive to user-based settings such as grid density and tuning parameters, as well as to aspects of the real data (eg. sampling design) that may result in an arbitrary choice of the outcome and lead to over-interpretations. I know the authors recommend to explore several combinations of regularization parameters and then compare FEEMS results with clustering/differentiation patterns based on approaches like ADMIXTURE or FST distances in order to support the results, nevertheless it is still difficult to grasp what's the best strategy to assess if a fitted graph is over-fitting the observed data or instead is pointing out to a real area of, let's say, low effective migration rate. Sentences like: "…, while setting up the tuning parameter ɑ to a value that we found that worked for multiple data applications.…" (lines: 225-226) or "it is helpful to look more closely at particular solutions that find balance between spatial homogeneity and complexity…" (lines: 243-245) are confusing and make difficult the choice of the final regularization parameters."

"The grid design is another arbitrary aspect of the method whose influence on the identification of regions of low or high migration isn't clear. Imagining one has the computational power to construct a very dense grid, is it worth doing it once there is observed data in 1% of the nodes? Is there a good relationship between density and number of sampled points? Does it affect the outcome?"

"I think the points above would be clearer if the authors would provide for instance, step-by-step guidelines to help future users of FEEMS and referring to specific examples in the manuscript in order to more clearly justify their parameter choice (eg ɑ = 50 line 226)."

We thank the reviewers for these helpful comments in regards to selecting the hyper-parameters of the penalty and agree that an automated selection procedure would help improve the interpretability and usability of FEEMS. To this end, we have made a number of updates to our modeling approach that have allowed for fully automated hyper-parameter selection and have proven to work well in practice. These developments were yielded through a simple but effective update of the parameterization of our penalty that allows for a cross-validation approach over just the smoothness parameter lambda alone. Specifically, we utilize the solution of the edge weights fitted under a homogenous migration model to penalize differences in neighboring edge weights on both the linear and log scale relative to a homogeneous fitted parameter which we call w0. This natural parameterization of the penalty and pre-estimation of w0 under a simple homogenous null model, allowed us to effectively remove the alpha parameter from the original penalty, allowing an one-dimensional cross-validation algorithm for selecting lambda in a computationally efficient and reliable manner. For more details please see the updated “Overview of FEEMS” in the Results section and “Penalty description” in the Materials and methods section.

With this new penalty in hand we used leave-one-out cross-validation to select the smoothness parameter lambda. In the cross-validation algorithm, for each grid value of lambda, we held-out an individual observed node (population) on the graph and then predicted underlying allele frequencies at these held-out nodes under our fitted spatial model from the rest of the training-set nodes (see the new section “Leave-one-out cross-validation to select tuning parameters” in Materials and methods for details). We found leave-one-out cross-validation over lambda to provide satisfactory results that recovered true migration histories in coalescent simulations and aligned with biological expectations in real datasets.

We have updated all of the text with descriptions of the new penalty and have reanalyzed and reproduced the figures using leave-one-out cross-validation to select lambda. This greatly simplifies the text and we hope this helps to alleviate the comments posed by the reviewers in regards to hyper-parameter selection. Figure 3 now shows fitted FEEMS visualizations across the grid points of lambda that were used in leave-one-out cross-validation. We added a new panel Figure 3E which shows the cross-validation error for the full grid and highlights the visualized maps for a subset of lambda values. In Figure 4, we now display the solution that achieves the minimum cross-validation error. All coalescent simulation figures have been updated as well. In general, the interpretation and results have not changed using this new penalty and cross-validation procedure but the ease of use and clarity of FEEMS has greatly improved.

We also thank the reviewers for the suggestion on using the BIC for model selection under the L1 penalty. As mentioned above, for our new penalty – which is still the L2 distance between migration weights on neighboring edges, we have found leave-one-out cross-validation to perform well. Because cross-validation, in principle, can work for both the L1 and L2 penalties we prefer to use a method that works more "universally" for any penalty whether it induces exact sparsity or not. We also prefer the L2 penalty for other statistical and computational reasons and expand upon this point in our response in the next section.

In terms of concerns about the grid density, we do not have new solutions for this problem which is also a caveat in the original EEMS method (see Petkova et al. 2016 discussion); however we highlight this issue more prominently with a new paragraph in the discussion.

2. Please clarify the modeling decision and its comparison to the L1 approach. For instance, isn't smoothing over nodes vs edges really just the same thing with different penalties? The authors form a lifted graph, where edges are now nodes and they penalize differences between neighboring edges. In the L1 penalty case with a squared error loss, the fused lasso / total variation penalty on neighboring edges is equivalent to linear trend filtering on the nodes. The choice then of linear trend filtering could have been replaced with a higher order trend filtering step to achieve the smoothness that the authors seem to say is lacking in the L1 model.

In a sense, yes, “smoothing over nodes vs edges is the same thing with different penalties”; however the penalties differ in key ways. In the original parameterization of EEMS, each node was given a parameter and the edge weights were deterministically computed as the average of adjacent connected nodes. While this node-level parameterization reduces the number of parameters needed to be estimated, the edge-level parameterization has two advantages in our view:

1. Each edge is free to take a unique value and that allows for a wider range of anisotropic migration scenarios to be modeled (e.g. spatially homogeneous anisotropy as in Figure 2 right hand column). That was the main driver for the decision for this new smoothing scheme.

2. By assigning each edge to be the average of connected nodes, a form of implicit spatial regularization is imposed because multiple edges connected to the same node would average that node’s parameter value. We found it more natural to separate out the regularization from the parameterization of the model. This preference led us to adding a smoothness penalty on edge-level parameters.

We thank the reviewer for the suggestion of using higher order trend filtering. We tested a L1 penalization approach on the edge weights and find it often fails to give satisfactory results (e.g. Supplementary Materials “Smooth penalty with L1 norm’’ and Supplementary Figure 18), primarily because the L1 penalty is too locally adaptive to regions with many unobserved locations. In the regime where there is a high degree of missingness in the graph, the global consideration of smoothing the unobserved locations seems to be necessary and this is the primary driver for using the L2 penalty. We believe that using the higher order trend filtering with the L1 penalty may suffer a similar issue. The smoothness of the L2 penalty also allowed us to employ a quasi-newton algorithm for optimizing our objective function which decreased our runtime more than 10 fold when compared to first order methods such as proximal gradient descent and ADMM when using the L1 formulation of our objective function. In addition we were able to utilize a widely used and tested implementation of L-BFGS in scipy which worked well out of the box and had few algorithmic parameters to tune. Given we observed satisfactory results with the L2 penalty in addition to the fast convergence and runtime of the quasi-newton algorithm we preferred it over the L1 penalization approach.

3. Are the data points we observe actually sampled at random? Is some sort of latent confounding likely? For example, maybe wolves migrate based on the season and the scientists collecting the data only look in one spot in one particular time of year?

For (a), in our model we treat the geographic locations of each sample as fixed but the distribution of genetic data as random. We find the method is relatively robust to non-random sampling, as long as it is not so sparse as to lose the key signals of differentiation in the data (e.g. Supplementary Figure 2).

Regarding (b), we have expanded the discussion with a paragraph that addresses how a form of confounding between seasonal migration and the sample collection process could be problematic in some datasets. Specifically, we discuss variation in the wolf migratory behavior to illustrate the point, and the suggestion helped us add nuance to our discussion of the results.

4. Please clarify the simulation results in Supp Fig19, panels I and J. Without any data points in the orange regions for panel I, the model somehow infers that there is a band of different edge weights. How? In the 1d case, it's as if someone showed you: [5, 5, 5, missing, missing, missing, 5 ,5 5], and you come back and told me [5, 5, 5, -3, -3, -3, 5, 5, 5]. Is this possible?

It is possible – and there are a few ways to understand where the signal for the inference derives from. (1) Consider that in sampled regions one can “learn” a relationship between geographic and genetic distance, and then with that in hand recognize that the observed covariances across a gap of unobserved locations are too large or too small relative to what is seen in the observed data. Regarding the reviewer’s analogy, let’s define a spatial position of each node as its index in the 1-D array. Let's further suppose we collected a dataset of samples from nodes 1,2,3 and 7,8,9. First, assume we see that the observed allele frequency covariance within pairs of nodes 1,2,3 and 7,8,9 decays some constant amount for every unit of geographic distance. Then, if the covariance between pairs nodes 1,2,3 and 7,8,9 is lower than what would be expected given the geographic distance between them, it would provide a signal that the migration rates should be higher for edges connecting 1,2,3 and 7,8,9 and lower for edges connecting 4,5,6. In contrast, if pairwise covariances between all of the observed nodes decayed over geographic space at a constant level, we would infer homogenous migration rates for edges connecting both the missing and non-missing nodes. Our likelihood uses the migration rate parameters to match expected covariances at all nodes to observed covariances, whereas our penalty encourages the migration rates to be smooth which helps inference in regions with high levels of missing nodes. As a simple one-locus example, consider the allele frequency data observed at the same 1-D array of 9 populations: [0.02,0.01,0.03, missing, missing, missing, 0.99,0.98,0.985]. From such allele frequency data, one can infer migration rates are unlikely to be homogeneous, and moreover, likely are lower in the region with no observed data. (2) In a population genetic model, the key determinant of variation in genetic data observed today is where and when the ancestors of the sampled data “coalesce” (i.e. have common ancestry). The genetic ancestors of the present sample can occupy locations where there is no data today, and when they coalesce with each other will be impacted by local migration rates at unsampled locations. That is to say, migration rates at locations where there are no samples today can still impact the genetic data observed. Surprisingly, this implies that, to some extent, one can learn about migration rates even outside the convex hull of the sampled points.

5. At present, there is no real quantitative assessment of how good the FEEMS solutions are relative to the EEMS solution. This should be provided.

We now include a new supplemental figure 22 with the observed vs fitted dissimilarities output by EEMS when applied to the North American gray wolf dataset. This can be compared to the analogous results for FEEMS already presented in supplemental figure 14. The results show that EEMS fits the data relatively well, though there is a collection of points that seem to be poorly fit as illustrated by a systematically higher fitted dissimilarity to what is observed. FEEMS, in contrast, does not seem to have the same cluster of poor fitting samples.

6. It would be useful to provide another example of a heterogeneous migration scenario where the reduction in migration is less than one order of magnitude in order to give an idea to the user of how the method performs in a less heterogeneous scenario (ie the lower bound).

We now include such an example. To assess the performance of FEEMS in a less heterogeneous migration scenario, we applied FEEMS to coalescent simulations where the migration in the center of the habitat (spatial grid) was only reduced by 25% relative to the edges. In Supplementary Figure 21, we can see that FEEMS is still able to recover this reduction of gene-flow but the output visualization is slightly noisier than the strong heterogeneous migration simulations in Figure 2. Please refer to a new supplementary section titled “Coalescent simulations with weak migration” for a detailed discussion of the results.

https://doi.org/10.7554/eLife.61927.sa2

Article and author information

Author details

  1. Joseph Marcus

    Department of Human Genetics, University of Chicago, Chicago, United States
    Contribution
    Conceptualization, Software, Formal analysis, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Contributed equally with
    Wooseok Ha
    For correspondence
    jhmarcus@uchicago.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0923-9881
  2. Wooseok Ha

    Department of Statistics, University of California, Berkeley, Berkeley, United States
    Contribution
    Software, Formal analysis, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Contributed equally with
    Joseph Marcus
    For correspondence
    haywse@berkeley.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9069-854X
  3. Rina Foygel Barber

    Department of Statistics, University of Chicago, Chicago, United States
    Contribution
    Formal analysis, Supervision, Methodology, Writing - review and editing
    Contributed equally with
    John Novembre
    For correspondence
    rina@uchicago.edu
    Competing interests
    No competing interests declared
  4. John Novembre

    1. Department of Human Genetics, University of Chicago, Chicago, United States
    2. Department of Ecology and Evolution, University of Chicago, Chicago, United States
    Contribution
    Conceptualization, Formal analysis, Supervision, Visualization, Writing - review and editing
    Contributed equally with
    Rina Foygel Barber
    For correspondence
    jnovembre@uchicago.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5345-0214

Funding

National Science Foundation (DGE-1746045)

  • Joseph Marcus

National Institute of General Medical Sciences (T32GM007197)

  • Joseph Marcus

National Institute of General Medical Sciences (R01GM132383)

  • John Novembre

National Science Foundation (TRIPODS Program)

  • Wooseok Ha

University of California Berkeley (Institute for Data Science)

  • Wooseok Ha

National Science Foundation (DMS-1654076)

  • Rina Foygel Barber

Office of Naval Research (N00014-20-1-2337)

  • Rina Foygel Barber

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Rena Schweizer for helping us download and process the gray wolf dataset used in the paper, Ben Peter for providing feedback and code for helping to construct the discrete global grids and preparing the human genetic dataset, and Hussein Al-Asadi, Peter Carbonetto, Dan Rice for helpful conversations about the optimization and modeling approach. We also acknowledge helpful feedback from Arjun Biddanda, Anna Di Rienzo, Matthew Stephens, the Stephens Lab, the Novembre Lab, and the University of Chicago 4th floor Cummings Life Science Center computational biology community. This study was supported in part by the National Science Foundation via fellowship DGE-1746045 and the National Institute of General Medical Sciences via training grant T32GM007197 to JHM and R01GM132383 to JN. WH was partially supported by the NSF via the TRIPODS program and by the Berkeley Institute for Data Science. RFB was supported by the National Science Foundation via grant DMS–1654076, and by the Office of Naval Research via grant N00014-20-1-2337.

Senior and Reviewing Editor

  1. George H Perry, Pennsylvania State University, United States

Reviewers

  1. Isabel Alves, University of Nantes, France
  2. Wesley Tansey

Publication history

  1. Received: August 8, 2020
  2. Accepted: June 7, 2021
  3. Version of Record published: July 30, 2021 (version 1)

Copyright

© 2021, Marcus et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 384
    Page views
  • 37
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Computational and Systems Biology
    2. Neuroscience
    Davide Spalla et al.
    Research Article

    Episodic memory has a dynamic nature: when we recall past episodes, we retrieve not only their content, but also their temporal structure. The phenomenon of replay, in the hippocampus of mammals, offers a remarkable example of this temporal dynamics. However, most quantitative models of memory treat memories as static configurations, neglecting the temporal unfolding of the retrieval process. Here, we introduce a continuous attractor network model with a memory-dependent asymmetric component in the synaptic connectivity, which spontaneously breaks the equilibrium of the memory configurations and produces dynamic retrieval. The detailed analysis of the model with analytical calculations and numerical simulations shows that it can robustly retrieve multiple dynamical memories, and that this feature is largely independent of the details of its implementation. By calculating the storage capacity, we show that the dynamic component does not impair memory capacity, and can even enhance it in certain regimes.

    1. Computational and Systems Biology
    2. Neuroscience
    Dhruva Raman, Timothy O'Leary
    Research Article

    Synaptic connections in many brain circuits fluctuate, exhibiting substantial turnover and remodelling over hours to days. Surprisingly, experiments show that most of this flux in connectivity persists in the absence of learning or known plasticity signals. How can neural circuits retain learned information despite a large proportion of ongoing and potentially disruptive synaptic changes? We address this question from first principles by analysing how much compensatory plasticity would be required to optimally counteract ongoing fluctuations, regardless of whether fluctuations are random or systematic. Remarkably, we find that the answer is largely independent of plasticity mechanisms and circuit architectures: compensatory plasticity should be at most equal in magnitude to fluctuations, and often less, in direct agreement with previously unexplained experimental observations. Moreover, our analysis shows that a high proportion of learning-independent synaptic change is consistent with plasticity mechanisms that accurately compute error gradients.