{'fraction_non_alphanumeric': 0.03430945957706031, 'fraction_numerical': 0.011095640066138716, 'mean_word_length': 5.125399786780384, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
KLFDAPC: a supervised machine learning approach for spatial genetic structure analysis

Abstract Geographic patterns of human genetic variation provide important insights into human evolution and disease. A commonly used tool to detect and describe them is principal component analysis (PCA) or the supervised linear discriminant analysis of principal components (DAPC). However, genetic features produced from both approaches could fail to correctly characterize population structure for complex scenarios involving admixture. In this study, we introduce Kernel Local Fisher Discriminant Analysis of Principal Components (KLFDAPC), a supervised non-linear approach for inferring individual geographic genetic structure that could rectify the limitations of these approaches by preserving the multimodal space of samples. We tested the power of KLFDAPC to infer population structure and to predict individual geographic origin using neural networks. Simulation results showed that KLFDAPC has higher discriminatory power than PCA and DAPC. The application of our method to empirical European and East Asian genome-wide genetic datasets indicated that the first two reduced features of KLFDAPC correctly recapitulated the geography of individuals and significantly improved the accuracy of predicting individual geographic origin when compared to PCA and DAPC. Therefore, KLFDAPC can be useful for geographic ancestry inference, design of genome scans and correction for spatial stratification in GWAS that link genes to adaptation or disease susceptibility.


Introduction
The genetic differentiation and substructure of human populations are impacted by spatially heterogeneous landscapes [1,2], social stratification [3,4], as well as culture [5]. For a long time, an interesting debate in population genetics is whether continuous clines or discrete clusters can better characterize human genetic variation [6][7][8][9]. However, without doubts, human population genetic structure exhibits a strong spatial pattern due to population history. On the global scale, this spatial pattern has been described by the isolation-by-distance model, where genetic differentiation between populations increases with increasing geographic distance, as a result of within-population genetic drift and reduced exchange of migrants between populations [10]. Recent studies, mostly at a continental scale (i.e. Europe, Asia), have shown that genetic variation significantly aligns with geography and exhibits spatial patterns that can be inferred by principal component analysis (PCA) [11][12][13] and model-based analyses [14,15]. Some studies have reported that the geographical spread of alleles favoured by natural selection contribute to local adaptation [16]. On the other hand, alleles underlying human complex diseases such as cancer, schizophrenia and heart disease also exhibit geographic patterns [17,18]. Therefore, successful detection of the genetic structure and correct inference of the individual geographic origin will be helpful for applications to personalized medicine, anthropology and forensics.
To date, several non-parametric approaches have been developed to make inferences about the genetic structure of populations and detect loci under selection. PCA is one of the most widely used approaches for these purposes [19]. The link with population structure was demonstrated by a series of studies [11,[20][21][22], which gradually showed that the proportion of the variance explained by the first principal component (PC) computed from the genome-wide single-nucleotide variation (SNP) genotype matrix is highly correlated with the fixation index, F ST . Furthermore, the linear superposition of PC maps has been used to infer the human geographic origin for various present-day and ancient individuals [11,18,23,24]. However, recent studies reported that under some complex scenarios, PCs are not sufficiently informative to represent population structure, as PCs are linear combinations of the variants without consideration of potential non-linear relationship [25,26]. It has also been suggested that it might be more robust to use non-linear functions of the top PCs, rather than more PCs, to capture non-linear spatial trends [27].
Despite its widespread use in population genetics, the spatial genetic structures represented by the PCs are, to some extent, not discernible between populations because PCs are a summary of the overall variance lumping together between-and within-population variation [28]. In contrast, Fisher linear discriminant analysis (LDA) [29] can maximize between-group variance while simultaneously minimizing within-group variance. In order to take advantage of this property and fulfil the assumption that variables submitted to LDA are perfectly uncorrelated, Jombart et al. proposed DAPC (Discriminant Analysis of Principal Components) [28], a hybrid statistical technique for dimensionality reduction that combines LDA and PCA. DAPC is statistically validated for linear inference and has been successfully applied to study population structure [28,30]. Nevertheless, individual scores in a population determined by LDA may be subject to bias as LDA assumes equal variance for all populations and weighs individuals in a population using the centroid of the genetic components of that population [31]. This property typically merges samples that might be from multiple populations into a single population. As it is the case with LDA, DAPC does not allow for within-group sub-structuring that may arise through migration or non-random mating [32].
Here, we propose a new method to overcome this limitation, Kernel Local Fisher Discriminant Analysis of Principal Components (KLFDAPC), which follows the same principle as DAPC but uses Kernel Local Fisher Discriminant Analysis (KLFDA) [32] instead of LDA. KLFDA is a more general approach for discriminant analysis that allows not only for within-group sub-structuring (multimodality) but also for non-linear associations among samples (individuals) within groups [32,33]. Therefore, our method combines non-linear and multimodal feature extraction of KLFDA and the dimension reduction of PCA, which helps overcome some of the limitations of PCA and DAPC.
We compared the performance of KLFDAPC for population structure inference and individual geography prediction with those of PCA and DAPC by applying all three methods to both simulated and empirical datasets.
The implementation of our method is freely available in the R package KLFDAPC at https://xinghuq.github.io/ KLFDAPC/.

Materials and methods
KLFDAPC is aimed at overcoming limitations of other popular dimensionality reduction techniques used to infer the genetic structure of populations. These include the presence of hidden genetic structure and non-linear genetic associations between samples. In principle, this could be achieved using KLFDA, an extension of Fisher discriminant analysis that preserves within-class local structure by evaluating the within-and between-class scatter in a local manner and incorporates non-linear associations using the kernel transformation technique [32,34]. However, KLFDA works well only for small datasets because the kernel transformation faces two key problems, (i) a heavy computational cost and (ii) large diagonals in the kernel matrix if the number of variables is greater than the number of samples (P > > n) [35], which is typically the case for genetic data comprising millions of genetic features (loci). As a solution to this problem, and following the example of DAPC [28], we propose to introduce an initial dimensionality reduction step that captures much of the variance present in the original data. This can be achieved using PCA. Thus, the method we propose integrates dimensionality reduction of PCA and the non-linear feature extraction of KLFDA, making the KLFDAPC scalable to genome-wide variation data.
KLFDAPC can be applied not only to genotype matrices but also to many other types of datasets, such as phenotypic traits and species counts. In what follows, we describe the steps required to implement it in general.

KLFDAPC formulation
The first step in the implementation of KLFDAPC is to obtain the PCs. The detailed description of this step is presented in Supplementary Methods available online at http://bib.oxfordjournals.org/. The next step is to conduct the KLFDA analysis using the first P PCs; exploration and guidelines to choosing the value of P is explained below (see section 'Tuning KLFDAPC parameters'). Here, we focus on describing in detail the formulation of KLFDA. As opposed to the PCA step, in KLFDA, we have to take into account the populations where the individuals were sampled in order to now consider the partitioning of genetic variation into its within-and between-population components. Thus, each individual has a population label y i ∈ (1, 2, . . . , c). Some individuals in a population could be recent migrants, and therefore, the population labels for these individuals might not represent their true source population. KLFDA uses a measure of local affinity (see below) that preserves the within population multimodality while maximizing the between population difference. Therefore, individuals in the same labelled group but actually from different populations can still be embedded and separated appropriately.
The main objective of KLFDA is to estimate the local Fisher transformation matrix, T LFDA , using the withinpopulation scatter matrix S (w) ∈ R n and the betweenpopulation scatter matrix S (b) ∈ R n , and then carry out a generalized eigenvalue decomposition. More details about the kernel local Fisher discriminant analysis formulation can be found in Refs. [32,33]. Below, we briefly recap the main steps for implementing KLFDA, which include (i) computing the kernel matrix M, and affinity matrix A; (ii) defining the local Fisher transformation matrix T LFDA in terms of within-and between-population scatter matrices S (w) and S (b) and (iii) computing an analytical form of the transformation matrix by solving a generalized eigenvalue problem.

Computing the kernel matrix
Once the data have been reduced to P PCs, KLFDA first transforms the PC scores into a kernel matrix M via non-linear mapping using a kernel function [36]. In this study, we use a Gaussian kernel, also known as radial basis function kernel. Let with σ determining the width of the Gaussian kernel [37]. This is a parameter that needs to be tuned, a step that is described in the section 'Tuning KLFDAPC parameters'. The kernel matrix can be viewed as a n-dimensional genetic distance matrix between pairs of individuals where each individual has a population label y i (1, 2, . . . , c). From the kernel distance matrix, we estimated the genetic affinities between individuals and then used them to calculate the within-and between-population weights.

Computing the affinity matrix
Here, the affinity matrix A i,j between individual i and individual j is computed using the k-nearest neighbour search with the local scaling method [38]. Let m i = (M ik ) and m j = (M jk ) be n-dimensional vectors of kernel distances between each one of these individuals and all other individuals calculated from Eq. 1. Let NN K i be the set of K-nearest neighbours of individual i under the Euclidean distance, where K is the neighbourhood size. If m i ∈ NN K j and m j ∈ NN K i , i and j are identified as neighbours; otherwise, they are non-neighbours. The elements of the affinity matrix A are given by [38] where σ i represents the local scaling of the data samples around m i , which is determined using where m K i is the vector of kernel distances for the Kth nearest neighbour of i. A i,j ∈ [0, 1] with small values indicating individuals have a low genetic affinity (i.e. are genetically far apart), and larger values indicating high affinity (genetically close individuals).

Calculating the local Fisher transformation matrix T LFDA
As described above, m i and m j are n-dimensional vectors of genetic distances for the i-th and j-th individuals, respectively. Let n y represent the sample size for population y so that n = c y=1 n y . Furthermore, let A As opposed to linear discriminant analysis (LDA) in which the within-group scatter and the between-group scatter are obtained using the group centroids and their overall average, here the scatter matrices in S (w) and S (b) are weighted by the affinities. In this case, genetically distant individuals within a population have less influence on S (w) and S (b) .
We define the local mixture scatter matrix as S (m) ≡ S (w) + S (b) ; thus, Therefore, i,j m i m j T (10) Equation (10) can be expressed in matrix form as where i,j . Using S (w) and S (b) , the local Fisher transformation matrix T LFDA can be defined as [33], T LFDA is the ratio of between-population (S (b) ) and within-population (S (w) ) variances, also known as the F-statistic in LDA, which is used to find the best transformation matrix to maximize Fisher's criterion [29].

Solution of the eigenvalue decomposition problem to obtain T LFDA
Noting that M is a symmetric matrix, T LFDA is obtained by solving [32,33], where In practice, Eq. (13) cannot be solved because L (w) is always singular. Therefore, Sugiyama [32] proposes regularizing ML (w) M and solving instead.
where I is the identity matrix and ε is a small constant used to regularize the within population distances to provide a more stable matrix.
Using neural networks to evaluate performance of dimensionality reduction methods and to assign individuals to their geographic origin The three methods we compared are aimed at selecting and combining the input variables into a reduced number of features that capture most of the genetic structure information present in the original dataset. Each method produces a distinct set of reduced features that cannot be directly compared. However, each set contains information about the geographic origin of individuals in the sample, and therefore, they can be used to assign individuals to source populations or geographic coordinates. Therefore, in order to compare the performance of the three methods, we implemented an artificial neuronal network that uses the top two or three features as predictors of each individual geographic or population origin and then calculates the accuracy of each method to assign individuals of known origin. To achieve this goal, we use synthetic data generated through a simulation study covering a wide range of population structure scenarios (see below). Based on the above-described rationale, we implemented a neural network ( Figure 1) and used it as a classifier to assign individuals to populations or as a regression to predict the individuals' geographic origin (latitude and longitude). Neural networks have been well explained in a series of studies [39][40][41][42][43][44][45]. A typical single hidden layer neural network consist of an input layer, a hidden layer and an output layer, with nodes in the hidden layers transforming the information between layers using a non-linear activation function (see Supplementary Information). The neural network is optimized based on a loss function that measures the fit of the predicted output to the true value. In the case of classification, we constructed a single-layer neural network with a logistic activation function to assign individuals to populations and used Shannon entropy as the loss function. In the case of regression problems, we used three hidden layers and a logistic activation function with the mean squared error as the loss function. More details about the neural networks and its tuning and fitting are given in Supplementary Methods available online at http://bib.oxfordjournals.org/.

Simulation study
To compare the performance of KLFDAPC over the existing commonly used approaches (PCA, DAPC), we simulated four scenarios that differ in spatial structure: island model, stepping stone model, hierarchical island model and hierarchical stepping stone model using the coalescent-based simulator fastsimcoal2 [46,47]. Under all scenarios, we simulated 44 independent chromosomes with 100 Kb DNA sequences per chromosome with a constant mutation rate of u = 1 × 10 −8 Figure 1. A neural network model for assigning individual membership and predicting the individual geographic coordinates. This framework is based on training a supervised neural network on the reduced genetic features from a dimensionality reduction technique (such as PCA, DAPC and KLFDAPC) given population labels or individual geographic coordinates. The reduced feature matrix (n × d, n is sample size and d is the number of reduced features) obtained from the genetic data are used as the predictor variables (A). If the population labels are provided (B), they are used as the response variable to carry out classification training through neural network (C). The individuals are assigned to the corresponding populations with an optimal neural network model. If the individual geographic coordinates are provided (B), the geographic coordinates are used as the response variable to carry out the regression training with neural network (C). An optimal neural network model is found and trained to predict the individual geographic coordinates. Finally, the accuracy of the reduced features for assigning individuals to correct populations or for predicting individual geographic coordinates is assessed (D) from the optimal neural network model. per bp per generation, which is typical of humans [48], virus [49], yeast and nematodes [50], and a recombination rate of r = 1 × 10 −8 per bp per generation, which is typical of mammals such as humans, but also of plasmids [51], bacteria [52] and human pathogens [53]. In the case of the non-hierarchical scenarios (island model and stepping-stone model), we assumed migration rates between populations to be 0.001, which leads to Nm = 1 allowing for the maintenance of polymorphism in each local population thanks to the inf lux of migrants (c.f. [54]). In the case of the hierarchical island and hierarchical stepping stone models, migration rate between pairs of populations within regions was also 0.001. However, migration between populations from different regions was set to 0.0001 (Nm = 0.1) to generate a strong hierarchical subdivision whereby migration from another region cannot, on its own, counteract the effect of genetic drift in any local population. The four scenarios and the respective parameter values used in the simulations are presented in Supplementary Table S1 available online at http://bib.oxfordjournals.org/.
We carried out 10 independent simulations for each scenario and sampled 200 individuals from each population. In total, we obtained 3200 individuals from 16 populations under each spatial scenario. Each scenario generated more than 27 000 polymorphic sites. We removed monomorphic SNPs and filtered the SNPs with a MAF > 0.05 and randomly selected 10 000 sites (biallelic) for downstream analysis.
For each scenario, we first carried out a PCA on the genotype matrix. We then found the number of PC axes and σ that maximized discriminatory power as explained below (see Supplementary Methods available online at http://bib.oxfordjournals.org/). DAPC was implemented using lda function from MASS package [55], which is initially employed by dapc function in adegenet package [56]. DAPC and KLFDAPC were conducted using the source population names as the group labels.

Tuning KLFDAPC parameters
Both DAPC and KLFDAPC require to find the optimal number of PCs as input. We chose the number of PCs that has the highest cumulative discriminatory power (accuracy) for population assignment (see above). More precisely, and as described above, we implemented a neural network classifier that in this particular step used PCs as predictive variables to assign individuals to populations under the various spatial scenarios we explored. Then, we evaluated the discriminatory power of the classifier as the number of PCs increased. As Supplementary Figure S2 available online at http:// bib.oxfordjournals.org/ shows, the cumulative discriminatory power reached an asymptote as the number of PCs increased from 5 to 20. We found that 20 PCs could correctly discriminate all 16 populations under all scenarios with an accuracy of 1.
Once the optimal number of PCs was found, we then applied a similar procedure to tune σ . In this case, we carried out KLFDAPC analyses of the simulated data based on the 20 top PCs and increasing values of σ (= 0.2, 05, 1, 2, 5). We then used the two and three first reduced features obtained from these analyses as predictive variables of the above described neural network classifier that assigned individuals to populations. When using only two reduced features (Supplementary Figure S3 available online at http://bib.oxfordjournals.org/), the effect of σ is not straightforward. Overall, the power decreases as σ increases but under the island model, it reaches a minimum at σ = 2 and then increases again, while under the hierarchical island model there is a steady and rapid decrease in power. On the other hand, under the stepping stone model, the decrease is very slow and the power remains above 90% under all σ values. The decrease in power is less clear under the hierarchical stepping stone model, but here again it remains above 90% under all σ values. Note, however, that when using three reduced features, there is a much clearer pattern of steady decrease in power as σ increases (Supplementary Figure S4 available online at http://bib.oxfordjournals.org/).
Based on this tuning process, we present most KLFDAPC results based on the top 20 PCs and σ = 0.5 but we also present results for varying σ to visually explore how this parameter affects the clustering of samples in geographic and reduced feature (kernel-induced) space.

Comparing methods performance using simulated data
As mentioned before, we implemented a neural network to assign individuals to source populations based on the top two or three reduced features obtained from each dimensionality reduction method. We then estimated the neural network accuracy and Cohen's Kappa coefficient (κ) to assign individuals to the labelled populations using the information provided by the reduced features (Supplementary Methods available online at http://bib.oxfordjournals.org/).
The use of DAPC and KLFDAPC assumes that population labels assigned to individuals correspond to discrete demographic units. However, in practice, individuals sharing the same habitat patch may represent genetic mixtures, which can introduce errors in the assignment of individuals to populations using the neural network. To compare the performance of DAPC and KLFDAPC under this particularly difficult scenario, we generated synthetic data under a hierarchical island model and then created genetically mixed regions consisting of populations from two different regions (see Figure 4A). The DAPC and KLFDAPC analyses were carried out using the first 20 PCs, and in the case of KLFDAPC, a sigma value of 5.

Application to real datasets
We tested the performance of all three approaches to predict the geographic locations of individuals using two datasets, the European populations from POPRES datasets [57] (dbGaP accession number phs000145.v4.p2) and the Han Chinese populations from CONVERGE data [58] (http://www.ebi.ac.uk/ena/data/view/PRJNA289433). The details on data quality control can be found in Supplementary Methods available online at http://bib.oxford journals.org/. To assess the performance of the three approaches (PCA, DAPC and KLFDAPC) in inferring the geographic origin, the predictive performance (R 2 observed values versus the predicted values) between different methods was assessed using model resampling with neural network regression. We also used standard methods to compare the predictive power of PCA, DAPC and KLFDAPC: (i) correlation analysis and (ii) Procrustes analysis. Details of testing using each metric can be found in Supplementary Methods available online at http://bib.oxfordjournals.org/.

Discriminatory power
We assessed the discriminatory power of KLFDAPC for population delineation using the simulated scenarios (Methods and Supplementary Figure S1 available online at http://bib.oxfordjournals.org/), including two non-hierarchical spatial models (the classic island model, Supplementary Figure S1A available online at http://bib.oxfordjournals.org/, and stepping stone model, Supplementary Figure S1C available online at http://bib.oxfordjournals.org/) and two hierarchical spatial models (the hierarchical island model, Supplementary Figure S1B available online at http://bib.oxfordjournals.org/, and hierarchical stepping stone model, Supplementary Figure S1D available online at http://bib.oxfordjournals.org/). We first carried out PCA on the genotype matrix of 3200 sampled individuals under all simulated scenarios. Based on the tuning step described above, we retained the first 20 PCs to conduct the DAPC and KLFDAPC analyses. In addition, we set the KLFDAPC parameter σ to 0.5 (but also present some results using other values). In a first step, we present 2D and 3D plots as the representation of population genetic structure (Figure 2 and Supplementary Figure S5 available online at http://bib.oxfordjournals.org/). We then tested the predictive power of the three approaches to assign individuals to the sampled populations using neural networks ( Figure 3).
All methods successfully discriminated between regions under the hierarchical spatial scenarios (hierarchical island model, Figure 2B, F and J, and hierarchical stepping stone model, Figure 2D, H and L). However, PCA and DAPC both failed to clearly delineate local populations under all four scenarios (Figure 2A-H). In contrast, KLFDAPC clearly distinguished genetic stratification among local populations under the stepping stone model based on the first two reduced features ( Figure 2K and L) and under the hierarchical stepping stone model based on the first three reduced features (Supplementary Figure S5 available online at http://bib.oxfordjournals.org/). In this latter case, the first reduced feature discriminated between the two higher level regions in the hierarchy while the second and third reduced features together discriminated among populations within regions. The second reduced feature discriminated among populations within region 1 while the third feature discriminated among populations within region 2. Overall, KLFDAPC performed better in identifying the populations under the isolation-bydistance models (the stepping-stone and hierarchical stepping-stone model) than under the island models To quantitatively compare the performance of each method (PCA, DAPC and KLFDAPC) in describing genetic structuring, we implemented an artificial neural network that used the first three reduced features obtained from each method to assign individuals to populations (see Methods and Supplementary Information). We then assessed the classification accuracy of each set of reduced features in terms of classification accuracy and Cohen's Kappa coefficient [59]. Consistent with the graphical representation of the spatial structures, the discriminatory accuracy and Cohen's Kappa coefficient (κ) for KLFDAPC were much higher than those achieved by PCA and DAPC under all scenarios (Figure 3). Note that the accuracy and κ of all methods with three reduced features improved over those obtained when only two axes were used (Supplementary Figure S6 available online at http://bib.oxfordjournals.org/) but the strongest improvement was observed for KLFDAPC. We therefore recommend considering more features (i.e. the first three features in Supplementary Figure S5 available online at http://bib.oxfordjournals.org/) to fully characterize population structure under complex scenarios.
As Supplementary Figure S7 illustrates, the pattern of local genetic aggregation is sensitive to the parameter σ . KLFDAPC introduces non-linear genotypic associations using a Gaussian kernel in which σ controls the strength of dispersal/aggregation of the local structure. Lowering the σ values of KLFDAPC could increase the discrimination of discrete clusters, thus increasing the ability of KLFDAPC to delineate distinct aggregates (Supplementary Figs S4 and S7 available online at http://bib.oxfordjournals.org/). In summary, KLFDAPC outperformed PCA and DAPC in discriminating population genetic structure. Furthermore, in the case of stepping stone models, KLFDAPC was able to characterize a spectrum of genetic structure from continuous genetic gradients to discrete clusters with appropriate kernel parameter values (i.e. σ in Gaussian kernel). Therefore, we recommend users to vary kernel parameter values to explore how it inf luences the results.

The effects of hidden substructure on the delineation of regions under hierarchically structured scenarios
An important problem when using supervised learning is the effect of group mislabelling due to hidden substructure. This can happen, for example when sampling takes place in wintering or feeding areas that can receive migrants from several different regions. To investigate this issue, we considered a scenario where four breeding grounds contributed each to two different feeding grounds (see Figure 4A). Thus, each feeding ground consisted of genetic mixtures of two distinct genetic clusters. Figure 4B shows that DAPC was unable to group individuals according to the region where they bred. On the other hand, KLFDAPC correctly grouped together individuals from the same breeding region. This difference is due to the different way in which the two methods calculate the within-class scatter matrix (i.e. the distance between the position of each sample in multidimensional space and the average position of the class). More precisely, DAPC simply uses the class centroid while KLFDAPC takes into account the genetic affinity between samples.

Analysis of POPRES data
We tested the performance of PCA, DAPC and KLFDAPC to predict the geographic locations of European individuals using POPRES data. The first two PCs of the POPRES data only accounted for 0.29% and 0.15% of the total SNP variance, respectively. These two PCs were remarkably aligned with the map of Europe ( Figure 5A and D; Table 1, PC1 versus longitude, R = 0.872, PC2 versus latitude, R = 0.873), as previously reported [11]. DAPC and KLFDAPC analyses using the top 20 PCs also provided accurate geographic representations of the genetic samples ( Figure 5B, C, E and F). Moreover, DAPC and KLFDAPC improved the alignment of genetic samples with their locations. For example, DAPC and KLFDAPC rectified the projected geographic locations between Turkey (TR) and Albania (AL) by bringing them close to each other, while also locating PL (Poland) to its correct position and IR and UK samples closer to their correct location ( Figure 5; Supplementary Table S2   available online at http://bib.oxfordjournals.org/). KLF-DAPC uses a Gaussian kernel to take into account non-linear associations between samples, where σ determines the decay in the association [32]. Increasing σ from 1 to 5 induced a gradual change of the projected locations between Cyprus (CY) and south-east European countries, as well as Russia (RU) and north-east European countries (Supplementary Figure S9 available online at http://bib.oxfordjournals.org/). KLFDAPC performed better at predicting the geographic locations of POPRES individuals than PCA and DAPC (Tables 1 and 2; Supplementary Table S2 available online at http://bib.oxfordjournals.org/). Overall, lower σ values tend to aggregate individuals into compact clusters, while high σ values make the individuals more scattered.
When predicting the individual origin via a neural network model using the first two reduced features as the predictor variables, KLFDAPC performed the best among the three methods to predict the individual longitude and the latitude (Tables 1 and 2). DAPC and PCA showed similar power in predicting the individual geographic locations (Tables 1 and 2). Overall, KLFDAPC showed

Analysis of CONVERGE data
We also assessed the efficacy of our method to infer the individual geographic origin for a large Han Chinese population from the CONVERGE dataset. The p values marked with bold indicate statistically significant differences between two methods at the p < 0.05 level. Upper diagonal: estimates of the difference. Lower diagonal: P-value for H0: difference = 0. R 2 (true value versus predicted value) was estimated after 5-fold cross-validations, repeated five times. Comparison between models was tested using a pairwise t-test between 100 resamples. P-value adjustment: Bonferroni.
CONVERGE data consist of individuals from 24 out of 33 administrative divisions (19 provinces, 4 municipalities and 1 autonomous region) across China. We used the individual-level birthplace information at the province level to denote the geographic origin of each sample [13]. PCA performed poorly in recapitulating the geography of individuals, as PC1 corresponded poorly to both latitude and longitude (Supplementary Table S3 available online at http://bib.oxfordjournals.org/, Figure 6A and D; PC1 versus longitude R = 0.0508, PC1 versus latitude R = −0.351). However, we can still observe a significant North-South gradient (Supplementary Table S3 available online at http://bib.oxfordjournals.org/, Figure 6A  Due to the limitation of sample collection (samples are mainly collected from eastern, coastal China provinces), there was a poor correlation between the first reduced genetic feature and longitude. We note that PC did worse at predicting longitude while DAPC and KLFDAPC improved predictive accuracy but still performed poorly at predicting longitude (Table 3). All methods failed to accurately predict individual longitude using neural networks (Table 3). However, they still performed well in predicting the individuals' latitude ( Table 3). The predictive power analysed using neural networks showed that PCA did worst in predicting the individual latitude among these three approaches (R 2 = 0.727, Table 3). Even though DAPC (R 2 = 0.740) did better than PCA (P = 5.97E−07), KLFDAPC (σ = 0.5, R 2 = 0.767) performed significantly better than both PCA (P = 2.2E−16) and DAPC (P = 2.2E−16) (Tables 3 and 4). In addition, the predictive power R 2 obtained from neural networks is higher than the conventional linear correlation coefficient (Table 3; Supplementary Table S3 available online at http://bib.oxfordjournals.org/).
In summary, KLFDAPC with a σ value of 0.5 outperforms PCA and DAPC in all aspects when predicting the individual geographic origins of Han Chinese people using the CONVERGE dataset (Tables 3  and 4; Supplementary Table S3 available online at http://bib.oxfordjournals.org/), suggesting that genetic features produced by KLFDAPC seem to be a better surrogate for geographic coordinates than PCs.

Discussion
The availability of large genomic databases has pushed researchers to put aside model-based methods in favour of non-parametric approaches, such as PCA [22,60] and DAPC [28]. In this study, we introduced KLFDAPC, a non-linear approach for inferring population genetic structure and individual geographic origin. Using a neural network with KLFDAPC reduced features as predictive variables, we tested the performance of  KLFDAPC for inference of individual population membership and geographic origin. We showed that KLFDAPC outperformed both PCA and DAPC in population structure discrimination and in predicting individual geographic origin using simulated scenarios and empirical population datasets (Figs 2 and 3). Analyses of the POPRES dataset showed that all three methods retrieved a strong correspondence between genetic structure and geography, but KLFDAPC outperformed both other methods in inferring the individual geographic Comparison between models was tested using a pairwise t-test between 100 resamples. P-value adjustment: Bonferroni. All models failed to predict the longitude because of the small variance between individuals and poor correlation between features and longitude.
locations (Tables 1 and 2; Supplementary Table S2 available online at http://bib.oxfordjournals.org/). When applying the three methods to the Han Chinese population, PCA exhibited poor performance in characterizing the individual geography. Both DAPC and KLFDAPC provided a much better alignment between genetic structure and geography and remarkably improved the predictive accuracy of individual geographic origin compared to PCA ( Figure 6). Again, KLFADPC outperformed both PCA and DAPC in predicting the individual geography in the CONVERGE dataset (Tables 3 and 4). Overall, our study highlights the remarkable performance of KLFDAPC in identifying population genetic structure and in predicting individual geographic origin. We thus propose that KLFDAPC, which extracts the non-linear genetic features and also allows for within-population hidden structuring, may be a useful alternative to PCA and DAPC for many population genetics studies.
There are model-based alternatives to the machine learning methods we considered in our study. However, they generally require users to define a non-linear function modelling the slope of allele frequency across geographic locations [14,15]. Any pre-specified parametric function is unlikely to sufficiently capture complex geographic patterns in genetic variation, such as multiple modes or peaks in the allele frequency surface as reported by Yang et al. [14]. In contrast, KLFDAPC is a non-parametric method that can incorporate both non-linear genetic associations between individuals and hidden sub-structuring within populations. It does not require a non-linear function being defined but only needs to choose an appropriate kernel.
One particular feature of KLFDAPC that differentiates it from PCA and DAPC is the need to specify the value of the kernel parameter, for example σ in the Gaussian kernel used in our study. Our simulation study suggests that values between 0.2 and 5 provide optimal results under most scenarios. However, we consider that it is useful to vary this parameter to explore genetic structuring at different spatial scales. Our simulation study showed that low values of σ help to clearly delineate different regions in hierarchically structured scenarios, while larger values tend to highlight within-region structuring. Also, real-data analyses showed that increasing σ tends to produce more scattered or continuous patterns (Supplementary Figs S9 and S10 available online at http://bib.oxfordjournals.org/). Thus, by tuning σ , it is possible to maximize the power to predict individual geographic origin and, therefore, to better account for population stratification and spatial effects in GWAS analyses.
We found that, in analyses of the CONVERGE data, PCA tended to provide a poor indication of the correct geographic origin of individuals. A previous study of the population structure of Han Chinese people based on the CONVERGE dataset filtered ∼30% of individuals on the basis of poorly imputed genotype in order to reveal strong correlation between the first two PCs and longitude and latitude [13]. In our study, we were able to demonstrate much stronger correlation between reduced genetic features and geography using KLFDAPC compared to PCA, while retaining all 10 461 individuals in the analysis. This result suggests that KLFDAPC is more robust to varying data quality issues such as missing or poorly imputed genotypes than conventional PCA. In particular, missing genotypes is a common problem in genomic studies. The k-nearest neighbour algorithm is often used for imputing missing genotypes [61][62][63]. Note that KLFDAPC uses the k-nearest neighbour algorithm (KNN) to calculate the local genetic affinity, which could strongly decrease the influence of genotyping errors on inference. Therefore, KLFDAPC overcomes this artefact experienced by PCA by preserving the non-linear genetic features and local genetic affinity.
A common problem with supervised or populationbased methods such as LDA and F ST is the requirement for pre-grouping. Pre-grouping might be subjective or arbitrary, as we rarely know if some individuals in a group might have immigrated recently from other groups or some individuals within a group have unknown origin, which introduce bias to the inference [64] and could mask important evolutionary processes such as migration and cross-breeding [65]. Our results showed that DAPC suffers from this problem as it minimizes the within-group variation based on the group means and can thus lead to erroneous assignment of individual to populations ( Figure 4). Thus, extensive genetic mixing or admixing can represent a great challenge for DAPC and lead to inaccurate representations of genetic structure.
Unlike DAPC and other within-group minimizing/averaging approaches, KLFDAPC is less affected by pregrouping ( Figure 4). The main reason is that KLFDAPC computes the within-and between-population affinity based on the KNN weight, minimizing the influence of 'local outliers' (i.e. migrants) on the inference of population structure. Therefore, KLFDAPC can preserve the multimodal genetic structure within populations and overcomes the problems caused by the group-mean (or group-centroid)-based method.
Existing studies propose to use neural network approaches for population assignment [66,67] and prediction of geographic origin [68] using genetic data. Guinand et al. [66] show that assignment tests based on a neural network classifier with one hidden layer generally outperforms the other methods. In Supplementary Methods available online at http://bib.oxfordjournals.org/, we show that a KLFDAPC-based single-layer neural network largely outperforms the method of Guinand et al. [66]. The most likely explanation is that when using individuals' genotypes as input, a single hidden layer may not be capable of extracting all the information present in the raw data. Using instead the reduced features of KLFDAPC overcomes this limitation. Another way of improving the accuracy of individual assignment based on genotype data is proposed by Battey et al. [68], who used a neural network with several hidden layers to predict the individual locations based on allele counts and known locations.
Besides predicting the individual geographic origin, we also implemented a genome scan approach to identify genomic regions involved in local adaptation based on KLFDAPC (https://xinghuq.github.io/KLFDAPC/ articles/Genome_scan_KLFDAPC.html). Many ecological, evolutionary and medical datasets are complex and may exhibit hidden genetic structuring. We expect that KLFDAPC would be very helpful in these particular situations. First, population structure integrating multimodal structure within large populations could correct spurious results in inferring individual geographic origin (c.f. Han Chinese example in our study). The genetic structure with a non-linear outlier identification model, such as neural network (as opposed to linear regression models in pcadapt [69] and LFMM [70]), would provide a better characterization of genomic regions under natural selection or involved in adaptation. Therefore, KLFDAPC could improve the power of detecting loci under selection or involved in local adaptation. Also, KLFDAPC could also be used to correct for stratification in genome-wide association studies, which so far have done so using PCA.
Machine learning algorithms, from simple general linear regression [71], PCA [22,60], to random forest [72], extreme gradient boosting [73], as well as neural networks [74], have enabled us to capture the systematic signatures of biological or genetic patterns from genomic samples, allowing for the association of genes to phenotypes/diseases and facilitating molecular-based medical applications [75][76][77]. KLFDAPC represents a new addition to the population genomics toolbox but it is also potentially applicable to other Omics data throughout the biological sciences, including applications in medicine and agriculture.

Key Points
• Most species' genetic diversity exhibits geographic patterns and several methods have been proposed to characterize it. Through a simulation study and real-data analyses, the limitations of the commonly used methods, PCA and DAPC, for spatial genetic data analysis were revealed. • A supervised machine learning method, KLFDAPC, is introduced to rectify the limitations of PCA and DAPC by capturing non-linear information and preserving the multimodal space of samples. • KLFDAPC outperformed PCA and DAPC in discriminatory power and in predicting the geographic origin of individuals. • KLFDAPC can be useful for geographic ancestry inference, design of genome scans and correction for spatial stratification in GWAS. • KLFDAPC is freely available at https://xinghuq.github.io/ KLFDAPC/index.html.

Supplementary Data
Supplementary data are available online at https:// academic.oup.com/bib.

Author contribution
X.Q. and O.E.G. designed the study. C.K.W.C. reviewed the design, acquired the datasets and provided computational resources. X.Q. compiled the package, carried out the analyses and interpreted results with input from C.W.K.C. and O.E.G. X.Q. and O.E.G. wrote the manuscript with the input from C.K.W.C. All authors contributed to editing and revisions of the manuscript.

Data and code availability
The input files and scripts for generating simulation are available at https://github.com/xinghuq/KLFDAPC/ tree/sm/Simulation_files. The POPRES data with request approval #90291-1 is accessible via dbGaP Study accession number phs000145.v4.p2. The CONVERGE data are available at the European Nucleotide Archive (ENA)