Deep Learning for Population Genetic Inference — TEI-tagged view

Sara SheehanpersName, Yun S. SongpersName (2016date). “Deep Learning for Population Genetic Inference.” PLOS Computational Biologytitle·j.

Object of study — a deep-learning framework for jointly inferring demographyterm·key and selectionterm·dom from genomic data

Research question — how to infer complex population geneticterm·key models whose likelihood cannot be computed

Contribution — evidence that deep learningterm·mth works for population geneticterm·key inference and finds informative features

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Deep Learningterm·mth for Population Geneticterm·key Inference Abstract Given genomic variation data from multiple individuals, computing the likelihood of complex population geneticterm·key models is often infeasible. To circumvent this problem, we introduce a novel likelihood-freeterm·key inference framework by applying deep learningterm·mth, a powerful modern technique in machine learning. Deep learningterm·mth makes use of multilayer neural networks to learn a feature-based function from the input (e.g., hundreds of correlated summary statisticsterm·key of data) to the output (e.g., population geneticterm·key parameters of interest). We demonstrate that deep learningterm·mth can be effectively employed for population geneticterm·key inference and learning informative features of data. As a concrete application, we focus on the challenging problem of jointly inferring natural selectionterm·key and demographyterm·key (in the form of a population size change history). Our method is able to separate the global nature of demographyterm·key from the local nature of selectionterm·dom, without sequential steps for these two factors. Studying demographyterm·key and selectionterm·dom jointly is motivated by Drosophila , where pervasive selectionterm·dom confounds demographic analysis. We apply our method to 197 African Drosophila melanogasterterm·key genomes from Zambia to infer both their overall demographyterm·key, and regions of their genome under selectionterm·dom. We find many regions of the genome that have experienced hard sweepsterm·key, and fewer under selectionterm·dom on standing variation (soft sweepterm·key) or balancing selectionterm·key. Interestingly, we find that soft sweepsterm·key and balancing selectionterm·key occur more frequently closer to the centromere of each chromosome. In addition, our demographic inference suggests that previously estimated bottlenecks for African Drosophila melanogasterterm·key are too extreme.

Introduction With the advent of large-scale whole-genome variation data, population geneticists are currently interested in considering increasingly more complex models. However, statistical inference in this setting is a challenging task, as computing the likelihood of a complex population geneticterm·key model is a difficult problem both theoretically and computationally.

In this paper, we introduce a novel likelihood-freeterm·key inference framework for population genomics by applying deep learningterm·mth , which is an active area of machine learning research. To our knowledge, deep learningterm·mth has not been employed in population genomics before. A recent survey article [ 1 ] provides an accessible introduction to deep learningterm·mth, and we provide a high-level description below. Our general goal in this paper is to demonstrate the potential of deep learningterm·mth as a powerful framework for population geneticterm·key analysis that can allow accurate inference of previously intractable models.

As a concrete example, we consider models of non-equilibrium demographyterm·key and natural selectionterm·key, for which multi-locus full-likelihood computation is prohibitive, and apply deep learningterm·mth to the challenging problem of jointly inferring demographic history and selectionterm·dom (see [ 2 ] for a recent review of this topic). One reason why this problem is difficult is that demographyterm·key (for example, a bottleneck in population size) and selectionterm·dom can leave similar signals in the genome. Untangling the two factors directly has rarely been attempted; most methods that estimate selectionterm·dom try to demonstrate robustness to demographic scenarios, rather than estimating demographic parameters jointly with selectionterm·dom. The demographic models considered in this paper are restricted to a single population with time-varying effective population sizeterm·var, but the overall framework presented here can be applied to more general demographic models. We anticipate many algorithmic, modeling, and application directions that could emerge from this proof of concept.

Our focus on the joint inference of demographyterm·key and selectionterm·dom is motivated by Drosophila , where previous demographic estimates [ 3 ] may have been confounded by pervasive selectionterm·dom. The reverse has occurred as well, with selectionterm·dom estimates being confounded by demographyterm·key [ 4 ]. See [ 5 ] for a more thorough discussion of the role selectionterm·dom plays in the Drosophila genome. To test our method, we simulate data under a variety of realistic demographies for an African population of Drosophila melanogasterterm·key . For each demographic history, we simulate many regions under different selectionterm·dom parameters. We then apply our tailored deep learningterm·mth method using hundreds of potentially informative summary statisticsterm·key. From the output of this trained network, we demonstrate that parameters can be accurately inferred for test datasets, and interpret which statistics are making the biggest contributions. After training our deep network, we also apply it to analyze 197 Drosophila melanogasterterm·key genomes from Zambia, Africa [ 6 ] to learn about their effective population sizeterm·var change history and selective landscape.

A summary of related works Several machine learning methods have been developed for selectionterm·dom only, often focusing on classifying the genome into neutral versus selected regions. Examples include methods based on support vector machinesterm·mth (SVMsterm·mth) [ 7 – 9 ] or boosting [ 10 – 12 ]. Often these methods demonstrate robustness to different demographic scenarios, but do not explicitly infer demographyterm·key.

Many methods have been developed to infer ancestral population size changes, including PSMC [ 13 ], diCal [ 14 , 15 ], and MSMC [ 16 ]. These methods make the underlying assumption that selectionterm·dom does not significantly bias the results. This is perhaps somewhat true for humans, but is not a valid assumption for many species such as Drosophila , for which selectionterm·dom seems ubiquitous throughout the genome.

Few previous works have addressed both population size changes and selectionterm·dom. Galtier et al. [ 17 ] developed a likelihood-based method for distinguishing between a bottleneck and selectionterm·dom. They applied their method to Drosophila data to conclude that a series of selective sweeps was more likely than a bottleneck, but did not explicitly infer parameters for both selectionterm·dom and demographyterm·key. In their method, Galtier et al. assumed that demographic eventsterm·dom affect the entire genome, whereas selectionterm·dom is a local phenomenon. In contrast, Gossmann et al. [ 18 ] estimated the effective population sizeterm·var locally along the genome, and reported that it is correlated with the density of selected sites. To make our results as easily interpretable as evolutionary events as possible, we estimate global effective population sizeterm·var changes.

Approximate Bayesian Computationterm·mth (ABCterm·mth) [ 19 , 20 ] is a likelihood-free inference methodterm·mth based on simulating datasets and comparing their summary statisticsterm·key. (A more detailed description of the framework is provided below.) This approach has been used to study various complex population geneticterm·key models (e.g., [ 21 – 23 ]) for which likelihood computation is prohibitive. Partly due to several influential theoretical works [ 24 – 31 ], the popularity of ABCterm·mth has grown rapidly over the past decade. ABCterm·mth’s main advantages are that it is easy to use and is able to output a posterior distribution. There are a few challenges, however: 1) ABCterm·mth suffers from the “curse of dimensionality,” with decreasing accuracy and stability as the number of summary statisticsterm·key grows [ 32 ]. 2) ABCterm·mth typically uses a rejection algorithm, so the simulated datasets are not used optimally, and often a very large number are required. 3) It is difficult to interpret the output of ABCterm·mth in terms of which statistics were the most informative for each parameter. 4) Although ABCterm·mth can be used very effectively for continuous parameter inference, it would be difficult to use ABCterm·mth for a complex scenario that involved joint inference of continuous parameters and categorical distributions. To our knowledge, ABCterm·mth has thus not been previously applied to the problem of jointly inferring demographyterm·key (using continuous effective population sizesterm·var) and selectionterm·dom (using discrete classification). The work of Bazin et al. [ 33 ] estimates demographic and selectionterm·dom parameters using a hierarchical ABCterm·mth approach, though not in an explicit joint framework.

We see our deep learningterm·mth method as a complimentary approach to, rather than a direct replacement of, other likelihood-free inference methodsterm·mth such as ABCterm·mth or SVMterm·mth. Deep learningterm·mth has several advantages such as making full use of datasets, elegantly handling correlations between summary statisticsterm·key, and producing interpretable features of the data. It also produces an “inverse” function from the input summary statisticsterm·key to the population geneticterm·key parameters of interest. However, in contrast to ABCterm·mth, it does not provide a posterior distribution.

A brief introduction to deep learningterm·mth Deep learningterm·mth has its beginnings in neural networks, which were originally inspired by the way neurons are connected in the brain [ 34 ]. Neural networks have been studied for over 60 years and a huge body of literature exists on the topic. A neural network is typically used to learn a complex function between inputs (data) and outputs ( response variables ) in the absence of a model, so it can be thought of as a generalized regression framework. While standard regression and classification methods involve fitting linear combinations of fixed basis functions, a neural network tries to learn basis functions (usually non-linear) appropriate for the data. A neural network architecture consists of multiple layers of computational units (nodes), with connections between the layers but not between nodes within a layer. Within a layer, each node computes a transformation (usually non-linear) of the outputs from the previous layer. S1 Fig illustrates a simple feed-forward neural network with a single hidden layer. The phrase “deep learningterm·mth” refers to algorithms for learning deep neural network architectures with many hidden layers.

The universal approximation theorem [ 35 , 36 ] states that any continuous function on compact subsets of R n can be uniformly approximated by a feed-forward neural network with a single hidden layer, provided that the number of nodes in the hidden layer is sufficiently large and the transformation (called the activation function) associated with each node satisfies some mild conditions. However, it can be challenging to learn the weights of such a network and to interpret the hidden layer. So as learning problems became more complex, it was desirable to train networks with more hidden layers. Since their introduction over 30 years ago, deep architectures have proved adept at modeling multiple levels of abstraction. However, they were notoriously difficult to train since their objective functions are non-convex and highly non-linear, and the level of non-linearity increases with the number of layers in the network.

A major breakthrough was made in 2006 when Hinton and Salakhutdinov [ 37 ] showed that a deep feed-forward neural network can be trained effectively by first performing unsupervised “pretraining” one layer at a time, followed by supervised fine-tuning using a gradient-descent algorithm called back-propagation [ 38 ]. (Simply put, pretraining provides a good initialization point for non-convex optimization.) They applied their learning algorithm to dimensionality reduction of images and achieved substantially better results than PCA-based methods.

Following the work of Hinton and Salakhutdinov, deep learningterm·mth has been applied to various challenging problems in computer science over the last several years, making groundbreaking progress. Deep learningterm·mth broke long-standing records for accuracy that had been set by approaches based on hand-coded rules. Well-known examples include automatic speech recognition (transforming spoken words into typed text) [ 39 , 40 ] and computer vision (automatically classifying images into different categories and tagging objects/individuals in photos) [ 41 ].

Many variations have been developed, including dropout , which attempts to learn better and more robust features of the data (see [ 42 , 43 ]). Deep learningterm·mth has also been applied to problems in neuroscience [ 44 ] and computational biology [ 45 – 48 ], but has not been used for population geneticsterm·key before. We demonstrate here that deep learningterm·mth can accommodate more complex models than existing techniques, and can provide a complimentary approach to existing likelihood-freeterm·key inference in population geneticsterm·key. For example, deep learningterm·mth could be used to select optimal statistics for ABCterm·mth (as demonstrated in [ 49 ]), or another method.

ABCterm·mth background Rejection-based methods have been used since the late 1990’s (see [ 19 , 50 ]) to estimate population geneticterm·key parameters when the likelihood is difficult to compute. Early improvements to ABCterm·mth quickly helped make it a popular method for a variety of scenarios (see [ 20 ] for a good introduction). ABCterm·mth works by simulating many datasets under a prior for the parameters of interest. Then these datasets are reduced to a vector of summary statisticsterm·key that are ideally informative for the parameters. The summary statisticsterm·key that are closest to the summary statisticsterm·key for the target dataset are retained, and the corresponding parameters used to estimate the desired posterior distributions. ABCterm·mth has many advantages, including ease of use and the output of a posterior distribution. However, there are a few shortcomings, which various works have sought to address.

In its simplest form, ABCterm·mth does not naturally handle correlated or weakly informative summary statisticsterm·key, which can add noise to the data. To tackle this problem, many methods for dimensionality reduction or wisely selecting summary statisticsterm·key have been proposed (see [ 26 , 29 , 31 , 51 , 52 ], and [ 53 ] for a dimensionality reduction comparison). However, simple reductions cannot always learn subtle relationships between the data and the parameters. Expert pruning of statistics helps some methods, but given the lack of sufficient statistics, valuable information can be eliminated. Blum and François [ 28 ] suggested performing a dimensionality reduction step on the summary statisticsterm·key via a neural network similar to S1 Fig . This reduction is similar in spirit to the work presented here, although there are many algorithmic and application differences.

A second issue with ABCterm·mth is the rejection step, which does not make optimal use of the datasets that are not retained. The more statistics and parameters used, the more datasets must be simulated and rejected to properly explore the space, making the interaction between these two issues problematic. To address the problem of rejecting datasets, different weighting approaches have been proposed (see [ 28 ] for a good example of how the estimation error changes as fewer datasets are rejected). The idea is to keep more datasets, but then weight each retained dataset by its distance to the target dataset. However, few approaches utilize all the datasets in this way, and the most popular implementation of ABCterm·mth ( ABCtoolbox , [ 54 ]) typically still rejects most of the simulated datasets by default.

One final issue with ABCterm·mth is the black-box nature of the output. Given the distances between the simulated datasets and the target dataset, and the posterior, there is no clear way to tell which statistics were the most informative.

Results We first describe the simulated dataset we generated to investigate joint inference of demographic history and selectionterm·dom. In the Methods section, we describe how to transform each region (as well as each region of the real data) into a set of 345 summary statisticsterm·key, and also detail our deep learningterm·mth algorithm. In what follows, we present the results of our method in a variety of different contexts.

Simulating data To create a simulated dataset that is appropriate for our scenario of interest, we first define the response variables we would like to estimate. For simplicity, we consider piecewise-constant effective population sizeterm·var histories with three epochs: a recent effective population sizeterm·var N 1 , a bottleneck size N 2 , and an ancient size N 3 . Further, we define a genomic region as belonging to 4 different selectionterm·dom classes: no selectionterm·dom (neutral), positive directional selectionterm·dom (hard sweepterm·key), selectionterm·dom on standing variation (soft sweepterm·key), and balancing selectionterm·key. See [ 55 ] for a more complete analysis of the different types of selectionterm·dom in Drosophila . To accurately reflect that demographyterm·key affects the entire genome while selectionterm·dom affects a particular region, we simulate many genomic regions under the same demographic history, but the selectionterm·dom class for each region is chosen independently. See Fig 1 for a simplified illustration of the data.

Input data for the demographyterm·key and selectionterm·dom scenario. For each demographic history (bottleneck), we simulated many different genomic regions. Each region can either have no selectionterm·dom (neutral), one site with an de novo mutation under positive selectionterm·dom (hard sweepterm·key), one site under balancing selectionterm·key, or one standing variant under positive selectionterm·dom (soft sweepterm·key).

To simulate data, we used the program msms [ 56 ]. To make our simulated data as close to the real data as possible, we simulated n = 100 haplotypes, and correspondingly downsampled the Zambia Drosophila melanogasterterm·key dataset to match. We repeated the following procedure 2500 times. First we selected three population sizes for the demographic model, then simulated 160 regions with these sizes, 40 for each selectionterm·dom scenario. Each region was 100 kb, with the selected site (if present) occurring randomly in the middle 20 kb of the region. We used a baseline effective population sizeterm·var N ref = 100,000, a per-base, per-generation mutation rate μ = 8.4 × 10 −9 [ 57 ], and a per-base, per-generation recombination rate r equal to μ , as inferred by PSMC [ 13 ]. This ratio of recombination to mutation rate is somewhat lower than what has been inferred in other works (see [ 58 , 59 ]), so we also simulated a testing dataset with a higher recombination rate, described below. However, it is likely that there is significant variation in this ratio across the genome [ 60 ]. We used a generation time of 10 generations per year. Based on the population size change results of PSMC, we chose the time change-points for demographic history to be t 1 = 0.5 and t 2 = 5 in coalescent units. Scaled effective population sizeterm·var parameters λ i : = N i / N ref and their prior distributions are below: Recent effective population sizeterm·var scaling factor: λ 1 ∼ Unif(3, 14).

Bottleneck effective population sizeterm·var scaling factor: λ 2 ∼ Unif(0.5, 6).

Ancient effective population sizeterm·var scaling factor: λ 3 ∼ Unif(2, 10).

For the selectionterm·dom classes, the different types are shown below: Class 0, neutral : no selectionterm·dom, neutral region.

Class 1, hard sweepterm·key : positive selectionterm·dom on a de novo mutation (i.e., hard sweepterm·key). For the selectionterm·dom coefficient, we used s ∈ {0.01, 0.02, 0.05, 0.1}, with 10 regions for each value of s . The onset of selectionterm·dom was chosen to be 0.005 in coalescent units, which provided sweeps at a variety of different stages of completion at the present time. We discarded datasets where the frequency of the selected allele was 0 (early loss of the beneficial mutation due to drift), but for a large fraction of datasets, the beneficial allele had not yet fixed.

Class 2, soft sweepterm·key : positive selectionterm·dom on standing variation (i.e., soft sweepterm·key). The selectionterm·dom coefficients and selectionterm·dom start time were chosen as in the hard sweepterm·key scenario, but now the initial frequency of the beneficial allele was chosen to be 0.001.

Class 3, balancing : heterozygote advantage, balancing selectionterm·key. The selectionterm·dom start time and selectionterm·dom coefficients were chosen in the same fashion as Class 1.

Using this strategy, we simulated 2500 different demographic histories, with 160 regions for each one, for a total of 400,000 datasets. To build 160 regions for each demographyterm·key, we simulated 40 datasets for each of the selectionterm·dom classes above.

Results for demographyterm·key and selectionterm·dom on simulated data Of the 400,000 total datasets, we used 75% for training and left out 25% for testing. In Table 1 , we show the effective population sizeterm·var results for a network with 3 hidden layers respectively with 25, 25, and 10 nodes (5 layers total including input and output). The best results were found when we averaged the statistics for all of the datasets (160 “regions”) with the same demographyterm·key, then ran this average through the trained network.

Deep learningterm·mth results for population sizes under the scenario with demographyterm·key and selectionterm·dom, using a network with 3 hidden layers of sizes 25, 25, and 10. We evaluate the results in three ways, using the relative error of the estimates: | N est − N true |/ N true . First, for each test demographyterm·key, we average the statistics for each dataset, and then run these values through the training network. Second, for each test demographyterm·key, we run the datasets through the network one by one, then average the predictions. Finally, we perform the second procedure, but with only the regions we classified as neutral. We note that the most ancient size ( N 3 ) is always the least accurately estimated.

Table 2 shows the results of our three different size prediction methods (see Methods ) for an example demographyterm·key. For each parameter, we include a range between the 2.5th and 97.5th quantiles. To further investigate uncertainty in our estimates, we also include a violin plot of the three population sizes, shown in Fig 2A . In this example (and more generally), our point estimate of the most ancient size N 3 is the least accurate, but the most recent size N 1 has the largest uncertainty. When considering a single dataset, there is not always a clear winner among our three prediction methods. Overall, the average statistic method is usually the most accurate, followed by the neutral regions method.

Results for an example demographyterm·key, with the true population sizes shown in the first row. The rows below each dotted line are the three different size predictions for this dataset. Average stat : average the statistics for each region within the same demographyterm·key, then run these averages through the trained network. Final : run the statistics for each region through the trained network, then average the results. Neutral regions : only average the results for the regions we predict as neutral.

Violin plots of population size estimates for (A) an example simulated dataset and (B) the real dataset. In both cases, we have the largest uncertainty for the most recent population size. For the simulated dataset, the true population sizes are shown in blue, and we note that our point estimate of the most ancient population size is the least accurate.

To analyze the selectionterm·dom results, we calculated a confusion matrix in Table 3 to show which datasets of each class were classified correctly, or misclassified as belonging to a different class. The rows of a confusion matrix represent the true class of each dataset, and the columns represent the predicted or called class. An ideal confusion matrix would have all 1’s (100%) on the diagonal, and 0’s off the diagonal. Our most frequent error is classifying hard sweepterm·key datasets as neutral (row 2, column 1 of the confusion matrix). We hypothesized that this was either because selectionterm·dom occurred anciently and quickly, or because selectionterm·dom occurred recently and the sweep was not yet complete. To test this, we examined the results conditional on the frequency of the selected allele at present. The results shown in Table 4 suggest that in fact many of the sweeps were not complete, which is why the regions sometimes appeared neutral. Regardless, this type of false negative error is arguably preferable to falsely calling many regions to be under selectionterm·dom when they are in fact neutral.

Confusion matrix for the selectionterm·dom predictions, in the demographyterm·key and selectionterm·dom scenario. Each row represents the datasets that truly belong to each selectionterm·dom class. Each column represents the datasets that were actually classified as each selectionterm·dom class. Ideally we would like all 1’s down the diagonal, and 0’s in the off-diagonal entries. The largest number in each row is shown in boldface. We can see that neutral datasets are the easiest to classify, and sometimes regions under selectionterm·dom (hard sweepsterm·key in particular) look neutral as well (first column). The overall percentage of misclassified datasets was 6.2%.

Hard sweepterm·key results, broken down by the present-day frequency ( f ) of the selected allele (which we would not know for a real dataset). We defined “low” frequency as f ∈ [0, 0.3), “moderate” frequency as f ∈ [0.3, 0.7), and “high” frequency as f ∈ [0.7, 1]. We can see that if the frequency of the selected allele is low, the region is often misclassified as neutral, since the selective sweep is not yet complete. However, if the frequency is moderate or high, the dataset is usually classified correctly as a hard sweepterm·key.

We also wanted to test the impact of unsupervised pretraining using autoencoders (see Methods ). Tables 5 and 6 compare the results for a randomly initialized network to a network initialized using autoencoders. These results demonstrate that pretraining is very effective. Due to the non-convex nature of the optimization problem, random initialization is most likely finding a poor local minima for the cost function.

Relative error on the test dataset, for a deep network with 6 hidden layers. For the results in the first row, the weights of the entire network were initialized randomly, then optimized. In the second row, the weights were initialized using autoencoders for each layer. The positive impact of unsupervised pretraining is clear; random initialization causes the optimization to get stuck in a local minima.

Confusion matrix for the selectionterm·dom predictions, compared between random initialization (top) and autoencoder initialization (bottom), for a deep network with 6 hidden layers. Again, ideally we would like all 1’s down the diagonal, and 0’s in the off-diagonal entries. The largest number in each row is shown in boldface. When the network is initialized randomly, almost every dataset is classified as neutral; the network has not really learned anything meaningful from the input data. The overall percentage of misclassification is 74.8% for random initialization, while it is only 6.1% for autoencoder initialization.

Lastly, we wanted to assess the impact of training the deep network on data that was created under different conditions than the testing data. In particular, we simulated testing data where the ratio of the recombination rate to the mutation rate was equal to 4 (instead of 1 as in our previous results). The results of this experiment are shown in S1 Table . This rather large difference in the testing data negatively impact the results, but much more so for selectionterm·dom than for the effective population sizesterm·var. This is perhaps due to the dependence of selectionterm·dom inference on LDterm·mth and haplotype statistics, which would be more affected by recombination rate variation.

We also simulated testing data with bottlenecks that were more extreme than the training data, to assess the impact of a misspecified range of demographic models. Specifically, the most recent population size was outside the range of the training data for all the testing datasets. The results are shown in S2 Table . In contrast to the previous testing scenario, here the population size results are somewhat negatively impacted (largely for the most recent size), but the selectionterm·dom results are only impacted slightly. In general, identifying model misspecification for continuous population geneticterm·key parameters can be done as follows: if the predicted value of a normalized continuous parameter is outside the range [0,1], then the testing/real data parameter is outside the range of the training/simulated data parameter. This simple criterion can help users identify when their simulated data needs to be modified to provide accurate results for their real data.

Results on real Drosophila melanogasterterm·key data We then ran the real Drosophila melanogasterterm·key data through our trained network just like any other test dataset. We considered only the chromosome arms 2L, 2R, 3L, and 3R in this study. We partitioned each chromosome arm into 20 kb windows and ran our method on five consecutive windows at a time, sliding by 20 kb after each run. Classification was performed on the middle 20 kb window. The population size predictions for some of these regions (roughly 14%) were outside the range of the simulated data, so we discarded them from all further analysis. These regions tended to be near the telomeres and centromeres, where assembly is often poor.

For the demographic history, the population size results are shown in Table 7 . The first row (average statistic method) is our best estimate, which we also plot and compare with other histories in Fig 3 . Our history is close to the PSMC result, although more resolution would be needed for a proper comparison. We based our time change points ( t 1 and t 2 ) on PSMC, so it would be surprising if these histories were notably inconsistent. The expansion after the bottleneck is roughly consistent with previous results citing the range expansion of Drosophila melanogasterterm·key (out of sub-Saharan Africa) as beginning around 15,000 years ago [ 61 ]. This range expansion likely led to an effective population sizeterm·var increase like the one we infer. Our expansion time is also consistent with bottleneck estimates in non-African Drosophila melanogasterterm·key populations, such as the Netherlands population investigated by Thornton and Andolfatto [ 58 ]. In Fig 2B , we include a violin plot of the population size estimates from each region. We have the largest uncertainty for the most recent population size.

Population size results for African Drosophila from Zambia, rounded to the nearest hundred. The first row (predictions based on averaging the statistics for each region in the Drosophila genome) represents our best estimate of the population sizes.

Comparison between demographic histories. On the left in green is PSMC [ 13 ], run on the entire genome for a subset of the data ( n = 20). In blue is our history from the first line of Table 7 . On the right in red is the history from Duchen et al. [ 3 ] (note the change in the y -axis scale). The difference between our estimates and these previous estimates is likely due in large part to their assumption of a very short (1000 generations) bottleneck. Given this timing constraint, only a very severe bottleneck could fit the data. A more gradual bottleneck seems more realistic, although we do not have a simple explanation for why there was a bottleneck at all.

The number of windows classified as Neutral, Hard Sweepterm·key, Soft Sweepterm·key, and Balancing Selectionterm·key are 1101, 2455, 179, 404, respectively. See S3 Table for further details. If we restrict our analysis to regions classified with probability greater than 0.9999, then we find 42 hard sweepsterm·key, 34 soft sweepsterm·key, and 17 regions under balancing selectionterm·key. In S4 Table , we include a table of these high-confidence windows, along with which genes are found in each window. We also include a plot, Fig 4 , of where the selected regions fall throughout the genome (restricted to regions with a probability at least 0.95). Interestingly, soft sweepsterm·key and balancing selectionterm·key seem to occur more frequently closer to the centromere of each chromosome. There is also an excess of hard sweepsterm·key on chromosome arm 2L.

Plot of the selected regions we infer for African Drosophila melanogasterterm·key . Each of the 4 subplots is a chromosome arm (2L, 2R, 3L and 3R). We restricted the plot to selected sites with a probability of at least 0.95.

Upon examining the genes in the top regions classified to be under selectionterm·dom, we find several notable results. In the hard sweepterm·key category, a few of the top regions harbor many genes involved in chromatin assembly or disassembly. Further, there are two regions each containing exactly one gene, where the gene has a known function. These are Fic domain-containing protein (detection of light and visual behavior), and charybde (negative regulation of growth).

In the soft sweepterm·key category, we find many genes related to transcription, and several related to pheromone detection (chemosensory proteins). We also find four regions each containing exactly one gene, with that gene having a known function. These are gooseberry (segment polarity determination), steppke (positive regulation of growth), Krüppel (proper body formation), and Accessory gland protein 36DE (sperm storage, [ 62 ]). Fly embryos with a mutant Krüppel (German for “cripple”) gene have a “gap” phenotype with the middle body segments missing [ 63 ]. We also identify three interesting genes in regions with other genes: seminal fluid proteins 36F and 38D (see [ 64 ]), and deadlock , which is involved in transposable element silencing, which preserves genome integrity [ 65 ].

Finally in the balancing selectionterm·key category, one interesting gene we find is nervana 3 , which is involved in the sensory perception of sound. We also find balancing selectionterm·key regions containing the genes cycle (circadian rhythms) and Dihydropterin deaminase (eye pigment), although there are other genes within these putatively selected regions.

Most informative statistics One of the advantages of deep learningterm·mth is that it outputs an optimal weight on each edge connecting two nodes of the network, starting with the input statistics and ending with the output variables. By investigating this network, we are able to determine which statistics are the “best” or most helpful for estimating our parameters of interest. Here we use two different methods, one using permutation testing and the other using a perturbation approach described in S1 Algorithm .

To perform permutation testing, we randomly permuted the values of each statistic across all the test datasets and then measured the resulting decrease in accuracy. For each output (3 population sizes and selectionterm·dom category), we then retained the 25 statistics that caused the largest accuracy decreases. The results are shown in the 4-way Venn diagram in Fig 5 , which highlights statistics common to each subset of outputs. For each statistic name in the diagram, close , mid , and far represent the genomic subregion where the statistic was calculated. The numbers after each colon refer to the position of the statistic within its distribution or order; for the folded site frequency spectrum (SFS) statistics, it is the number of minor alleles. Several statistics are informative for all the parameters, including the raw number of segregating sites (“S”) in each subregion and the number of singletons (“SFS, far: 1”). There is generally more overlap between the statistics for selectionterm·dom and N 1 or N 2 , than with the statistics for selectionterm·dom and N 3 , possibly because selectionterm·dom occurs on a more recent timescale. Tajima’s D is also most important for selectionterm·dom (which it was designed to detect) and N 1 . Identity-by-state (IBS) statistics seem to play an increasingly more important role as the population sizes become more ancient. Perhaps unsurprisingly, LDterm·mth statistics are most important for selectionterm·dom. Interestingly, the SFS is not as informative as one might have anticipated. However, the number of singletons is very informative, especially for selectionterm·dom and N 1 , which could have been foreseen given that recent events shape the terminal branches of a genealogy, and thus the singletons. The distances between segregating sites (“BET” statistics) do not generally seem very helpful. Neither does H2, the frequency of the second most common haplotype. It is comparing H1 to H2 (via the H12 statistic) that is helpful.

A Venn diagram of most informative statistics for each output variable ( N 1 , N 2 , N 3 , and selectionterm·dom). For each variable, the top 25 statistics were chosen using permutation testing. The Venn diagram captures statistics common to each subset of output variables, with notable less informative statistics shown in the lower right. Close, mid, and far represent the genomic region where the statistic was calculated. The numbers after each colon refer to the position of the statistic within its distribution or order. For the SFS statistics, it is number of minor alleles. For each region, there are 50 SFS statistics, 16 BET statistics (distribution between segregating sites), 30 IBS statistics, and 16 LDterm·mth statistics.

The results ( S2 Fig ) from the perturbation method share some similarities and have some differences from the Venn diagram in Fig 5 . An interesting set of four statistics is informative for all response variables, including very long IBS tracts close to the selected site (“IBS, close: 30”) and the H12 statistic. Additionally, “LDterm·mth, mid: 4” represents low LDterm·mth between sites close to the selected site and sites mid-range from the selected site (likewise for LDterm·mth “LDterm·mth, far: 4”). Many low LDterm·mth pairs could signal a lack of selectionterm·dom, and vice versa. IBS statistics are generally not as helpful for selectionterm·dom, but extremely informative for the population sizes, especially N 2 and N 3 , as also observed in the permutation method. Bottlenecks can have a significant impact on IBS tract distributions, so it makes sense that N 2 (the bottleneck size) is the most reliant on IBS statistics. The number of segregating sites S is generally quite informative, especially for selectionterm·dom.

It is interesting that the list of notable less-informative statistics is very similar for both permutation and perturbation methods. Overall, roughly the same set of informative statistics were identified by each method, but the precise position of each statistic on the Venn diagram shifted.

Comparison with ABCtoolbox Although ABCterm·mth is not well suited for our scenario of interest and deep learningterm·mth is a complementary method, we wanted to find a scenario where we could compare the performance of these two methods. To this effect, we restricted the analysis to estimating (continuous) demographic parameters only. We used the popular ABCtoolbox [ 54 ], using the same training and testing datasets as for deep learningterm·mth. For ABCterm·mth, the training data represents the data simulated under the prior distributions (uniform in our case), and each test dataset was compared with the training data separately. We retained 5% of the training datasets, and used half of these retained datasets for posterior density estimation. Overall, we used 75% of the datasets for training and 25% for testing.

We tested two scenarios, one with the full set of summary statisticsterm·key (345 total), and the other with a reduced set of summary statisticsterm·key (100 total). For the reduced set of summary statisticsterm·key, we chose statistics which seemed to be informative: the number of segregating sites, Tajima’s D , the first 15 entries of site frequency spectrum, H1, and the distribution of distances between segregating sites. It is interesting that the informative statistics identified above do not match this set of “intuitively” chosen statistics. The results are shown in Table 8 , which suggest that deep learningterm·mth produces more accurate estimates of the recent population size ( N 1 ) than does ABCtoolbox, whereas they have comparable accuracies for more ancient sizes.

A comparison between ABCtoolbox and deep learningterm·mth for demographyterm·key only. Out of 1000 demographies (160,000 datasets total), 75% were used for training and 25% for testing. In this scenario, deep learningterm·mth generally outperforms ABCtoolbox, as measured by the relative error: | N est − N true |/ N true . There is generally more improvement using deep learningterm·mth when the number of statistics is larger.

Regularization of the network weights In deep learningterm·mth, one hyper-parameter that should be investigated closely is the regularization parameter λ in the cost function, also called the weight decay parameter. (See Methods for details). If λ is set to be too high, large weights will be penalized too much, and interesting features of the data cannot be learned well. But if λ is set too low, the weights tend to display runaway behavior. Due to this balance, a validation procedure is typically used to find the right λ. In our case, the additional runtime of simulating more data and performing more training would be too computationally expensive, but we do provide a small validation study in S3 Fig Runtime In terms of runtime, the majority is spent simulating the data. During training, most of the runtime is spent fine-tuning the deep network, which requires computing the cost function and derivatives many times. To speed up this computation, our deep learningterm·mth implementation is parallelized across datasets, since each dataset adds to the cost function independently. This significantly improves the training time for deep learningterm·mth, which can be run in a few days on our simulated dataset (400,000 regions) with a modest number of hidden nodes/layers. Once the training is completed, an arbitrary number of datasets can be tested more or less instantaneously. In contrast, each of the “training” datasets for ABCterm·mth must be examined for each test dataset. This takes several weeks for a dataset of this size (which is why we tested ABCterm·mth on a subset of the data), although it could be parallelized across the test datasets. We include a runtime comparison on the reduced dataset in Table 9 .

Approximate runtime results. Simulating the data (which is the same for both methods) and computing the summary statisticsterm·key is the most time consuming part of the process, although it is highly parallelized. ABCtoolbox could also be parallelized, but we did not implement that here. For the “filtered statistics” row, the 345 statistics have been downsampled to 100 for all the datasets.

Discussion In this paper, we have sought to demonstrate the potential of deep learningterm·mth as a new inference framework for population genomic analysis. This investigation is still in its infancy, and there are many directions of future work. One important advantage of deep learningterm·mth is that it provides a way to distinguish informative summary statisticsterm·key from uninformative ones. Here we have presented two methods for extracting informative statistics given a trained deep network. Other methods are possible, and it is an open question which one would produce the most accurate results. It would be interesting to down-sample statistics using a variety of methods, then compare the results. Overall, learning more about how various summary statisticsterm·key relate to parameters would be useful for population geneticsterm·key going forward.

The prospect of using deep learningterm·mth to classify regions as neutral or selected is very appealing for subsequent demographic inference. There are other machine learning methods that perform such classification, but they are generally limited to two classes (selected or neutral). One exception is a study in humans [ 66 ], which classifies genomic regions as neutral or under positive, negative, or balancing selectionterm·key. Although their approach does not jointly infer selectionterm·dom and demographyterm·key, it would be interesting to see their method used on Drosophila .

We infer many hard sweepsterm·key in African Drosophila , which is perhaps expected given their large effective population sizeterm·var. However, when restricting our analysis to selected regions with high confidence, the numbers of hard sweepsterm·key and soft sweepsterm·key are comparable. It would be interesting to analyze our results in the context of a simulation study by Schrider et al. [ 67 ], which found that regions classified as soft sweepsterm·key are often truly the “shoulders” of hard sweepsterm·key. This possibility is worth investigating, as the signatures of soft sweepsterm·key and soft shoulders are extremely similar. In our simulations, both hard and soft sweepsterm·key are relatively recent, which makes them easier to detect. It would be useful to incorporate datasets with a wider range of selectionterm·dom onset times, which could capture regions of the real data that are currently being classified as neutral. Right now we have 4 different selectionterm·dom classes which are all forms of positive selectionterm·dom. To incorporate background or purifying selectionterm·dom, a 5th class of (weakly) deleterious mutations could be added.

One aspect of the simulations that warrants future study is the effect of the ratio of recombination to mutation rate, ρ / θ . From running PSMC on 20 haplotypes from Zambia, we inferred ρ / θ = 1, whereas other studies have found ρ / θ = 6.8 [ 59 ] and ρ / θ = 10 [ 58 ] for a population from Zimbabwe. Both these works also show significant variation in this ratio across different Drosophila populations. Andolfatto and Przeworski [ 60 ] also find significant variation in this ratio across the genome. As noted by Haddrill et al. [ 59 ], recent bottlenecks can make this ratio appear lower than it actually is, which is possibly the case for our dataset. As shown in the Supporting Information, misspecification of this ratio can lead to poor selectionterm·dom classification results. In general, assuming a lower recombination rate than reality leads to conservative results, since sites are assumed to be less independent [ 58 – 60 ]. Our takeaway from this analysis is that ρ / θ should be accurately estimated using a more robust method than PSMC or incorporated as a model parameter, with datasets simulated under a variety of ratios. It would also be interesting to investigate variation in the recombination rate across the genome, which could be incorporated into future simulations.

Deep learningterm·mth can make efficient use of even a limited number of simulated datasets. In this vein, it would be interesting to use an approach such as ABCterm·mth MCMC [ 24 ] to simulate data, and then use deep learningterm·mth on these simulated datasets. Alternatively, deep learningterm·mth could be used to select informative statistics for a subsequent method such as ABCterm·mth [ 49 ]. Blum and François [ 28 ] used a single-layer neural network to perform (non-linear) regression to obtain a posterior mean, which they used to adjust the parameters of the posterior sample for a more accurate posterior distribution. It is unclear whether deep networks (where the theory is much more complicated) could be used in the same way, although there is the concept of a residual which could be used for parameter adjustment in an rejection-based framework. Such hybrid approaches could be a fruitful area of further exploration, potentially extending the generality and usefulness of ABCterm·mth.

We would also like to apply deep learningterm·mth to a wider variety of scenarios in populations genetics. Population structure and splits would be examples, although these scenarios would most likely require very different sets of summary statisticsterm·key.

On the computer science side, deep learningterm·mth has almost exclusively been used for classification, not continuous parameter inference. It would be interesting to see the type of continuous parameter inference presented here used in other fields and applications.

Finally, machine learning methods have been criticized for their “black-box” nature. In some sense they throw away a lot of the coalescent modeling that we know to be realistic, although this is included somewhat in the expert summary statisticsterm·key of the data. It would be advantageous to somehow combine the strengths of coalescent theory and the strengths of machine learning to create a robust method for population geneticterm·key inference.

Methods Deep learningterm·mth details In this section we provide the theory behind training deep networks. The notation in this section follows that of [ 68 ]. Let x ( i ) be the vector of summary statisticsterm·key for dataset i , and y ( i ) the vector of response variables that dataset i was simulated under. If we have m such datasets, then together {( x (1) , y (1) ), …, ( x ( m ) , y ( m ) )} form the training data that will be used to learn the function from summary statisticsterm·key to response variables. Deep neural networks are a way to express this type of complex, non-linear function. The first layer of the network is the input data, the next layers are the “hidden layers” of the network, and the final layer represents the network’s prediction of the response variable.

Cost function for a deep network An example deep network is shown in Fig 6 . The collection of weights between layer ℓ and layer ℓ +1 is denoted W ( ℓ ) = ( w j k ( ℓ ) ) , where w j k ( ℓ ) is the weight associated with the connection between node j in layer ℓ and node k in layer ℓ +1. The biases for layer ℓ is denoted b ( ℓ ) . The total number of layers (including the input and output layers) is denoted L , and the number of hidden nodes in layer ℓ is denoted u ℓ . The main goal is to learn the weights that best describe the function between the inputs and the outputs.

An example of a deep neural network with two hidden layers. The first layer is the input data (each dataset has 5 statistics), and the last layer predicts the 2 response variables. The last node in each input layer (+1) represents the bias term. Here the number of layers L = 4, and the number of nodes (computational units) in each layer is u 1 = 5, u 2 = 3, u 3 = 3, and u 4 = 2 (these counts exclude the biases).

To learn this function, we first describe how the values of the hidden nodes are computed, given a trial weight vector. The value of hidden node j in layer ℓ ≥ 2 is denoted a j ( ℓ ) , and is defined by a j ( ℓ ) = f z j ( ℓ ) and z j ( ℓ ) = W j ( ℓ - 1 ) · a ( ℓ - 1 ) + b j ( ℓ - 1 ) , (1) where W j ( ℓ - 1 ) is the j th column of the weight matrix W ( ℓ - 1 ) (i.e., all the weights going into node j of layer ℓ − 1), a ( ℓ - 1 ) = ( a k ( ℓ - 1 ) ) is a vector of the values of all the nodes in layer ℓ − 1 (with a ( 1 ) = x , the input vector), and f ( z ) : = 1 1 + exp ( - z ) (2) is the activation function . Here we use a logistic function, but other functions can be used. Another common activation function is the hyperbolic tangent function.

Hence, given the input data and a set of weights, we can feed forward to learn the values of all hidden nodes, and a prediction of the response variables. These predictions are usually denoted by h W , b ( x ( i ) ) for our hypothesis for dataset i , based on all the weights W = ( W (1) , …, W ( L -1) ) and biases b = ( b (1) , …, b ( L -1) ). We will discuss different ways to compute the hypothesis function later on. To find the best weights, we define a loss function based on the square of the l 2 -norm between this hypothesis and the true response variables. This loss function is given by J ( W , b ) = 1 m ∑ i = 1 m 1 2 ‖ h W , b ( x ( i ) ) - y ( i ) ‖ 2 . (3) The goal of deep learningterm·mth is to find the weights ( W and b ) that minimize this loss function. To efficiently find these optimal weights, we can use backpropagation to find the gradient. The intuition behind this approach is that once we have found the hypothesis, we then want to see how much each of the weights contributed to any differences between the hypothesis and the truth. Therefore we start at the last hidden layer, see how much each of those weights contributed, then work our way backwards, using the gradient of the previous layer to compute the gradient of the next layer. For this procedure we need to compute the partial derivatives of the cost function with respect to each weight.

Consider a single training example x with associated output y . The cost for this dataset is a term in Eq (3) and is denoted by J ( W , b , x , y ) : = 1 2 ‖ h W , b ( x ) - y ‖ 2 . Then, define δ j ( ℓ ) : = ∂ J ( W , b , x , y ) ∂ z j ( ℓ ) , where z j ( ℓ ) , defined in Eq (1) , is the input to the activation function for node j in layer ℓ . We first consider δ j ( L ) for the final layer. Noting that the j th entry of h W , b ( x ) is f ( z j ( L ) ) , we get δ j ( L ) = [ f ( z j ( L ) ) - y j ] f ′ z j ( L ) . Based on this initialization, we can recursively compute all the δ variables: δ j ( ℓ ) = ∑ k = 1 u ℓ + 1 w j k ( ℓ ) δ k ( ℓ + 1 ) f ′ z j ( ℓ ) . Now we can use the δ variables to recursively compute the partial derivatives for one dataset: ∂ J ( W , b , x , y ) ∂ w j k ( ℓ ) = a j ( ℓ ) δ k ( ℓ + 1 ) , ∂ J ( W , b , x , y ) ∂ b k ( ℓ ) = δ k ( ℓ + 1 ) . Finally, putting all the datasets together we get ∂ J ( W , b ) ∂ w j k ( ℓ ) = ∑ i = 1 m ∂ J ( W , b , x ( i ) , y ( i ) ) ∂ w j k ( ℓ ) and ∂ J ( W , b ) ∂ b k ( ℓ ) = ∑ i = 1 m ∂ J ( W , b , x ( i ) , y ( i ) ) ∂ b k ( ℓ ) . Since we can compute the derivatives using this backpropagation algorithm, we can use the LBFGS optimization routine (as implemented in [ 69 ]) to find the weights that minimize the cost function.

Unsupervised pretraining using autoencoders It is possible to train a deep network by attempting to minimize the cost function described above directly, but in practice, this proved difficult due to the high-dimensionality and non-convexity of the optimization problem. Initializing the weights randomly before training resulted in poor local minima. Hinton and Salakhutdinov [ 37 ] sought to initialize the weights in a more informed way, using an unsupervised pretraining routine. Unsupervised training ignores the output (often called the “labels”) and attempts to learn as much as possible about the structure of the data on its own. PCA is an example of unsupervised learning. In Hinton and Salakhutdinov, the unsupervised pretraining step uses an autoencoder to try to learn the best function from the data to itself, after it has gone through a dimensionality reduction step (this can be thought of as trying to compress the data, then reconstruct it with minimal loss). Autoencoding provides a way to initialize the weights of a deep network that will ideally be close to optimal for the supervised learning step as well. See S4 Fig for a diagram of an autoencoder.

Training an autoencoder is an optimization procedure in itself. As before, let W (1) be the vector of weights connecting the input x to the hidden layer a , and x (2) be the vector of weights connecting a to the output layer x ^ , which in this case should be as close as possible to the original input data. We again typically use the logistic function shown in Eq (2) as our activation function f , so we can compute the output using: a j = f ( W j ( 1 ) · x + b j ( 1 ) ) and x ^ k = f ( W k ( 2 ) · a + b k ( 2 ) ) . If a linear activation function is used instead of a logistic function, the hidden layer becomes the principle components of the data. This makes dimensionality reduction with an autoencoder similar in spirit to PCA, which has been used frequently in genetic analysis (see [ 70 ] for an example). However, the non-linear nature of an autoencoder has been shown to reconstruct complex data more accurately than PCA. Using backpropagation as we did before, we can minimize the following autoencoder cost function using all m input datasets: A ( W , b ) = 1 m ∑ i = 1 m 1 2 ‖ x ^ ( i ) - x ( i ) ‖ 2 . The resulting weights W * (1) and b * (1) will then be used to initialize the weights between the first and second layers of our deep network. The weights W * (2) > and b * (2) are discarded. To initialize the rest of the weights, we can repeat the autoencoder procedure, but this time we will use the hidden layer a * = ( a j * ) , where a j * = f ( W j * ( 1 ) · x + b j * ( 1 ) ) , as our input data, and feed it through the next hidden layer. In this way we can use “stacked” autoencoders to initialize all the weights of the deep network. Finally, the supervised training procedure described in the previous section can be used to fine-tune the weights to obtain the best function from the inputs to the response variables.

When the number of hidden nodes is large, we would like to constrain an autoencoder such that only a fraction of the hidden nodes are “firing” at any given time. This corresponds to the idea that only a subset of the neurons in our brains are firing at once, depending on the input stimulus. To create a similar phenomenon for an autoencoder, we create a sparsity constraint that ensures the activation of most of the nodes is close to 0, and the activation of a small fraction, ρ , of nodes is close to 1. Let ρ ^ j be the average activation of the hidden node j : ρ ^ j = 1 m ∑ i = 1 m a j ( x ( i ) ) , where a j ( x ( i ) ) is the value of the j th hidden node when activated with dataset x ( i ) . To ensure sparsity, we would like ρ ^ j to be close to ρ , our desired fraction of active nodes. This can be accomplished by minimizing the KL divergence: ∑ j = 1 u 2 KL ( ρ ‖ ρ ^ j ) = ∑ j = 1 u 2 ρ log ρ ρ ^ j + ( 1 - ρ ) log 1 - ρ 1 - ρ ^ j , where u 2 is the number of nodes in the hidden layer of the autoencoder. We multiply this term by a sparsity weight β .

In addition, a regularization term is included, which prevents the magnitude of the weights from becoming too large. To accomplish this, we add a penalty to the cost function that is the sum of the squares of all weights (excluding the biases), weighted by a well-chosen constant λ, which is often called the weight decay parameter . Including both sparsity and regularization, our final autoencoder cost becomes: A λ ( W , b ) = 1 m ∑ i = 1 m 1 2 ‖ x ^ ( i ) - x ( i ) ‖ 2 + β ∑ j = 1 u 2 KL ( ρ ‖ ρ ^ j ) + λ 2 ∑ ℓ = 1 2 ∑ j = 1 u ℓ ∑ k = 1 u ℓ + 1 ( w j k ( ℓ ) 2 .

We also regularize the weights on the last layer during fine-tuning, so our deep learningterm·mth cost function becomes: J λ ( W , b ) = 1 m ∑ i = 1 m 1 2 ‖ h W , b ( x ( i ) ) - y ( i ) ‖ 2 + λ 2 ∑ j = 1 u L - 1 ∑ k = 1 u L w j k ( L - 1 ) 2 .

The final layer: Continuous variable vs. classification In population geneticsterm·key, often we want to estimate continuous response variables. To compute our hypothesis for a response variable, based on a set of weights, we could use a logistic activation function, Eq (2) , as we did for the other layers. However, the logistic function is more suitable for binary classification. Instead, we use a linear activation function, so in the case of a single response variable, our hypothesis for dataset i becomes h W , b linear ( x ( i ) ) = [ a ( L - 1 ) ( x ( i ) ) ] T W ( L - 1 ) + b ( L - 1 ) . In other words, it is the dot product of the activations of the final hidden layer and the weights that connect the final hidden layer to the response variables.

For such classification results, if we had two classes, we could use logistic regression to find the probability a dataset was assigned to each class. With more than two classes, we can extend this concept and use softmax regression to assign a probability to each class. If we have K classes labeled 1, …, K , we can define our hypothesis as follows h W , b softmax ( x ( i ) ) = [ p ( y ( i ) = 2 ∣ x ( i ) ; W , b ) p ( y ( i ) = 2 ∣ x ( i ) ; W , b ) ⋮ p ( y ( i ) = K ∣ x ( i ) ; W , b ) ] = 1 Z [ exp { W 1 ( L − 1 ) · a ( L − 1 ) ( x ( i ) ) + b 1 ( L − 1 ) } exp { W 2 ( L − 1 ) · a ( L − 1 ) ( x ( i ) ) + b 2 ( L − 1 ) } ⋮ exp { W K ( L − 1 ) · a ( L − 1 ) ( x ( i ) ) + b K ( L − 1 ) } ] , where Z is the sum of all the entries, so that our probabilities sum to 1. Using this formulation, we can define our classification cost function: J λ softmax ( W , b ) = - 1 m ∑ i = 1 m ∑ j = 1 K 1 { y ( i ) = j } log p ( y ( i ) = j | x ( i ) ; W , b ) + λ 2 ∑ s = 1 u L - 1 ∑ t = 1 u L W s t ( L - 1 ) 2 . Intuitively, we can think about this cost function as making the log probability of the correct class as close to 0 as possible.

Transforming input data into summary statisticsterm·key For many deep learningterm·mth applications, the raw data can be used directly (the pixels of an image, for example). Unfortunately, we currently cannot input raw genomic data into a deep learningterm·mth method. Similarly to ABCterm·mth, we need to transform the data into summary statisticsterm·key that are potentially informative about the parameters of interest. Unlike ABCterm·mth, however, deep learningterm·mth should not be negatively affected by correlated or uninformative summary statisticsterm·key. Thus we sought to include a large number of potentially informative summary statisticsterm·key of the data. To account for the impact of selectionterm·dom, we divided each 100 kb region into three smaller regions: 1) close to the selected site (40–60 kb), 2) mid-range from the selected site (20–40 kb and 60–80 kb), and 3) far from the selected site (0–20 kb and 80–100 kb). These regions are based off of the simulation scenario in Peter et al. [ 23 ], and shown more explicitly in Fig 7 . Within each of these three regions, we calculated the following statistics, except where noted.

Regions used for computing the statistics, which are based off of Peter et al. [ 23 ]. Note that the selected site was chosen randomly within region 1.

For all the statistics described below, n is the haploid sample size. In the case of simulated data, n = 100. For the real data, we had 197 samples; within each 100 kb region we sorted the samples by missing data, then retained the 100 most complete samples, except where noted.

Number of segregating sites within each smaller region, S . Since all of the statistics must be in [0, 1], we normalized by S max = 5000 (any S > S max was truncated): 3 statistics.

Tajima’s D statistic [ 71 ], computed as follows D = π - S / a 1 Var ^ ( π - S / a 1 ) , where π is the average number of pairwise differences between two samples, and a 1 = ∑ i = 1 n - 1 1 / i . We normalized by D min = −3.0 and D max = 3.0, again truncating the rare statistic outside this range: 3 statistics.

Folded site frequency spectrum (SFS): η i is the number of segregating sites where the minor allele occurs i times out of n samples, for i = 1, 2, …, ⌊ n /2⌋. For the real data, for each segregating site, enough samples were included to obtain 100 with non-missing data. If that was not possible for a given site, the site was not included. To normalize the SFS, we divided each η i by the sum of the entries, which gives us the probability of observing i minor alleles, given the site was segregating: 50 ⋅ 3 = 150 statistics.

Length distribution between segregating sites: let B k be the number of bases between the k and k +1 segregating sites. To compute the distribution of these lengths, we define J bins and count the number of lengths that fall into each bin: d j BET = | k ∈ { 1 , 2 , ⋯ , S - 1 } : B k ∈ bin j | S . We choose J = 16 equally spaced bins, the first starting at 0 and the last starting at 300: 16 ⋅ 3 = 48 statistics.

Identity-by-state (IBS) tract length distribution: for each pair of samples, and IBS tract is a contiguous genomic region where the samples are identical at every base (delimited by bases where they differ). For all pairs, let L be the set of IBS tract lengths. In a similar fashion to the length distribution statistics, we define M bins, and count the number of IBS tract lengths that fall into each bin: d m IBS = | ℓ ∈ L : ℓ ∈ bin m | | L | . We choose M = 30 equally spaced bins, the first starting at 0 and the last starting at 5000: 30 ⋅ 3 = 90 statistics.

Linkage disequilibriumterm·mth (LDterm·mth) distribution: LDterm·mth is a measure of the correlation between two segregating sites. For example, if there was no recombination between two sites, their alleles would be highly correlated and LDterm·mth would be high in magnitude. If there was infinite recombination between two sites, they would be independent and have LDterm·mth close to 0. For two loci, let A , a be the alleles for the first site, and B , b be the alleles for the second site. Let p A be the frequency of allele A , p B be the frequency of allele B , and p AB be the frequency of haplotype AB . Then the linkage disequilibriumterm·mth is computed by D A B = p A B - p A p B . Here we compute the LDterm·mth between pairs of sites, where one site is in the “selected” region, and the other is in each of the three regions (including the selected region). Then we create an LDterm·mth distribution similar to the IBS distribution above, using 16 bins, with the first bin ending at −0.05, and the last bin starting at 0.2: 16 ⋅ 3 = 48 statistics.

H1, H12, and H2 statistics, as described in Garud et al. [ 72 ]. These statistics help to distinguish between hard and soft sweepsterm·key, and are calculated in the selected (middle) region only: 3 statistics.

This gives us a total of 345 statistics.

Deep learningterm·mth for demographyterm·key and selectionterm·dom Deep learningterm·mth model To modify our deep learningterm·mth method to accommodate this type of inference problem, during training we have an outer-loop that changes the demographyterm·key as necessary, and an inner loop that accounts for differences in selectionterm·dom for each region. During testing, we compare three different methods for predicting the effective population sizesterm·var: (1) Average stat : we average the statistics for each region within the dataset (corresponding to the same demographyterm·key), then run this average through the trained network. (2) Final : we run the statistics for each region through the trained network separately, then average the results. (3) Neutral regions : we use the same method as (2), but instead of using all the regions, we only use the ones we predict as being neutral. We also include a range between the 2.5th and 97.5th quantiles. To predict the selectionterm·dom class, we obtain a probability distribution over the classes for each region, then select the class with the highest probability.

One final complication is that we estimate continuous response variables for the population sizes, but consider selectionterm·dom to be a discrete response variable. This involves a linear activation function for the population sizes and a softmax classifier for selectionterm·dom. A diagram of our deep learningterm·mth method is shown in Fig 8 .

Our deep learningterm·mth framework for effective population sizeterm·var changes and selectionterm·dom.

Analysis of informative statistics One advantage of deep learningterm·mth is that the weights of the optimal network give us an interpretable link between the summary statisticsterm·key and the output variables. However, from these weights, it is not immediately obvious which statistics are the most informative or “best” for a given output variable. If we only wanted to use only a small subset of the statistics, which ones would give us the best results? To answer this question, we use both a permutation testing approach and a perturbation method described in S1 Algorithm .

Software availability An open-source software package (called evoNet) that implements deep learningterm·mth algorithms for population geneticterm·key inference is available at https://sourceforge.net/projects/evonet .

Supporting Information In this scenario, the testing data is simulated with a recombination rate that is 4 times higher than that of the training data. The effective population sizeterm·var results (top table) are still generally accurate, but selectionterm·dom (bottom table) is harder to predict. Neutral regions are predicted correctly, but soft sweepsterm·key are also often classified as neutral. The classifier has difficulty distinguishing between hard sweepsterm·key and balancing selectionterm·key. The top table should be compared to Table 1 and the bottom table should be compared to Table 3 .

(TIFF) Click here for additional data file.

In this scenario, the testing data is simulated with bottleneck parameters that are more severe than the training data. This has a slight negative impact on the population size results (largely on the most recent size which was outside the training range), but has little effect on the selectionterm·dom results. The overall percentage of misclassified regions is 7.8%. The top table should be compared to Table 1 and the bottom table should be compared to Table 3 .

(TIFF) Click here for additional data file.

(XLS) Click here for additional data file.

(XLS) Click here for additional data file.

The single hidden layer serves to learn informative combinations of the inputs, remove correlations, and typically reduce the dimension of the data. After the optimal weight on each connecting arrow is learned through labeled training data, unlabeled data can be fed through the network to estimate the response variables.

(TIFF) Click here for additional data file.

For each variable, the top 25 statistics were chosen, according to the procedure in S1 Algorithm . The Venn diagram captures statistics common to each subset of output variables, with notable less informative statistics shown in the lower right. Close, mid, and far represent the genomic region where the statistic was calculated. The numbers after each colon refer to the position of the statistic within its distribution or order. For the SFS statistics, it is number of minor alleles. For each region, there are 50 SFS statistics, 16 BET statistics (distribution between segregating sites), 30 IBS statistics, and 16 LDterm·mth statistics.

(TIFF) Click here for additional data file.

The x-axis shows increasing values of λ, and the y-axis shows the error on the validation dataset. The curve shows a characteristic shape with low and high λ producing poorer results than an intermediate value. For these hidden layers sizes and this dataset, λ ^ = 0 . 0001 was optimal.

(TIFF) Click here for additional data file.

The input data ( x ) is projected into a (usually) lower dimension ( a ), then reconstructed ( x ^ ). The weights of an autoencoder are optimized such that the difference between the reconstructed data and the original data is minimal.

(TIFF) Click here for additional data file.

(TIFF) Click here for additional data file.