Detecting adaptive introgression in human evolution using convolutional neural networks

Abstract
Introduction
Results
Discussion
Materials and methods
Appendix 1
Appendix 2
Appendix 3
Appendix 4
Appendix 5
Data availability
References
Article and author information
Metrics

Abstract

Studies in a variety of species have shown evidence for positively selected variants introduced into a population via introgression from another, distantly related population—a process known as adaptive introgression. However, there are few explicit frameworks for jointly modelling introgression and positive selection, in order to detect these variants using genomic sequence data. Here, we develop an approach based on convolutional neural networks (CNNs). CNNs do not require the specification of an analytical model of allele frequency dynamics and have outperformed alternative methods for classification and parameter estimation tasks in various areas of population genetics. Thus, they are potentially well suited to the identification of adaptive introgression. Using simulations, we trained CNNs on genotype matrices derived from genomes sampled from the donor population, the recipient population and a related non-introgressed population, in order to distinguish regions of the genome evolving under adaptive introgression from those evolving neutrally or experiencing selective sweeps. Our CNN architecture exhibits 95% accuracy on simulated data, even when the genomes are unphased, and accuracy decreases only moderately in the presence of heterosis. As a proof of concept, we applied our trained CNNs to human genomic datasets—both phased and unphased—to detect candidates for adaptive introgression that shaped our evolutionary history.

Introduction

Ancient DNA studies have shown that human evolution during the Pleistocene was characterised by numerous episodes of interbreeding between distantly related groups (Green et al., 2010; Reich et al., 2010; Meyer et al., 2012; Prüfer et al., 2017; Kuhlwilm et al., 2016). We now know, for example, that considerable portions of the modern human gene pool derive from Neanderthals and Denisovans (Green et al., 2010; Reich et al., 2010; Prüfer et al., 2014). In the past few years, several methods have been developed to identify regions of present-day or ancient human genomes containing haplotypes that were introgressed from other groups of hominins. These include methods based on probabilistic models (Sankararaman et al., 2014; Steinrücken et al., 2018; Racimo et al., 2017a), on summary statistics (Vernot and Akey, 2014; Vernot et al., 2016; Racimo et al., 2017b; Durvasula and Sankararaman, 2019) and on ancestral recombination graph reconstructions (Kuhlwilm et al., 2016; Hubisz et al., 2020; Speidel et al., 2019). Presumably, some of the introgressed material may have had fitness consequences in the recipient populations. While recent evidence suggests that a large proportion of Neanderthal ancestry was likely negatively selected (Harris and Nielsen, 2016; Juric et al., 2016), there is also support for positive selection on a smaller proportion of the genome—a phenomenon known as adaptive introgression (AI) (Whitney et al., 2006; Hawks and Cochran, 2006; Racimo et al., 2015).

Genomic evidence for AI has been found in numerous other species, including butterflies (Pardo-Diaz et al., 2012; Enciso-Romero et al., 2017), mosquitoes (Norris et al., 2015), hares (Jones et al., 2018), poplars (Suarez-Gonzalez et al., 2016), and monkeyflowers (Hendrick et al., 2016). A particularly striking example is AI in dogs, which appears to show strong parallels to AI in humans when occupying the same environmental niches. For example, a variant of the gene EPAS1 has been shown to have introgressed from an archaic human population into the ancestors of Tibetans, and subsequently risen in frequency in the latter population, as a consequence of positive selection to high altitude (Huerta-Sánchez et al., 2014). A different high-frequency EPAS1 variant is also uniquely found in Tibetan Mastiffs, and appears to also have introgressed into this gene pool via admixture with a different species, in this case Tibetan wolves (Miao et al., 2017), likely due to the same selective pressures.

To detect AI, researchers can look for regions of the genome with a particularly high frequency of introgressed haplotypes from a donor species or population into a recipient species or population. These haplotypes are often detected assuming neutrality of archaic alleles since the introgression event (Vernot et al., 2016; Vernot and Akey, 2014; Sankararaman et al., 2016). Other studies have designed statistics that are sensitive to characteristic patterns left by AI, using simulations incorporating both admixture and selection (Gittelman et al., 2016; Racimo et al., 2017b). More recently, Setter et al., 2020 developed a likelihood framework to look for local alterations to the site frequency spectrum that are consistent with adaptive introgression, using only data from the recipient species. The main challenge that these studies face is that it is hard to jointly model selection from material introduced via admixture (Racimo et al., 2015).

To overcome the need to compress data into summary statistics (which might miss important features), deep learning techniques are increasingly becoming a popular solution to address problems in population genetics. These problems include the inference of demographic histories (Sheehan and Song, 2016; Flagel et al., 2019; Villanea and Schraiber, 2019; Mondal et al., 2019; Sanchez et al., 2021), admixture (Blischak et al., 2021), recombination (Chan et al., 2018; Flagel et al., 2019; Adrion et al., 2020b), and natural selection (Schrider and Kern, 2018; Sheehan and Song, 2016; Torada et al., 2019; Isildak et al., 2021). Deep learning is a branch of machine learning that relies on algorithms structured as multi-layered networks, which are trained using known relationships between the input data and the desired output. They can be used for classification, prediction, or data compression (Aggarwal, 2018). Among the techniques in this field, convolutional neural networks (CNNs) are a family of methods originally designed for image recognition and segmentation (LeCun and Bengio, 1995; Krizhevsky et al., 2012), which have been recently applied to population genetic data (Chan et al., 2018; Flagel et al., 2019; Torada et al., 2019; Isildak et al., 2021; Blischak et al., 2021; Sanchez et al., 2021). A CNN can learn complex spatial patterns from large datasets that may be informative for classification or prediction, using a series of linear operations known as convolutions, to compress the data into features that are useful for inference.

Despite the recent advances in deep learning for population genetics, no dedicated attempts have been made to identify AI from population genomic data. Here, we develop a deep learning method called genomatnn that jointly models archaic admixture and positive selection, in order to identify regions of the genome under adaptive introgression. We trained a CNN to learn relevant features directly from a genotype matrix at a candidate region, containing data from the donor population, the recipient population and a unadmixed outgroup. The method has >88% precision to detect AI and is effective on both ancient and recently selected introgressed haplotypes. We then applied our method to population genomic datasets where the donor population is either Neanderthals or Denisovans and the recipient populations are Europeans or Melanesians, respectively. In each case, we used the Yoruba population as a unadmixed outgroup and we were able to both recover previously identified AI regions and unveil new candidates for AI in human history.

Results

A CNN for detecting adaptive introgression

We assume we have sequence data from multiple populations: the donor population and the recipient population in an admixture event, as well as an unadmixed population that is a sister group to the recipient (Figure 1). Our method relies on partitioning the genomes into windows, which we chose to be 100 kbp in size. For each window, we constructed an $n \times m$ matrix, where $n$ corresponds to the number of haplotypes (or diploid genotypes, for unphased data) and $m$ correspond to a set of equally sized bins along the genomic length of the window. Each matrix entry contains the count of minor alleles in an individual’s haplotype (or diploid genotype) in a given bin. Within each population, we sorted these pseudo-haplotypes (or genotypes) according to similarity to the donor population, and concatenated the matrices for each of the populations into a single pseudo-genotype matrix (Figure 1).

Figure 1 with 2 supplements see all

Download asset Open asset

A schematic overview of how genomatnn detects adaptive introgression.

We first simulate a demographic history that includes introgression, such as Demographic Model A1 shown in (A), using the SLiM engine in stdpopsim. Parameter values for this model are given in Appendix 3—table 1. Three distinct scenarios are simulated for a given demographic model: neutral mutations only, a sweep in the recipient population, and adaptive introgression. The tree sequence file from each simulation is converted into a genotype matrix for input to the CNN. (B) shows a genotype matrix from an adaptive introgression simulation, where lighter pixels indicate a higher density of minor alleles, and haplotypes within each population are sorted left-to-right by similarity to the donor population (Nea). In this example, haplotype diversity is low in the recipient population (CEU), which closely resembles the donor (Nea). Thousands of simulations are produced for each simulation scenario, and their genotype matrices are used to train a binary-classification CNN (C). The CNN is trained to output Pr[AI], the probability that the input matrix corresponds to adaptive introgression. Finally, the trained CNN is applied to genotype matrices derived from a VCF/BCF file (D).

We designed a CNN (Figure 1) that takes this concatenated matrix as input to distinguish between adaptive introgression scenarios and other types of neutral or selection scenarios. The CNN was trained using simulations, and uses a series of convolution layers with successively smaller outputs, to extract increasingly higher level features of the genotype matrices—features which are simultaneously informative of introgression and selection. The CNN outputs the probability that the input matrix comes from a genomic region that underwent adaptive introgression. As our simulations used a wide range of selection coefficients and times of selection onset, the network does not assume these parameters are known a priori, and is able to detect complete or incomplete sweeps at any time after gene flow.

Our method has several innovative features relative to previous population genetic implementations of CNNs (described extensively in the Materials and methods section). For example, when loading the genotype matrices as input, we implemented an image resizing scheme that leads to fast training times, while avoiding the drawbacks of similar methods (Torada et al., 2019), by preserving inter-allele distances and thus the local density of segregating sites. Additionally, instead of using pooling layers, we used a 2 × 2 step size when performing convolutions. This has the same effect as pooling, in that the output size of the layer is smaller than the input, so the accuracy of the model is comparable to traditional implementations of CNNs, but it has a much lower computational burden (Springenberg et al., 2015).

Furthermore, we incorporated a framework to visualise the features of the input data that draw the most attention from the CNN, by plotting saliency maps (Simonyan et al., 2014). Saliency maps can help to understand which regions of the genotype matrix contribute the most towards the CNN prediction score.

We also provide downloadable pre-trained CNNs as well as a pipeline for training new CNNs (see Materials and methods). These interface with a new selection module that we designed and incorporated into the stdpopsim framework (Adrion et al., 2020a), using the forwards-in-time simulator SLiM (Haller and Messer, 2019). We believe this will facilitate the application of the method to other datasets, allowing users to modify its parameters according to the specific requirements of the biological system under study.

Performance on simulations

We aimed to assess the performance of our method on simulations. We performed simulations under two main demographic models:

Demographic Model A1: a three-population model including an African, a European and a Neanderthal population, with Neanderthal gene flow into Europeans (Figure 1 and Figure 1—figure supplement 1)
Demographic Model B: a more complex model, including an African, a Melanesian, a Neanderthal and a Denisovan population, with two pulses of Denisovan gene flow into Melanesians, plus Neanderthal gene flow into non-Africans, based on Jacobs et al., 2019 (Figure 1—figure supplement 2).

When training a CNN on Demographic Model A1 using phased data, we obtained a precision of 90.2% (the proportion of AI predictions that were AI simulations) and 97.9% negative predictive value (NPV; the proportion of ‘not-AI’ predictions that were either neutral or sweep simulations) (Figure 2). The network output higher probabilities for AI simulations with larger selection coefficients, and for older times of onset of selection. We also observed that the network falsely classified neutral simulations as AI more frequently than it falsely classified sweep simulations. When the CNN was trained on this same demographic model assuming genotypes were unphased, the results were very similar, with 88.1% precision and 98.7% NPV.

Figure 2 with 2 supplements see all

Download asset Open asset

CNN performance on validation simulations for Demographic Model A.

The CNN was trained using only AI simulations with selected mutation having allele frequency >0.25. (A) Confusion matrix. For the two prediction categories, either 'not AI' or AI, we show the proportion attributed to each of the true (simulated) scenarios. (B) Average CNN prediction for AI scenarios, binned by selection coefficient, $s$ , and time of onset of selection $T_{s e l}$ . (C) ROC curves, precision-recall curves and MCC-F₁ curves. The positive condition is AI. The negative conditions are shown using different line styles/colours. The circles indicate the point in ROC-space (respectively Precision-Recall-space, and MCC-F₁-space) when using the threshold Pr[AI]>0.5 for classifying a genotype matrix as AI. DFE: distribution of fitness effects. TP: true positives; FP: false positives; TN: true negatives; FN: false negatives; TPR: true positive rate; FPR: false positive rate; ROC: Receiver operating characteristics; MCC: Mathews correlation coefficient; F₁: harmonic mean of precision and recall.

When training a CNN on Demographic Model B (assuming unphased genotypes, as accurately phased data are not readily available for Melanesian genomes), we obtained 88.8% precision and 82.5% NPV (Figure 2—figure supplement 1). We note here that the network had greater precision when detecting AI derived from the more ancient pulse of Denisovan gene flow than the younger pulse.

Kim et al., 2018 and Zhang et al., 2020 recently suggested that introduced genetic material can mask deleterious recessive variation and produce a signal very similar to adaptive introgression. To assess whether heterosis following introgression affects the false positive rates in our CNN, we simulated a distribution of fitness effects (DFE) with recessive dominance for 70% of derived mutations (the rest were simulated as neutral), and found this only slightly increases the false positive rate (Figure 2).

MCC-F₁ curve

While precision, recall, and false positive rate are informative, these each consider only two of the four cells in a confusion matrix (true positives, true negatives, false positives, false negatives), and may produce a distorted view of performance with imbalanced datasets (Chicco, 2017). To obtain a more robust performance assessment, we plotted the Matthews correlation coefficient (MCC; Matthews, 1975) against F₁-score (the harmonic mean of precision and recall) for false-positive-rate thresholds from 0 to 100 (Figure 2, Figure 2—figure supplement 1, Figure 2—figure supplement 2), as recently suggested by Cao et al., 2020. MCC produces low scores unless a classifier has good performance in all four confusion matrix cells, and also accounts for class imbalance. In MCC-F₁ space, the point (1, 1) indicates perfect predictions, and values of 0.5 for the (unit-normalised) MCC indicate random guessing. These results confirm our earlier findings, that the CNN performance is excellent for Demographic Model A1 when considering either neutral and sweep simulations as the condition negative, and performance decreases slightly when DFE simulations are the negative condition (Figure 2). Furthermore, the CNN performance is not as good for Demographic Model B, but this is unlikely to be caused by using unphased genotypes (Figure 2—figure supplement 1 and Figure 2—figure supplement 2).

Comparison to other methods

We compared the performance of our CNN to VolcanoFinder (Setter et al., 2020), which scans genomes of the recipient population for patterns of diversity indicative of AI using a coalescent-based model of adaptive introgression (Figure 2—figure supplement 2). However, this method only incorporates information from a single population and is designed to detect 'ghost' introgression in cases where the source is not available. VolcanoFinder performed poorly for the demographic models considered here—in some cases, worse than guessing randomly. We also compared our CNN to an outlier-based approach for a range of summary statistics that are sensitive to AI (Racimo et al., 2017b). Our CNN is closest to a perfect MCC-F₁ score for Demographic model A1 and B, closely followed by the Q95(1%, 100%) and then U(1%, 20%, 100%) statistics developed in Racimo et al., 2017b.

Demographic model misspecification

We then tested robustness to demographic misspecification, by evaluating the CNN trained on Demographic Model A1 against simulations for two other demographic models (Figure 2—figure supplement 2). We considered weak misspecification, where the true demographic history is similar to Demographic Model A1 but also includes archaic admixture within Africa following Ragsdale and Gravel, 2019 (Demographic Model A2; Figure 1—figure supplement 1). This resulted in only a small performance reduction. We also considered strong misspecification, where the true demographic history is Demographic Model B. As there are more Melanesian individuals than European individuals in our simulations (because we aimed to mimic the real number of genomes available in our data analysis below), we downsampled the Melanesian genomes to match the number of European genomes, so as to perform a fair misspecification comparison. In this case, the performance of the CNN was noticeably worse than that of the summary statistics, but still better than VolcanoFinder. We note that the summary statistics performance decreased also, to match their performance for the correctly-specified assessments on Demographic Model B. Interestingly, we found that the Q95(1%, 100%) statistic was the most robust method for both cases of misspecification.

Network attention

To understand which features of the input matrices were used by the CNN to make its predictions, we constructed saliency maps (Simonyan et al., 2014). This technique works by computing the gradient of a network’s output with respect to a single input. Thus, highlighted regions from the saliency map indicate where small changes in the input matrix have a relatively large influence over the CNN output prediction.

We calculated an average saliency map for each simulation scenario (neutral, sweep, or AI), for a CNN trained on Demographic Model A1 (Figure 3). Our results show that when the network was presented with an AI matrix, it focused most of the attention on the Neanderthal and European haplotypes, rather than on the African haplotypes. In non-AI scenarios, the network focused sharply on the Neanderthal and left-most European haplotypes. The saliency maps also show a concentration of attention in the central region of the genomic window, closer to where the selected mutation was drawn (even though this mutation does not exist in neutral simulations, and was removed from sweep and AI simulations before constructing genotype matrices; see Methods). We also note a periodic vertical banding pattern in the saliency maps, corresponding to the filter width for the convolution layers.

Figure 3

Download asset Open asset

Saliency maps, showing the CNN’s attention across the input matrices for each simulated scenario, calculated for the CNN trained on Demographic Model A, filtered for beneficial allele frequency >0.25.

Each panel shows the average gradient over 300 input matrices encoding either neutral (top), sweep (middle), or AI (bottom) simulations. Pink/purple colours indicate larger gradients, where small changes in the genotype matrix have a relatively larger influence over the CNN’s prediction. Columns in the input matrix correspond to haplotypes from the populations labelled at the bottom.

Calibration

We implemented a score calibration scheme to account for the fact that our simulation categories (neutrality, sweep, and AI) will be highly imbalanced in real data applications (Guo et al., 2017; Kull et al., 2017). CNN classifiers sometimes produce improperly calibrated probabilities (Guo et al., 2017). In our case, this occurs because the proportion of each category is not known in reality, and thus does not match the simulated proportion. For this reason, we fitted our calibration procedure using training data resampled with various ratios of neutral:sweep:AI simulations (Figure 4). We tested different calibration methods by fitting the calibrator to the training dataset, and inspecting reliability plots and the sum of residuals on a validation dataset (see Materials and methods, Figure 4—figure supplement 1, Figure 4—figure supplement 2, Figure 4—figure supplement 3, Figure 4—figure supplement 4). When the probabilities are calibrated for even class ratios, Manhattan plots show a large number of high probability candidates across the genome, which obscure the strongest peaks (Figure 4). However, once calibrated for class ratios that are skewed towards the neutral class, strong candidates for AI become more apparent.

Figure 4 with 4 supplements see all

Download asset Open asset

Comparison of Manhattan plots using beta-calibrated output probabilities for different class ratios.

Each row indicates a single CNN, with equivalent data filtering. Each column indicates different class ratios used for calibration (Neutral:Sweep:AI). AF = Minimum beneficial allele frequency.

Candidates for Neanderthal adaptive introgression in European genomes

We applied our method to a combined genomic panel of archaic hominins (Prüfer et al., 2017; Meyer et al., 2012) and present-day humans (The 1000 Auton et al., 2015; Jacobs et al., 2019), to look for regions of the genome where Non-African humans show signatures of AI from archaic hominins. First, we looked for Neanderthal introgression into the ancestors of Northwestern Europeans (CEU panel), using Yoruba Africans (YRI panel) as the unadmixed sister population. To give the network the best chance of avoiding false positives, we tried two different beneficial-allele frequency cutoffs for training: 5% and 25% (Table 1 and Appendix 1—table 1).

Table 1

Top ranking gene candidates corresponding to Neanderthal AI in Europeans.

We show genes which overlap, or are within 100 kbp of, the 30 highest ranked 100 kbp intervals. Adjacent intervals have been merged. The CNN was trained using only AI simulations with selected mutation having allele frequency >0.25, and subsequently calibrated with resampled neutral:sweep:AI ratios of 1:0.1:0.02.

Chrom	Start	End	Genes
1	104500001	104600000
2	109360001	109460000	LIMS1; RANBP2; CCDC138; EDAR
2	160160001	160280000	TANC1; WDSUB1; BAZ2B
3	114480001	114620000	ZBTB20
4	54240001	54340000	SCFD2; FIP1L1; LNX1
5	39220001	39320000	FYB; C9; DAB2
6	28180001	28320000	ZSCAN16-AS1; ZSCAN16; ZKSCAN8; ZSCAN9; ZKSCAN4; NKAPL; PGBD1; ZSCAN31; ZKSCAN3; ZSCAN12; ZSCAN23
8	143440001	143560000	TSNARE1; BAI1
9	16700001	16820000	BNC2
12	85780001	85880000	ALX1
19	20220001	20380000	ZNF682; ZNF90; ZNF486
19	33580001	33740000	RHPN2; GPATCH1; WDR88; LRP3; SLC7A10
20	62100001	62280000	CHRNA4; KCNQ2; EEF1A2; PPDPF; PTK6; SRMS; C20orf195; HELZ2; GMEB2; STMN3; RTEL1; TNFRSF6B; ARFRP1; ZGPAT; LIME1; SLC2A4RG; ZBTB46
21	25840001	25940000

We focus here on describing the results from the 25% condition (Figure 5 and Appendix 4). We found several candidate genes for AI that have been reported before (Sankararaman et al., 2014; Vernot and Akey, 2014; Gittelman et al., 2016; Racimo et al., 2017b), including BNC2, KCNQ2/EEF1A2 WRD88/GPATCH1 and TANC1. Notably, the candidate region we identify on chromosome 2 around TANC1 extends farther downstream of this gene, also overlapping BAZ2B (Appendix 4—figure 3). This codes for a protein related to chromatin remodelling, and may have a role in transcriptional activation. Mutations in BAZ2B have recently been associated with neurodevelopmental disorders, including developmental delay, autism spectrum disorder and intellectual disability (Scott et al., 2020).

Figure 5

Download asset Open asset

Application of the trained CNN to the Vindija and Altai Neanderthals, and 1000 genomes populations YRI and CEU.

The CNN was applied to overlapping 100 kbp windows, moving along the chromosome in steps of size 20 kbp. The CNN was trained using only AI simulations with selected mutation having allele frequency > 25%, and subsequently calibrated with resampled neutral:sweep:AI ratios of 1:0.1:0.02.

Additionally, we found two novel candidates for AI that have not been previously reported, spanning the regions chr6:28.18 Mb–28.32 Mb (Appendix 4—figure 7) and chr20:62.1 Mb–62.28 Mb (Appendix 4—figure 13), including multiple genes encoding zinc finger proteins. UK-biobank PheWAS associations (Canela-Xandri et al., 2018) suggest both regions generally affect phenotypes related to blood, including platelet, erythrocyte and leukocyte counts (at the $p < 10^{- 8}$ association level, the chr6 region has 91 hits, while the chr20 region has 19, with 10 of these traits in common).

Candidates for Denisovan adaptive introgression in Melanesian genomes

We then looked for Denisovan AI in Melanesian genomes from the IGDP panel (Jacobs et al., 2019), also considering Yoruba Africans as the unadmixed sister group, and using two different beneficial-allele frequency cutoffs for training: 5% and 25% (Table 2 and Appendix 1—table 2).

Table 2

Top ranking gene candidates corresponding to Denisovan AI in Melanesians.

Chrom	Start	End	Genes
2	129960001	130060000
3	3740001	3840000	SUMF1; LRRN1
4	41980001	42080000	TMEM33; DCAF4L1; SLC30A9; BEND4
5	420001	520000	PDCD6; AHRR; C5orf55; EXOC3; CTD-2228K2.5; SLC9A3; CEP72
6	74640001	74740000
6	81960001	82060000
6	137920001	138120000	TNFAIP3
7	25100001	25200000	OSBPL3; CYCS; C7orf31; NPVF
7	38020001	38120000	EPDR1; NME8; SFRP4; STARD3NL
7	121160001	121260000
8	3040001	3140000	CSMD1
12	84640001	84740000
12	108240001	108340000	PRDM4; ASCL4
12	114020001	114280000	RBM19
14	61860001	61960000	PRKCH
14	63120001	63220000	KCNH5
14	96700001	96820000	BDKRB2; BDKRB1; ATG2B; GSKIP; AK7
15	55260001	55400000	RSL24D1; RAB27A
16	62600001	62700000
16	78360001	78460000	WWOX
18	22060001	22160000	OSBPL1A; IMPACT; HRH4
22	19040001	19140000	DGCR5; DGCR2; DGCR14; TSSK2; GSC2; SLC25A1; CLTCL1

Again, we focus on describing the results from the 25% condition (Figure 6 and Appendix 5) Among the top candidates, we found a previously reported candidate for AI in Melanesians: TNFAIP3 (Vernot et al., 2016; Gittelman et al., 2016). Denisovan substitutions carried by the introgressed haplotype in this gene have been found to enhance the immune response by tuning the phosphorylation of the encoded A20 protein, which is an immune response inhibitor (Zammit et al., 2019).

Figure 6

Download asset Open asset

Application of the trained CNN to the Altai Denisovan and Altai Neanderthal, 1000 genomes YRI populations, and IGDP Melanesians.

The CNN was applied to overlapping 100 kbp windows, moving along the chromosome in steps of size 20 kbp. The CNN was trained using only AI simulations with selected mutation having allele frequency > 25%, and subsequently calibrated with resampled neutral:sweep:AI ratios of 1:0.1:0.02.

We found evidence for Denisovan AI in Melanesians at several other candidate regions. A few of these regions (or contiguous regions) were previously reported by Sankararaman et al., 2016 but not extensively described, possibly because the previously reported sections of those regions deemed to be introgressed were intergenic. One of the regions with strong evidence for AI (chr7:25.1 Mb–25.2 Mb; Appendix 5—figure 8) overlaps the CYCS gene. This gene codes for cytochrome C: a small heme protein that plays a crucial role in the electron transport chain in mitochondria, and has been associated with various blood-related diseases, like thrombocytopenia (Morison et al., 2008; De Rocco et al., 2014; Uchiyama et al., 2018). Another top candidate region (chr12:108.24–108.34 Mb, Appendix 5—figure 13) is upstream of PRDM4 and ASCL4. The former gene codes for a transcription factor that may be involved in the nerve growth factor cell survival pathway and might play a role in tumour suppression (Yang and Huang, 1999). The latter gene codes for a different transcription factor that may be involved in skin development (Jonsson et al., 2004).

We detected signatures of Denisovan AI in a region in chromosome 3 near SUMF1 and LRNN1 (Appendix 5—figure 2), which was also identified in Jacobs et al., 2019. SUMF1 codes for an enzyme involved in the hydrolysis of sulfate esters, which has been associated with sulfatase deficiency (Cosma et al., 2003). LRNN1 encodes a protein involved in neuronal differentiation, which has been associated with neuroblastoma and Alzheimer’s disease (Bai et al., 2014; Hossain et al., 2012). Another candidate region is in chromosome 7 and is upstream of SFRP4 (Appendix 5—figure 9), which encodes a protein associated with diabetes (Mahdi et al., 2012) and Pyle’s disease (Kiper et al., 2016). Moreover, there is also a candidate region upstream of RAB27A, in chromosome 15 (Appendix 5—figure 18). Mutations in this gene cause Griscelli syndrome, which results in pigmentary dilution in the hair and skin, as well as melanosome accumulation in melanocytes (Ménasché et al., 2000). Finally, we found evidence for Denisovan AI in two nearby regions in chromosome 14 (Appendix 5—figure 15 and Appendix 5—figure 16). One of these overlaps with PRKCH—encoding a protein kinase associated with cerebral infarction (Kubo et al., 2007). The other overlaps with KCNH5—coding for a potassium channel that may be associated with epileptic encephalopathy (Veeramah et al., 2013).

Discussion

We have developed a new method to detect adaptive introgression along the genome using convolutional neural networks. The method has high precision when reporting candidate AI loci, and high negative predictive value when rejecting loci as not-AI: we obtain greater than 90% accuracy under a variety of different selection scenarios (Appendix 2—table 1), with low false positive rates.

As reported previously (Kim et al., 2018; Zhang et al., 2020), heterosis following introgression can produce patterns very similar to AI, and we found this can inflate false positive detection of AI by our CNN to a small extent. However, we simulated a DFE with recessive dominance for all mutations, which is not realistic in general, so our results in this regard represent a worst-case scenario. A possible future improvement would be to train the CNN on simulations incorporating heterosis. We did not attempt this here because realistic DFE simulations represent a substantial computational burden.

When the demographic model is correctly specified, our CNN performed better than using summary statistics, although the Q95(1%, 100%) statistic (Racimo et al., 2017b) also performed remarkably well. This statistic captures high-frequency derived alleles that are shared between the donor and recipient population, to the exclusion of a non-introgressed sister population—intuitively, these are the same features we expect our CNN to be learning. Because of its relative robustness to model misspecification, an outlier approach based on Q95 may be a better choice than our CNN when there is uncertainty regarding the demographic history of the study system. We also found that VolcanoFinder performed very poorly across all our tests, but this is arguably an unfair comparison because it only incorporates information from a single population, and Setter et al., 2020 themselves found that their method has low sensitivity when the donor population split from the recipient population recently (on the order of $N$ generations ago for Neanderthals/Denisovans and humans).

The CNN took approximately 15 min to train on one NVIDIA Tesla T4 GPU, which amounts to 60 CPU hours for an equivalent CPU-only training procedure. All data were loaded into memory, which required approximately 120 GB RAM during training. The computational bottleneck lay in the generation of SLiM forward simulations: 300,000 simulations took approximately 80 weeks of CPU time for each of demographic models A1 and B. In the future, considerable speedups could potentially be obtained by optimising the simulation step, perhaps by implementing an adaptive introgression simulation framework that takes advantage of the backwards coalescent (e.g. building on the work by Setter et al., 2020).

We applied the method to human data, to look for adaptive introgression from archaic humans into the ancestors of present-day human genomes. When looking for Neanderthal AI in European genomes, we found previously reported candidate genes (BNC2, WRD88/GPATCH1, KCNQ2/EEF1A2, TANC1/BAZ2B). We also recovered candidates for adaptive introgression from Denisovans by applying our method to unphased Melanesian genomes. The top candidates include TNFAIP3, which has been reported before, but also include other, novel regions, containing genes involved in blood diseases (CYCS), neurological diseases (PRKCH, KCNH5, LRNN1), metabolism (SFRP4, SUMF1), and skin development (ASCL4, RAB27A).

We note, however, that, as with previous methods, visual inspection of the haplotypes or genotypes of the top candidate regions remains a necessary criterion to accurately assess whether a region may have been under adaptive introgression. For example, in the scans we performed, we found a few candidate regions for Neanderthal AI in Europeans that are likely to be false positives, for example chr2:109360001–109460000 (Appendix 4—figure 2); chr4:54240001–54340000 (Appendix 4—figure 5); chr8:143440001–143560000 (Appendix 4—figure 8). These appear to be the result of shared ancestral variation between European and African populations, and yet are classified as having high probability of being under AI. These regions also appear to be generally low in diversity, which is possibly a result of data missingness rather than low diversity per se. Thus, our method allows for a rapid scan and prioritisation of potential targets, but these need to be further assessed with care. Inclusion of more complex selection scenarios, involving positive or balancing selection on ancestral variation, as well as linked selection, might serve to ameliorate the rate of false positives in the future. Furthermore, our simulation procedure does not model genotype errors or variation in data missingness. Not explicitly accounting for this may negatively impact the robustness of the minor allele density computation and the subsequent haplotype sorting procedure, and, in turn, affect the accuracy of the CNN.

Conversely, there may be regions under AI that are classified as highly probable by the CNN, but that did not appear in our top candidates. Validating a large number of candidates might be difficult, but one could imagine running a differently trained CNN (perhaps one better tailored to distinguish AI from more similar scenarios, like selection on shared ancestral variation) on the subset of the regions that are predicted to be AI using a lenient probability cutoff. One could also use our method more generally, to assess the impact of AI across the genome, by comparing the distribution of probability scores with those of simulation scenarios under different amounts of admixture and selection, though in that case one would need to train the CNN on a wider range of admixture rates and demographic models.

The performance of our method necessarily depends upon the demographic history of the populations involved. We found it more challenging to detect AI when the timing of gene flow is younger or the introgressing population is more diverged from the panel that is used to represent it. This is apparent when comparing results for the Neanderthal-into-European demographic scenario and the Denisovan-into-Melanesian demographic scenario. In the former, gene flow is older (∼55 kya versus ∼50 kya and ∼30 kya) (Sankararaman et al., 2016; Jacobs et al., 2019) and sequences are available for a population closely related to the putative source, which increases power. Furthermore, for the two putative pulses of Denisovan gene flow (Jacobs et al., 2019), we find our model has greater recall with AI for the more ancient pulse (94% versus 83.6%; Figure 2—figure supplement 1), likely because haplotypes from the older pulse have more time to rise in frequency. Similarly, recall is diminished when the onset of selection is more recent. We also found that distinguishing AI from a selective sweep (hard or soft), is relatively easier than distinguishing AI from neutral variation.

Our method requires sequencing data from the population from which the introgression event originated. This may be problematic in cases where the source of introgression may be distantly related to the population genomic panel that is used to represent it. Future work could involve developing a CNN that can detect adaptive introgression from a ghost (unsampled) population, for cases in which genomic data from the source is unavailable (e.g. see Setter et al., 2020).

The method can take either phased or unphased data as input. This flexibility allows for its application to a range of study systems in the future, in which phasing may not be financially or methodologically feasible. It does, however, require called genotypes and is therefore not yet suitable for genomes sequenced at low coverage. One could envision extending the framework developed here to low-coverage genomes by working with matrices of genotype likelihoods (Korneliussen et al., 2014) rather than matrices of genotypes or haplotypes. Flagel et al., 2019, for example, developed a CNN to infer recombination rates in tetraploids without genotype calls, using read pileup information. A CNN could learn the relationship between observed read counts or genotype likelihoods under a given adaptive introgression scenario and the model parameters that can generate that data, but we leave that to a future work.

Future studies could also address the fact that we must use simulations to train the network, which involves an implicit amount of supervision by the user. The range of parameters and models that are simulated during training are necessarily specified a priori, and misspecification can negatively affect CNN performance. Progress in this regard could involve the use of generative adversarial networks (GANs), which appears to be a fruitful way to address this. Indeed, recent work suggests that one can train a GAN to learn to generate realistic population genomic data for any population (Wang et al., 2020).

The attention analyses performed here allowed only a posteriori reasoning on how the network learned to predict AI, so further work is needed in this area. For instance, interpretability of neural networks can be assessed using symbolic metamodelling (Alaa and van der Schaar, 2019) with reinforcement learning algorithms deployed to identify the subset of most informative features of input data (Yoon et al., 2019). In this context, such approaches should be able to pinpoint the important characteristics of genomic data, and possibly derive more informative summary statistics to predict complex evolutionary events.

In summary, we have shown that CNNs are a powerful approach to detecting adaptive introgression and can recover both known and novel selection candidates that were introduced via admixture. As in previous applications to other problems in the field (Sheehan and Song, 2016; Flagel et al., 2019; Schrider and Kern, 2018; Villanea and Schraiber, 2019; Mondal et al., 2019; Torada et al., 2019; Isildak et al., 2021), this exemplifies how deep learning can serve as a very powerful tool for population genetic inference. This type of technique may thus be a useful resource for future studies aiming to unravel our past history and that of other species, as statistical methodologies and computational resources continue to improve.

Materials and methods

Simulations

Request a detailed protocol

For CNN training, we performed simulations under three scenarios: neutral mutations only; positive selection of a de novo mutation in the recipient population (selective sweep); and positive selection of a derived mutation that was transferred via gene flow from the donor population to the recipient population (adaptive introgression, AI). In the sweep and AI scenarios, the selection coefficient was drawn log-uniformly from between 0.0001 and 0.1 for Europeans and between 0.001 and 0.1 for Melanesians (in the latter case, very few selected alleles survive with a very small selection coefficient, so we narrowed the range to reduce computational burden). The uniformly distributed time of mutation was decoupled from the uniformly distributed time of selection onset, thus allowing for soft sweeps (Hermisson and Pennings, 2005). For the selective sweep scenario, the mutation and selection times could occur at any time older than 1 kya but more recent than the split between the recipient population and its unadmixed sister population, with the constraint that the mutation must be introduced before the onset of selection. For the AI scenario, a neutrally evolving mutation was introduced to the donor population any time more recent than the split between the donor and the ancestor of recipient and unadmixed sister population, but older than 1 kya before the introgression event. Then, this mutation was transmitted to the recipient population, whereupon selection could start to act on it at any time after introgression but before 1 kya.

We further evaluated our Demographic Model A1 CNNs using an additional 10,000 simulations that incorporated a DFE using the parameters estimated for Europeans in Kim et al., 2017 and used in Kim et al., 2018. We considered two mutation types: 30% neutral and 70% deleterious. The deleterious portion of introduced mutations had a selection coefficient drawn from a reflected gamma distribution with shape parameter 0.186, and expected value −0.01314833. We approximated the dominance scheme from Kim et al., 2018, using a fixed dominance coefficient for deleterious mutations of $0.5 / (1 - 7071.07 * E [s])$ where $E [s]$ is the expected value from the gamma distribution (i.e. all deleterious mutations were effectively recessive).

To incorporate selection, we implemented a new module in stdpopsim (Adrion et al., 2020a), which leverages the forwards-in-time simulator SLiM (Haller and Messer, 2019) for simulating selection. For consistency, we also used stdpopsim’s SLiM engine for neutral simulations. stdpopsim uses SLiM’s ability to output tree sequences (Haller et al., 2019; Kelleher et al., 2018), which retains complete information about the samples’ marginal genealogies. Further, stdpopsim recapitates the tree sequences (ensuring that all sampled lineages have a single common ancestor), and applies neutrally evolving mutations to the genealogies, using the coalescent framework of msprime (Kelleher et al., 2016).

We simulated 100 kbp regions, with a mutation rate of $1.29 \times 10^{- 8}$ per site per generation (Tian et al., 2019), an empirical recombination map drawn uniformly at random from the HapMapII genetic map (Frazer et al., 2007), and the selected mutation introduced at the region’s midpoint. For both the sweep scenario and the AI scenario, we used a rejection-sampling approach to condition on the selected allele’s frequency being ≥ in the recipient population at the end of the simulation. This was done by saving the simulation state prior to the introduction of the selected mutation (and saving again after successful transmission to the recipient population, for the AI scenario), then restoring simulations to the most recent save point if the mutation was lost, or if the allele frequency threshold was not met at the end of the simulation.

To speed up simulations, we applied a scaling factor of $Q = 10$ . Scaling divides population sizes ( $N$ ) and event times ( $T$ ) by $Q$ , and multiplies the mutation rate µ, recombination rate $r$ and selection coefficient $s$ by $Q$ , such that the population genetic parameters $θ = 4 N μ$ , $ρ = 4 N r$ , and $N s$ remain approximately invariant to the applied scaling factor (Haller and Messer, 2019). After simulating, we further filtered our AI scenario simulations to exclude those that ended with a minor beneficial allele frequency less than a specific cutoff. We tried two cutoffs—5% and 25%—and present results for both. Rejection sampling within SLiM was not possible at these higher thresholds, as simulations often had low probability of reaching the threshold, particularly for recently introduced mutations. We note that this post-simulation filtering alters the distributions of selection coefficients and times of selection onset used for CNN training.

To investigate Neanderthal gene flow into Europeans, we simulated an out-of-Africa demographic model with a single pulse of Neanderthal gene flow into Europeans but not into African Yoruba (Demographic Model A1, Figure 1—figure supplement 1), using a composite of previously published model parameters (Appendix 3—table 1). The number of samples to simulate for each population was chosen to match the YRI and CEU panels in the 1000 Genomes dataset (Auton et al., 2015), and the two high coverage Neanderthal genomes (Prüfer et al., 2014). The two simulated Neanderthals were sampled at times corresponding to the estimated ages of the samples as reported in Prüfer et al., 2017. To test model misspecification, we performed an additional 10,000 simulations per simulation scenario on a modified version of this model that also incorporates archaic admixture in Africa (Ragsdale and Gravel, 2019) (Demographic Model A2; Figure 1—figure supplement 1).

To investigate Denisovan gene flow into Melanesian populations, we simulated an out-of-Africa demographic history incorporating two pulses of Denisovan gene flow (Malaspinas et al., 2016; Jacobs et al., 2019) implemented as the PapuansOutOfAfrica_10J19 model in stdpopsim (Adrion et al., 2020a). For this demographic model we sampled a single Denisovan and a single Neanderthal (with sampling time of the latter corresponding to the Altai Neanderthal’s estimated age). The number of Melanesian samples was chosen to match a subset of the IGDP panel (Jacobs et al., 2019). The Baining population of New Britain was excluded at the request of the IGDP data access committee, and we also excluded first-degree relatives, resulting in a total of 139 Melanesian individuals used in the analysis. As this demographic model includes two pulses of Denisovan admixture, we simulated half of our AI simulations to correspond with gene flow from the first pulse, and half from the second pulse.

Conversion of simulations to genotype matrices

Request a detailed protocol

We converted the tree sequence files from the simulations into genotype matrices using the tskit Python API (Kelleher et al., 2016). Major alleles (those with sample frequency greater than 0.5 after merging all individuals) were encoded in the matrix as 0, while minor alleles were encoded as 1. In the event of equal counts for both alleles, the major allele was chosen at random. Only sites with a minor allele frequency >5% were retained. For sweep and AI simulations, we excluded the site of the selected mutation.

We note that different simulations result in different numbers of segregating sites, but a constraint for efficient CNN training is that each datum in a batch must have the same dimensions. Existing approaches to solve this problem are to use only a fixed number of segregating sites (Chan et al., 2018), to pad the matrix out to the maximum number of observed segregating sites (Flagel et al., 2019), or to use an image-resize function to constrain the size of the input data (Torada et al., 2019). Each approach discards spatial information about the local density of segregating sites, although this may be recovered by including an additional vector of inter-site distances as input to the network (Flagel et al., 2019).

To obtain the benefits of image resizing (fast training times for reduced sizes and easy application to genomic windows of a fixed size), while avoiding its drawbacks, we chose to resize our input matrices differently, and only along the dimension corresponding to sites. To resize the genomic window to have length $m$ , the window was partitioned into $m$ bins, and for each individual haplotype we counted the number of minor alleles observed per bin. Compared with interpolation-based resizing (Torada et al., 2019), binning is qualitatively similar, but preserves inter-allele distances and thus the local density of segregating sites. Furthermore, as we do not resize along the dimension corresponding to individuals, this also permits the use of permutation-invariant networks (Chan et al., 2018), although we do not pursue that network architecture here.

We report results for $m = 256$ , but also tried $m = 32$ , 64, and 128 bins. Preliminary results indicated greater training and validation accuracy for CNNs trained with more bins, around 1% difference between both 32 and 64, and 64 and 128, although only marginal improvement for 256 compared with 128 bins. When matching unphased data, we combined genotypes by summing minor allele counts between the chromosomes of each individual. We note that all data were treated as either phased, or unphased, and no mixed phasing was considered.

We then partitioned the resized genotype matrix into submatrices by population. Submatrices were ordered left-to-right according to the donor, recipient, and unadmixed populations respectively. For genotype matrices including both Neanderthals and Denisovans, we placed the non-donor archaic population to the left of the donor. To ensure that a non-permutation-invariant CNN could learn the structure in our data, we sorted the haplotypes (Flagel et al., 2019; Torada et al., 2019). The resized haplotypes/individuals within each submatrix were ordered left-to-right by decreasing similarity to the donor population, calculated as the Euclidean distance to the average minor-allele density of the donor population (analogous to a vector of the donor allele frequencies). An example (phased) genotype matrix image for an AI simulation is shown in Figure 1.

Conversion of empirical data to genotype matrices

Request a detailed protocol

Using bcftools (Li, 2011), we performed a locus-wise intersection of the following VCFs: 1000 Genomes (The 1000 Auton et al., 2015), IGDP (Jacobs et al., 2019), the high coverage Denisovan genome (Meyer et al., 2012), and the Altai and Vindija Neanderthal genomes (Prüfer et al., 2014). All VCFs corresponded to the GRCh37/hg19 reference sequence. Genotype matrices were constructed by parsing the output of bcftools query over 100 kbp windows, filtering out sites with sample allele frequency <5% or with more than 10% of genotypes missing, then excluding windows with fewer than 20 segregating sites. Each genotype matrix was then resized and sorted as described for simulations. When data were considered to be phased, as for the CEU/YRI populations, we also treated the Neanderthal genotypes as if they were phased according to REF/ALT columns in the VCF. While this is equivalent to random phasing, both high-coverage Neanderthal individuals are highly inbred, so this is unlikely to be problematic in practice.

CNN model architecture and training

Request a detailed protocol

We implemented the CNN model in Keras (Chollet, 2015), configured to use the Tensorflow backend (Abadi et al., 2015). To save disk space and memory, the input matrices were stored as eight bit integers rather than floating point numbers, and were not mean-centred or otherwise normalised prior to input into the network. We instead made the first layer of our network a batch normalisation layer, which simultaneously converts the input layer to floating point numbers and learns the best normalisation of the data for the network.

The CNN architecture (Figure 1) consists of $k$ convolution blocks each comprised of a batch normalisation layer followed by a 2D convolution layer with 2 × 2 stride, 16 filters of size 4 × 4, and leaky ReLU activation. The $k$ blocks are followed by a single fully-connected output node of size one, with sigmoid activation giving the probability Pr[AI]. We do not include pooling layers, as is common in a CNN architecture (e.g. Torada et al., 2019), and instead use a 2 × 2 stride size to reduce the output size of successive blocks (Springenberg et al., 2015). This is computationally cheaper and had no observable difference in network performance. We sought to maximise the depth of the network, but the size of the input matrix constrains the maximum number of blocks in the network due to successive halving of the dimensionality in each block. For $m = 256$ resizing bins, we used $k = 7$ blocks.

We partitioned 100,000 independent simulations for each of the three selection scenarios into training and validation sets (approximate 90%/10% split). The hyperparameters and network architecture were tuned on a smaller preliminary set of simulations that did not vary the selection coefficient or time of onset of selection, so we chose not to split the simulations into a third 'test’ set when evaluating the models trained on our final simulations. The model was trained for three epochs, with model weights updated after batches of 64, using the Adam optimiser and cross-entropy for the loss function. We evaluated model fit by inspecting loss and accuracy terms at the end of training (Appendix 2—table 1). Preliminary analyses indicated three epochs were sufficient for approximate convergence between training and validation metrics, but we did not observe divergence (likely indicating overfitting) even when training for additional epochs.

Comparison to other methods

Request a detailed protocol

We converted our simulated tree sequences to VolcanoFinder (Setter et al., 2020) input (a per-locus allele counts file, and a frequency-spectrum input file). One frequency-spectrum input file was created for each distinct demographic model, obtained by averaging over all neutral-scenario simulations. VolcanoFinder was run following examples in the manual, using 800 evenly-spaced genomic bins (-G 800), and taking the maximum value for the likelihood-ratio test statistic (LRT) as the summary statistic for that simulation. VolcanoFinder further requires a value for the divergence between the donor and recipient populations, which it can estimate by doing a grid search for the value which maximises the LRT. However, this is more computationally intensive than providing a value, so we obtained a value by grid search for a small sample of our simulations for each of Demographic Models A1 and B, and used the most frequently observed value (-D 0.001075 for Demographic Model A1 and A2, and -D 0.001465 for Demographic Model B).

Additional AI-related summary statistics were chosen based on Racimo et al., 2017b and calculated on the simulated tree sequences. The f_d statistic was implemented from the description in Martin et al., 2015 and the remaining statistics were implemented from their description in Racimo et al., 2017b. Summary statistics (including the VolcanoFinder LRT) were obtained for 600,000 simulations (100,000 for each of three simulation scenarios, for each of Demographic Models A1 and B). We calculated p-values by comparing each statistic to the null distribution that was obtained from the neutral simulation scenarios.

Calibration

Request a detailed protocol

For a well calibrated output, we expect proportion $x$ of the output probabilities with Pr[AI] $\sim x$ to be true positives. It has been noted elsewhere (Guo et al., 2017) that CNNs may produce improperly calibrated probabilities. However, even if the probabilities are calibrated with respect to the validation dataset (which has even class ratios), this is unlikely to hold for empirical data, as the relative ratios of AI versus not-AI windows in the genome are very skewed.

We tested three calibration methods: beta regression (Kull et al., 2017), isotonic regression (Chakravarti, 1989), and Platt, 1999 scaling. To calibrate our CNN output, we first resampled our training dataset to the desired class ratios. We then fit each calibrator to predict the true class in the resampled training dataset from the CNN prediction for the resampled training dataset. To assess the calibration procedure, we inspected reliability plots for our calibrated and uncalibrated predictions, as evaluated with a resampled validation dataset (Figure 4—figure supplement 1, Figure 4—figure supplement 2, Figure 4—figure supplement 3, Figure 4—figure supplement 4). We also checked if the sum of the residuals was normally distributed, following the approach of Turner et al., 2019. Both beta calibration and isotonic regression gave well-calibrated probabilities compared with uncalibrated model outputs, and for our predictions on empirical data we chose to apply beta calibration due to its relative simplicity.

Saliency maps

Request a detailed protocol

We computed average saliency maps, by averaging over a set of input-specific saliency maps that were calculated for a set of 300 simulated genotype matrices for each simulated scenario. The input-specific saliency maps were calculated using tf-keras-vis v0.5.5 (Kubota, 2020) configured to use ‘vanilla’ saliency maps. A sharper image was obtained by exchanging the CNN output layer’s sigmoid activation with linear activation, as recommended in the tf-keras-vis documentation. For the ‘vanilla’ saliency option, the image-specific class saliency is calculated by computing the gradient of a network’s output with respect to a single input. The exact details of how the saliency is calculated via propagation through a neural network can be found in Simonyan et al., 2014, who offer this interpretation: ‘[T]he magnitude of the derivative indicates which pixels need to be changed the least to affect the class score the most’.

Application of trained CNN to empirical datasets

Request a detailed protocol

We show Manhattan plots where each data point is a 100 kbp window that moves along the genome in steps of size 20 kbp. Gene annotations were extracted from the Ensembl release 87 GFF3 file (with GRCh37/hg19 coordinates), obtained via ensembl’s ftp server. We extracted the columns with source=‘ensembl_havana’ and type=‘gene’, and report the genes which intersected with the 30 top ranking CNN predictions or a 100 kbp flanking region. Adjacent regions were merged together prior to intersection, so that genes were reported only once.

Compute resources

Request a detailed protocol

All simulations and results reported here were obtained on an compute server with two Intel Xeon 6248 CPUs (80 cores total), 768 GB RAM, and five NVIDIA Tesla T4 GPUs. 300,000 SLiM simulations took approximately 80 weeks of CPU time for each of Demographic Model A1 and B. Each simulation executes independently, and is readily distributed across cores or compute nodes. This produced 450 GB of tree sequence files. The resized genotype matrices were compressed into a Zarr cache (Zarr Development Team, 2020) with size 2.8 GB, for faster loading. Training a single CNN on one GPU took approximately 15 min, or 60 CPU hours for an equivalent CPU-only training procedure. We did not attempt to optimise memory usage, and thus all data were loaded into memory, requiring approximately 120 GB RAM during training. Predicting AI for all genomic windows on an empirical dataset (22 single-chromosome BCF files) took 1 CPU hour. However, our prediction pipeline uses multiprocessing and efficiently scales to 80 cores.

Code availability

Request a detailed protocol

The source code for performing simulations, training and evaluating a CNN, and applying a CNN to empirical VCF data, were developed in a new Python application called genomatnn, available at Gower, 2021.

Appendix 1

Appendix 1—table 1

Top ranking gene candidates corresponding to Neanderthal AI in Europeans.

We show genes which overlap, or are within 100 kbp of, the 30 highest ranked 100 kbp intervals. Adjacent intervals have been merged. The CNN was trained using only AI simulations with selected mutation having allele frequency >5%, and subsequently calibrated with resampled neutral:sweep:AI ratios of 1:0.1:0.02 .c

Chrom	Start	End	Genes
1	39420001	39520000	RRAGC; MYCBP; GJA9; RHBDL2; AKIRIN1; NDUFS5; MACF1
2	159880001	160280000	TANC1; WDSUB1; BAZ2B
2	180060001	180160000	SESTD1
2	227800001	227900000	RHBDD1; COL4A4
2	238820001	238960000	LRRFIP1; RBM44; RAMP1; UBE2F; SCLY; ESPNL; KLHL30
3	114500001	114600000	ZBTB20
5	57960001	58060000	RAB3C
6	28160001	28380000	ZSCAN16-AS1; ZSCAN16; ZKSCAN8; ZSCAN9; ZKSCAN4; NKAPL; PGBD1; ZSCAN31; ZKSCAN3; ZSCAN12; ZSCAN23; GPX6
8	17060001	17160000	MICU3; ZDHHC2; CNOT7; VPS37A; MTMR7
8	91840001	91940000	TMEM64; NECAB1; TMEM55A
9	16700001	16860000	BNC2
10	11800001	11900000	ECHDC3; PROSER2; UPF2
11	37740001	37840000
19	20260001	20360000	ZNF90; ZNF486
19	33580001	33700000	RHPN2; GPATCH1; WDR88; LRP3; SLC7A10
20	14340001	14440000	MACROD2; FLRT3

Appendix 1—table 2

Top ranking gene candidates corresponding to Denisovan AI in Melanesians.

We show genes which overlap, or are within 100 kbp of, the 30 highest ranked 100 kbp intervals. Adjacent intervals have been merged. The CNN was trained using only AI simulations with selected mutation having allele frequency >5%, and subsequently calibrated with resampled neutral:sweep:AI ratios of 1:0.1:0.02.

Chrom	Start	End	Genes
1	2880001	2980000	ACTRT2; LINC00982; PRDM16
1	220080001	220180000	SLC30A10; EPRS; BPNT1; IARS2
2	221040001	221140000
3	15400001	15500000	SH3BP5; METTL6; EAF1; COLQ
4	41960001	42100000	TMEM33; DCAF4L1; SLC30A9; BEND4
5	135440001	135540000	TGFBI; SMAD5-AS1; SMAD5; TRPC7
6	81980001	82120000	FAM46A
7	121160001	121260000
9	95500001	95600000	IPPK; BICD2; ZNF484
10	59660001	59760000
12	80780001	80880000	OTOGL; PTPRQ
12	84620001	84740000
14	57620001	57760000	EXOC5; AP5M1; NAA30
17	29480001	29720000	NF1; OMG; EVI2B; EVI2A; RAB11FIP4
18	38180001	38320000
20	54340001	54440000

Appendix 2

Appendix 2—table 1

Loss and accuracy for CNNs after training for three epochs, as reported by Keras/Tensorflow, for the training and validation datasets.

Binary cross-entropy was used for the loss function.

Demographic model	Hyperparameters	Training		Validation
		Loss	Accuracy	Loss	Accuracy
A1	AF>0.05	0.1592	0.9458	0.1618	0.9468
A1	AF>0.25	0.1224	0.9585	0.1265	0.9578
A1	AF>0.25; unphased	0.1347	0.9537	0.1368	0.9530
B	AF>0.05; unphased	0.3415	0.8439	0.3441	0.8439
B	AF>0.25; unphased	0.3546	0.8372	0.3583	0.8376

Appendix 3

Appendix 3—table 1

Parameter values used for simulating Demographic Model A1.

A Demes-format YAML file for each demographic model is available from the genomatnn git repository.

Parameter	Description	Value	Units	Source
$N_{A n c}$	ancestral pop. size	18500		Kuhlwilm et al., 2016
$N_{N e a}$	Neanderthal pop. size	3400		Kuhlwilm et al., 2016
$N_{Y R I}$	YRI pop. size	27600		Kuhlwilm et al., 2016
$N_{C E U 0}$	CEU bottleneck pop. size	1080		Ragsdale and Gravel, 2019
$N_{C E U 1}$	CEU growth-start pop. size	1450		Ragsdale and Gravel, 2019
$N_{C E U 2}$	CEU current pop. size	13377
$r_{C E U}$	CEU growth rate	0.00202		Ragsdale and Gravel, 2019
$T_{C E U 2}$	CEU time at growth start	31.9	kya	Ragsdale and Gravel, 2019
$T_{0}$	Nea/other split time	550	kya	Prüfer et al., 2017
$T_{1}$	CEU/YRI split time	65.7	kya	Ragsdale and Gravel, 2019
$T_{2}$	time of Nea → CEU gene flow	55	kya	Prüfer et al., 2017
$g$	generation time	29	years	Prüfer et al., 2017
$α$	Nea → CEU admixture proportion	2.25		Prüfer et al., 2017
$T_{Altai}$	sampling time	115	kya	Prüfer et al., 2017
$T_{Vindija}$	sampling time	55	kya	Prüfer et al., 2017
$n_{Nean}$	sample size	2	diploid individuals
$n_{Afr}$	sample size	108	diploid individuals
$n_{Eur}$	sample size	99	diploid individuals
$s$	selection coefficient	$10^{Unif(-4,-1)}$
$T_{sel1}$	selection onset (sweep)	Unif(1, $T_{1}$ )	kya
$T_{mut1}$	mutation (sweep)	Unif( $T_{sel1}$ , $T_{1}$ )	kya
$T_{sel2}$	selection onset (AI)	Unif(1, $T_{2}$ )	kya
$T_{mut2}$	mutation (AI)	Unif( $T_{2}$ , $T_{0}$ )	kya

Appendix 4

Appendix 4—figure 1

Download asset Open asset

Haplotype plot for the candidate region chr1:104500001–104600000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 2

Download asset Open asset

Haplotype plot for the candidate region chr2:109360001–109460000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 3

Download asset Open asset

Haplotype plot for the candidate region chr2:160160001–160280000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 4

Download asset Open asset

Haplotype plot for the candidate region chr3:114480001–114620000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 5

Download asset Open asset

Haplotype plot for the candidate region chr4:54240001–54340000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 6

Download asset Open asset

Haplotype plot for the candidate region chr5:39220001–39320000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 7

Download asset Open asset

Haplotype plot for the candidate region chr6:28180001–28320000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 8

Download asset Open asset

Haplotype plot for the candidate region chr8:143440001–143560000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 9

Download asset Open asset

Haplotype plot for the candidate region chr9:16700001–16820000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 10

Download asset Open asset

Haplotype plot for the candidate region chr12:85780001–85880000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 11

Download asset Open asset

Haplotype plot for the candidate region chr19:20220001–20380000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 12

Download asset Open asset

Haplotype plot for the candidate region chr19:33580001–33740000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 13

Download asset Open asset

Haplotype plot for the candidate region chr20:62100001–62280000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 4—figure 14

Download asset Open asset

Haplotype plot for the candidate region chr21:25840001–25940000 in the Neanderthal-into-European AI scan.

Bright yellow indicates minor allele, dark blue indicates major allele. Haplotypes within populations are sorted left-to-right by similarity to Neanderthals.

Appendix 5

Appendix 5—figure 1

Download asset Open asset

Genotype plot for the candidate region chr2:129960001–130060000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 2

Download asset Open asset

Genotype plot for the candidate region chr3:3740001–3840000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 3

Download asset Open asset

Genotype plot for the candidate region chr4:41980001–42080000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 4

Download asset Open asset

Genotype plot for the candidate region chr5:420001–520000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 5

Download asset Open asset

Genotype plot for the candidate region chr6:74640001–74740000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 6

Download asset Open asset

Genotype plot for the candidate region chr6:81960001–82060000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 7

Download asset Open asset

Genotype plot for the candidate region chr6:137920001–138120000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 8

Download asset Open asset

Genotype plot for the candidate region chr7:25100001–25200000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 9

Download asset Open asset

Genotype plot for the candidate region chr7:38020001–38120000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 10

Download asset Open asset

Genotype plot for the candidate region chr7:121160001–121260000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 11

Download asset Open asset

Genotype plot for the candidate region chr8:3040001–3140000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 12

Download asset Open asset

Genotype plot for the candidate region chr12:84640001–84740000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 13

Download asset Open asset

Genotype plot for the candidate region chr12:108240001–108340000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 14

Download asset Open asset

Genotype plot for the candidate region chr12:114020001–114280000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 15

Download asset Open asset

Genotype plot for the candidate region chr14:61860001–61960000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 16

Download asset Open asset

Genotype plot for the candidate region chr14:63120001–63220000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 17

Download asset Open asset

Genotype plot for the candidate region chr14:96700001–96820000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 18

Download asset Open asset

Genotype plot for the candidate region chr15:55260001–55400000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 19

Download asset Open asset

Genotype plot for the candidate region chr16:62600001–62700000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 20

Download asset Open asset

Genotype plot for the candidate region chr16:78360001–78460000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 21

Download asset Open asset

Genotype plot for the candidate region chr18:22060001–22160000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Appendix 5—figure 22

Download asset Open asset

Genotype plot for the candidate region chr22:19040001–19140000 in the Denisovan-into-Melanesian AI scan.

Dark blue = homozygote major allele, light blue = heterozygote, yellow = homozygote minor allele. Genotypes within populations are sorted left-to-right by similarity to the Denisovan.

Data availability

Source code is available from https://github.com/grahamgower/genomatnn/.

References

Preprint
1. Abadi M
2. Agarwal A
3. Barham P
4. Brevdo E
5. Chen Z
6. Citro C
7. Corrado GS
8. Davis A
9. Dean J
10. Devin M
11. Ghemawat S
12. Goodfellow I
13. Harp A
14. Irving G
15. Isard M
16. Jia Y
17. Jozefowicz R
18. Kaiser L
19. Kudlur M
20. Levenberg J
21. Mané D
22. Monga R
23. Moore S
24. Murray D
25. Olah C
26. Schuster M
27. Shlens J
28. Steiner B
29. Sutskever I
30. Talwar K
31. Tucker P
32. Vanhoucke V
33. Vasudevan V
34. Viégas F
35. Vinyals O
36. Warden P
37. Wattenberg M
38. Wicke M
39. Yu Y
40. Zheng X
(2015) TensorFlow: large-scale machine learning on heterogeneous systems
arXiv.

https://arxiv.org/abs/1603.04467
- Google Scholar
1. Adrion JR
2. Cole CB
3. Dukler N
4. Galloway JG
5. Gladstein AL
6. Gower G
7. Kyriazis CC
8. Ragsdale AP
9. Tsambos G
10. Baumdicker F
11. Carlson J
12. Cartwright RA
13. Durvasula A
14. Gronau I
15. Kim BY
16. McKenzie P
17. Messer PW
18. Noskova E
19. Ortega-Del Vecchyo D
20. Racimo F
21. Struck TJ
22. Gravel S
23. Gutenkunst RN
24. Lohmueller KE
25. Ralph PL
26. Schrider DR
27. Siepel A
28. Kelleher J
29. Kern AD
(2020a) A community-maintained standard library of population genetic models
eLife 9:e54967.

https://doi.org/10.7554/eLife.54967
- PubMed
- Google Scholar
(2020b) Predicting the landscape of recombination using deep learning
Molecular Biology and Evolution 37:1790–1808.

https://doi.org/10.1093/molbev/msaa038
- PubMed
- Google Scholar
Book
1. Aggarwal CC
(2018) Neural Networks and Deep Learning
Springer.

https://doi.org/10.1007/978-3-319-94463-0
- Google Scholar
Book
1. Alaa AM
2. van der Schaar M
(2019)
Demystifying black-box models with symbolic metamodels

In: Wallach H, Larochelle H, Beygelzimer A, Alché-Buc F. d, Fox E, Garnett R, editors. Advances in Neural Information Processing Systems 32. Curran Associates, Inc. pp. 11304–11314.
- Google Scholar
(2015) A global reference for human genetic variation
Nature 526:68–74.

https://doi.org/10.1038/nature15393
- PubMed
- Google Scholar
1. Bai Z
2. Stamova B
3. Xu H
4. Ander BP
5. Wang J
6. Jickling GC
7. Zhan X
8. Liu D
9. Han G
10. Jin LW
11. DeCarli C
12. Lei H
13. Sharp FR
(2014) Distinctive RNA expression profiles in blood associated with Alzheimer disease after accounting for white matter hyperintensities
Alzheimer Disease and Associated Disorders 28:226–233.

https://doi.org/10.1097/WAD.0000000000000022
- PubMed
- Google Scholar
(2021) Chromosome-scale inference of hybrid speciation and admixture with convolutional neural networks
Molecular Ecology Resources 8:13355.

https://doi.org/10.1111/1755-0998.13355
- Google Scholar
(2018) An atlas of genetic associations in UK Biobank
Nature Genetics 50:1593–1599.

https://doi.org/10.1038/s41588-018-0248-z
- PubMed
- Google Scholar
Preprint
(2020) The MCC-F1 curve: a performance evaluation technique for binary classification
arXiv.

https://arxiv.org/abs/2006.11278
- Google Scholar
1. Chakravarti N
(1989) Isotonic median regression: a linear programming approach
Mathematics of Operations Research 14:303–308.

https://doi.org/10.1287/moor.14.2.303
- Google Scholar
Book
1. Chan J
2. Perrone V
3. Spence J
4. Jenkins P
5. Mathieson S
6. Song Y
(2018)
A likelihood-free inference framework for population genetic data using exchangeable neural networks

In: Bengio S, Wallach H, Larochelle H, Grauman K, Cesa-Bianchi N, Garnett R, editors. Advances in Neural Information Processing Systems 31. Curran Associates, Inc. pp. 8594–8605.
- Google Scholar
1. Chicco D
(2017) Ten quick tips for machine learning in computational biology
BioData Mining 10:35.

https://doi.org/10.1186/s13040-017-0155-3
- PubMed
- Google Scholar
Software
1. Chollet F
(2015) Keras
Keras.

https://keras.io
1. Cosma MP
2. Pepe S
3. Annunziata I
4. Newbold RF
5. Grompe M
6. Parenti G
7. Ballabio A
(2003) The multiple sulfatase deficiency gene encodes an essential and limiting factor for the activity of sulfatases
Cell 113:445–456.

https://doi.org/10.1016/s0092-8674(03)00348-9
- PubMed
- Google Scholar
1. De Rocco D
2. Cerqua C
3. Goffrini P
4. Russo G
5. Pastore A
6. Meloni F
7. Nicchia E
8. Moraes CT
9. Pecci A
10. Salviati L
11. Savoia A
(2014) Mutations of cytochrome c identified in patients with thrombocytopenia THC4 affect both apoptosis and cellular bioenergetics
Biochimica Et Biophysica Acta (BBA) - Molecular Basis of Disease 1842:269–274.

https://doi.org/10.1016/j.bbadis.2013.12.002
- PubMed
- Google Scholar
1. Durvasula A
2. Sankararaman S
(2019) A statistical model for reference-free inference of archaic local ancestry
PLOS Genetics 15:e1008175.

https://doi.org/10.1371/journal.pgen.1008175
- PubMed
- Google Scholar
(2017) Evolution of novel mimicry rings facilitated by adaptive introgression in tropical butterflies
Molecular Ecology 26:5160–5172.

https://doi.org/10.1111/mec.14277
- PubMed
- Google Scholar
(2019) The unreasonable effectiveness of convolutional neural networks in population genetic inference
Molecular Biology and Evolution 36:220–238.

https://doi.org/10.1093/molbev/msy224
- PubMed
- Google Scholar
1. Frazer KA
2. Ballinger DG
3. Cox DR
4. Hinds DA
5. Stuve LL
6. Gibbs RA
7. Belmont JW
8. Boudreau A
9. Hardenbol P
10. Leal SM
11. Pasternak S
12. Wheeler DA
13. Willis TD
14. Yu F
15. Yang H
16. Zeng C
17. Gao Y
18. Hu H
19. Hu W
20. Li C
21. Lin W
22. Liu S
23. Pan H
24. Tang X
25. Wang J
26. Wang W
27. Yu J
28. Zhang B
29. Zhang Q
30. Zhao H
31. Zhao H
32. Zhou J
33. Gabriel SB
34. Barry R
35. Blumenstiel B
36. Camargo A
37. Defelice M
38. Faggart M
39. Goyette M
40. Gupta S
41. Moore J
42. Nguyen H
43. Onofrio RC
44. Parkin M
45. Roy J
46. Stahl E
47. Winchester E
48. Ziaugra L
49. Altshuler D
50. Shen Y
51. Yao Z
52. Huang W
53. Chu X
54. He Y
55. Jin L
56. Liu Y
57. Shen Y
58. Sun W
59. Wang H
60. Wang Y
61. Wang Y
62. Xiong X
63. Xu L
64. Waye MM
65. Tsui SK
66. Xue H
67. Wong JT
68. Galver LM
69. Fan JB
70. Gunderson K
71. Murray SS
72. Oliphant AR
73. Chee MS
74. Montpetit A
75. Chagnon F
76. Ferretti V
77. Leboeuf M
78. Olivier JF
79. Phillips MS
80. Roumy S
81. Sallée C
82. Verner A
83. Hudson TJ
84. Kwok PY
85. Cai D
86. Koboldt DC
87. Miller RD
88. Pawlikowska L
89. Taillon-Miller P
90. Xiao M
91. Tsui LC
92. Mak W
93. Song YQ
94. Tam PK
95. Nakamura Y
96. Kawaguchi T
97. Kitamoto T
98. Morizono T
99. Nagashima A
100. Ohnishi Y
101. Sekine A
102. Tanaka T
103. Tsunoda T
104. Deloukas P
105. Bird CP
106. Delgado M
107. Dermitzakis ET
108. Gwilliam R
109. Hunt S
110. Morrison J
111. Powell D
112. Stranger BE
113. Whittaker P
114. Bentley DR
115. Daly MJ
116. de Bakker PI
117. Barrett J
118. Chretien YR
119. Maller J
120. McCarroll S
121. Patterson N
122. Pe'er I
123. Price A
124. Purcell S
125. Richter DJ
126. Sabeti P
127. Saxena R
128. Schaffner SF
129. Sham PC
130. Varilly P
131. Altshuler D
132. Stein LD
133. Krishnan L
134. Smith AV
135. Tello-Ruiz MK
136. Thorisson GA
137. Chakravarti A
138. Chen PE
139. Cutler DJ
140. Kashuk CS
141. Lin S
142. Abecasis GR
143. Guan W
144. Li Y
145. Munro HM
146. Qin ZS
147. Thomas DJ
148. McVean G
149. Auton A
150. Bottolo L
151. Cardin N
152. Eyheramendy S
153. Freeman C
154. Marchini J
155. Myers S
156. Spencer C
157. Stephens M
158. Donnelly P
159. Cardon LR
160. Clarke G
161. Evans DM
162. Morris AP
163. Weir BS
164. Tsunoda T
165. Mullikin JC
166. Sherry ST
167. Feolo M
168. Skol A
169. Zhang H
170. Zeng C
171. Zhao H
172. Matsuda I
173. Fukushima Y
174. Macer DR
175. Suda E
176. Rotimi CN
177. Adebamowo CA
178. Ajayi I
179. Aniagwu T
180. Marshall PA
181. Nkwodimmah C
182. Royal CD
183. Leppert MF
184. Dixon M
185. Peiffer A
186. Qiu R
187. Kent A
188. Kato K
189. Niikawa N
190. Adewole IF
191. Knoppers BM
192. Foster MW
193. Clayton EW
194. Watkin J
195. Gibbs RA
196. Belmont JW
197. Muzny D
198. Nazareth L
199. Sodergren E
200. Weinstock GM
201. Wheeler DA
202. Yakub I
203. Gabriel SB
204. Onofrio RC
205. Richter DJ
206. Ziaugra L
207. Birren BW
208. Daly MJ
209. Altshuler D
210. Wilson RK
211. Fulton LL
212. Rogers J
213. Burton J
214. Carter NP
215. Clee CM
216. Griffiths M
217. Jones MC
218. McLay K
219. Plumb RW
220. Ross MT
221. Sims SK
222. Willey DL
223. Chen Z
224. Han H
225. Kang L
226. Godbout M
227. Wallenburg JC
228. L'Archevêque P
229. Bellemare G
230. Saeki K
231. Wang H
232. An D
233. Fu H
234. Li Q
235. Wang Z
236. Wang R
237. Holden AL
238. Brooks LD
239. McEwen JE
240. Guyer MS
241. Wang VO
242. Peterson JL
243. Shi M
244. Spiegel J
245. Sung LM
246. Zacharia LF
247. Collins FS
248. Kennedy K
249. Jamieson R
250. Stewart J
251. International HapMap Consortium
(2007) A second generation human haplotype map of over 3.1 million SNPs
Nature 449:851–861.

https://doi.org/10.1038/nature06258
- PubMed
- Google Scholar
(2016) Archaic Hominin Admixture Facilitated Adaptation to Out-of-Africa Environments
Current Biology : CB 26:3375–3382.

https://doi.org/10.1016/j.cub.2016.10.041
- PubMed
- Google Scholar
Software
1. Gower G
(2021) Predicts adaptive introgression using a CNN trained on genotype matrices, version 7a51abd
GitHub.

https://github.com/grahamgower/genomatnn
1. Green RE
2. Krause J
3. Briggs AW
4. Maricic T
5. Stenzel U
6. Kircher M
7. Patterson N
8. Li H
9. Zhai W
10. Fritz MH
11. Hansen NF
12. Durand EY
13. Malaspinas AS
14. Jensen JD
15. Marques-Bonet T
16. Alkan C
17. Prüfer K
18. Meyer M
19. Burbano HA
20. Good JM
21. Schultz R
22. Aximu-Petri A
23. Butthof A
24. Höber B
25. Höffner B
26. Siegemund M
27. Weihmann A
28. Nusbaum C
29. Lander ES
30. Russ C
31. Novod N
32. Affourtit J
33. Egholm M
34. Verna C
35. Rudan P
36. Brajkovic D
37. Kucan Ž
38. Gušic I
39. Doronichev VB
40. Golovanova LV
41. Lalueza-Fox C
42. de la Rasilla M
43. Fortea J
44. Rosas A
45. Schmitz RW
46. Johnson PLF
47. Eichler EE
48. Falush D
49. Birney E
50. Mullikin JC
51. Slatkin M
52. Nielsen R
53. Kelso J
54. Lachmann M
55. Reich D
56. Pääbo S
(2010) A draft sequence of the Neandertal genome
Science 328:710–722.

https://doi.org/10.1126/science.1188021
- PubMed
- Google Scholar
Preprint
1. Guo C
2. Pleiss G
3. Sun Y
4. Weinberger KQ
(2017) On calibration of modern neural networks
arXiv.

https://arxiv.org/abs/1706.04599
- Google Scholar
(2019) Tree-sequence recording in SLiM opens new horizons for forward-time simulation of whole genomes
Molecular Ecology Resources 19:552–566.

https://doi.org/10.1111/1755-0998.12968
- PubMed
- Google Scholar
1. Haller BC
2. Messer PW
(2019) SLiM 3: Forward Genetic Simulations Beyond the Wright-Fisher Model
Molecular Biology and Evolution 36:632–637.

https://doi.org/10.1093/molbev/msy228
- PubMed
- Google Scholar
1. Harris K
2. Nielsen R
(2016) The Genetic Cost of Neanderthal Introgression
Genetics 203:881–891.

https://doi.org/10.1534/genetics.116.186890
- PubMed
- Google Scholar
1. Hawks J
2. Cochran G
(2006)
Dynamics of adaptive introgression from archaic to modern humans

PaleoAnthropology 2006:101–115.
- Google Scholar
(2016) The genetics of extreme microgeographic adaptation: an integrated approach identifies a major gene underlying leaf trichome divergence in Yellowstone Mimulus guttatus
Molecular Ecology 25:5647–5662.

https://doi.org/10.1111/mec.13753
- PubMed
- Google Scholar
1. Hermisson J
2. Pennings PS
(2005) Soft sweeps: molecular population genetics of adaptation from standing genetic variation
Genetics 169:2335–2352.

https://doi.org/10.1534/genetics.104.036947
- PubMed
- Google Scholar
(2012) NLRR1 enhances EGF-mediated MYCN induction in neuroblastoma and accelerates tumor growth in vivo
Cancer Research 72:4587–4596.

https://doi.org/10.1158/0008-5472.CAN-12-0943
- PubMed
- Google Scholar
(2020) Mapping gene flow between ancient hominins through demography-aware inference of the ancestral recombination graph
PLOS Genetics 16:e1008895.

https://doi.org/10.1371/journal.pgen.1008895
- PubMed
- Google Scholar
1. Huerta-Sánchez E
2. Jin X
3. Asan
4. Bianba Z
5. Peter BM
6. Vinckenbosch N
7. Liang Y
8. Yi X
9. He M
10. Somel M
11. Ni P
12. Wang B
13. Ou X
14. Huasang
15. Luosang J
16. Cuo ZX
17. Li K
18. Gao G
19. Yin Y
20. Wang W
21. Zhang X
22. Xu X
23. Yang H
24. Li Y
25. Wang J
26. Wang J
27. Nielsen R
(2014) Altitude adaptation in Tibetans caused by introgression of Denisovan-like DNA
Nature 512:194.

https://doi.org/10.1038/nature13408
- PubMed
- Google Scholar
(2021) Distinguishing between recent balancing selection and incomplete sweep using deep neural networks
Molecular Ecology Resources 1:13379.

https://doi.org/10.1111/1755-0998.13379
- Google Scholar
1. Jacobs GS
2. Hudjashov G
3. Saag L
4. Kusuma P
5. Darusallam CC
6. Lawson DJ
7. Mondal M
8. Pagani L
9. Ricaut FX
10. Stoneking M
11. Metspalu M
12. Sudoyo H
13. Lansing JS
14. Cox MP
(2019) Multiple deeply divergent denisovan ancestries in papuans
Cell 177:1010–1021.

https://doi.org/10.1016/j.cell.2019.02.035
- PubMed
- Google Scholar
1. Jones MR
2. Mills LS
3. Alves PC
4. Callahan CM
5. Alves JM
6. Lafferty DJR
7. Jiggins FM
8. Jensen JD
9. Melo-Ferreira J
10. Good JM
(2018) Adaptive introgression underlies polymorphic seasonal camouflage in snowshoe hares
Science 360:1355–1358.

https://doi.org/10.1126/science.aar5273
- PubMed
- Google Scholar
(2004) Hash4, a novel human achaete-scute homologue found in fetal skin
Genomics 84:859–866.

https://doi.org/10.1016/j.ygeno.2004.07.004
- PubMed
- Google Scholar
(2016) The Strength of Selection against Neanderthal Introgression
PLOS Genetics 12:e1006340.

https://doi.org/10.1371/journal.pgen.1006340
- PubMed
- Google Scholar
(2016) Efficient Coalescent Simulation and Genealogical Analysis for Large Sample Sizes
PLOS Computational Biology 12:e1004842.

https://doi.org/10.1371/journal.pcbi.1004842
- PubMed
- Google Scholar
(2018) Efficient pedigree recording for fast population genetics simulation
PLOS Computational Biology 14:e1006581.

https://doi.org/10.1371/journal.pcbi.1006581
- PubMed
- Google Scholar
(2017) Inference of the Distribution of Selection Coefficients for New Nonsynonymous Mutations Using Large Samples
Genetics 206:345–361.

https://doi.org/10.1534/genetics.116.197145
- PubMed
- Google Scholar
(2018) Deleterious variation shapes the genomic landscape of introgression
PLOS Genetics 14:e1007741.

https://doi.org/10.1371/journal.pgen.1007741
- PubMed
- Google Scholar
1. Kiper POS
2. Saito H
3. Gori F
4. Unger S
5. Hesse E
6. Yamana K
7. Kiviranta R
8. Solban N
9. Liu J
10. Brommage R
11. Boduroglu K
12. Bonafé L
13. Campos-Xavier B
14. Dikoglu E
15. Eastell R
16. Gossiel F
17. Harshman K
18. Nishimura G
19. Girisha KM
20. Stevenson BJ
21. Takita H
22. Rivolta C
23. Superti-Furga A
24. Baron R
(2016) Cortical-Bone fragility--insights from sFRP4 deficiency in Pyle's Disease
New England Journal of Medicine 374:2553–2562.

https://doi.org/10.1056/NEJMoa1509342
- PubMed
- Google Scholar
(2014) ANGSD: Analysis of Next Generation Sequencing Data
BMC Bioinformatics 15:356.

https://doi.org/10.1186/s12859-014-0356-4
- PubMed
- Google Scholar
Conference
(2012)
Imagenet classification with deep convolutional neural networks

Advances in Neural Information Processing Systems. pp. 84–90.
- Google Scholar
1. Kubo M
2. Hata J
3. Ninomiya T
4. Matsuda K
5. Yonemoto K
6. Nakano T
7. Matsushita T
8. Yamazaki K
9. Ohnishi Y
10. Saito S
11. Kitazono T
12. Ibayashi S
13. Sueishi K
14. Iida M
15. Nakamura Y
16. Kiyohara Y
(2007) A nonsynonymous SNP in PRKCH (protein kinase C eta) increases the risk of cerebral infarction
Nature Genetics 39:212–217.

https://doi.org/10.1038/ng1945
- PubMed
- Google Scholar
Software
1. Kubota Y
(2020) tf-keras-vis
tf-keras-vis.

https://github.com/keisen/tf-keras-vis
1. Kuhlwilm M
2. Gronau I
3. Hubisz MJ
4. de Filippo C
5. Prado-Martinez J
6. Kircher M
7. Fu Q
8. Burbano HA
9. Lalueza-Fox C
10. de la Rasilla M
11. Rosas A
12. Rudan P
13. Brajkovic D
14. Kucan Ž
15. Gušic I
16. Marques-Bonet T
17. Andrés AM
18. Viola B
19. Pääbo S
20. Meyer M
21. Siepel A
22. Castellano S
(2016) Ancient gene flow from early modern humans into Eastern Neanderthals
Nature 530:429–433.

https://doi.org/10.1038/nature16544
- PubMed
- Google Scholar
Conference
1. Kull M
2. Filho TS
3. Flach P
(2017)
Beta calibration: a well-founded and easily implemented improvement on logistic calibration for binary classifiers

Proceedings of the 20th International Conference on Artificial Intelligence and Statistics. pp. 623–631.
- Google Scholar
Book
1. LeCun Y
2. Bengio Y
(1995)
Convolutional networks for images, speech, and time series

In: Arbib M. A, editors. The Handbook of Brain Theory and Neural Networks. MIT Press. pp. 255–258.
- Google Scholar
1. Li H
(2011) A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data
Bioinformatics 27:2987–2993.

https://doi.org/10.1093/bioinformatics/btr509
- PubMed
- Google Scholar
1. Mahdi T
2. Hänzelmann S
3. Salehi A
4. Muhammed SJ
5. Reinbothe TM
6. Tang Y
7. Axelsson AS
8. Zhou Y
9. Jing X
10. Almgren P
11. Krus U
12. Taneera J
13. Blom AM
14. Lyssenko V
15. Esguerra JL
16. Hansson O
17. Eliasson L
18. Derry J
19. Zhang E
20. Wollheim CB
21. Groop L
22. Renström E
23. Rosengren AH
(2012) Secreted frizzled-related protein 4 reduces insulin secretion and is overexpressed in type 2 diabetes
Cell Metabolism 16:625–633.

https://doi.org/10.1016/j.cmet.2012.10.009
- PubMed
- Google Scholar
1. Malaspinas AS
2. Westaway MC
3. Muller C
4. Sousa VC
5. Lao O
6. Alves I
7. Bergström A
8. Athanasiadis G
9. Cheng JY
10. Crawford JE
11. Heupink TH
12. Macholdt E
13. Peischl S
14. Rasmussen S
15. Schiffels S
16. Subramanian S
17. Wright JL
18. Albrechtsen A
19. Barbieri C
20. Dupanloup I
21. Eriksson A
22. Margaryan A
23. Moltke I
24. Pugach I
25. Korneliussen TS
26. Levkivskyi IP
27. Moreno-Mayar JV
28. Ni S
29. Racimo F
30. Sikora M
31. Xue Y
32. Aghakhanian FA
33. Brucato N
34. Brunak S
35. Campos PF
36. Clark W
37. Ellingvåg S
38. Fourmile G
39. Gerbault P
40. Injie D
41. Koki G
42. Leavesley M
43. Logan B
44. Lynch A
45. Matisoo-Smith EA
46. McAllister PJ
47. Mentzer AJ
48. Metspalu M
49. Migliano AB
50. Murgha L
51. Phipps ME
52. Pomat W
53. Reynolds D
54. Ricaut FX
55. Siba P
56. Thomas MG
57. Wales T
58. Wall CM
59. Oppenheimer SJ
60. Tyler-Smith C
61. Durbin R
62. Dortch J
63. Manica A
64. Schierup MH
65. Foley RA
66. Lahr MM
67. Bowern C
68. Wall JD
69. Mailund T
70. Stoneking M
71. Nielsen R
72. Sandhu MS
73. Excoffier L
74. Lambert DM
75. Willerslev E
(2016) A genomic history of Aboriginal Australia
Nature 538:207–214.

https://doi.org/10.1038/nature18299
- PubMed
- Google Scholar
(2015) Evaluating the use of ABBA-BABA statistics to locate introgressed loci
Molecular Biology and Evolution 32:244–257.

https://doi.org/10.1093/molbev/msu269
- PubMed
- Google Scholar
1. Matthews BW
(1975) Comparison of the predicted and observed secondary structure of T4 phage lysozyme
Biochimica Et Biophysica Acta (BBA) - Protein Structure 405:442–451.

https://doi.org/10.1016/0005-2795(75)90109-9
- PubMed
- Google Scholar
1. Ménasché G
2. Pastural E
3. Feldmann J
4. Certain S
5. Ersoy F
6. Dupuis S
7. Wulffraat N
8. Bianchi D
9. Fischer A
10. Le Deist F
11. de Saint Basile G
(2000) Mutations in RAB27A cause Griscelli syndrome associated with haemophagocytic syndrome
Nature Genetics 25:173–176.

https://doi.org/10.1038/76024
- PubMed
- Google Scholar
1. Meyer M
2. Kircher M
3. Gansauge MT
4. Li H
5. Racimo F
6. Mallick S
7. Schraiber JG
8. Jay F
9. Prüfer K
10. de Filippo C
11. Sudmant PH
12. Alkan C
13. Fu Q
14. Do R
15. Rohland N
16. Tandon A
17. Siebauer M
18. Green RE
19. Bryc K
20. Briggs AW
21. Stenzel U
22. Dabney J
23. Shendure J
24. Kitzman J
25. Hammer MF
26. Shunkov MV
27. Derevianko AP
28. Patterson N
29. Andrés AM
30. Eichler EE
31. Slatkin M
32. Reich D
33. Kelso J
34. Pääbo S
(2012) A high-coverage genome sequence from an archaic Denisovan individual
Science 338:222–226.

https://doi.org/10.1126/science.1224344
- PubMed
- Google Scholar
1. Miao B
2. Wang Z
3. Li Y
(2017) Genomic analysis reveals hypoxia adaptation in the tibetan mastiff by introgression of the gray wolf from the tibetan plateau
Molecular Biology and Evolution 34:734–743.

https://doi.org/10.1093/molbev/msw274
- PubMed
- Google Scholar
(2019) Approximate Bayesian computation with deep learning supports a third archaic introgression in Asia and Oceania
Nature Communications 10:246.

https://doi.org/10.1038/s41467-018-08089-7
- PubMed
- Google Scholar
1. Morison IM
2. Cramer Bordé EM
3. Cheesman EJ
4. Cheong PL
5. Holyoake AJ
6. Fichelson S
7. Weeks RJ
8. Lo A
9. Davies SM
10. Wilbanks SM
11. Fagerlund RD
12. Ludgate MW
13. da Silva Tatley FM
14. Coker MS
15. Bockett NA
16. Hughes G
17. Pippig DA
18. Smith MP
19. Capron C
20. Ledgerwood EC
(2008) A mutation of human cytochrome c enhances the intrinsic apoptotic pathway but causes only thrombocytopenia
Nature Genetics 40:387–389.

https://doi.org/10.1038/ng.103
- PubMed
- Google Scholar
1. Norris LC
2. Main BJ
3. Lee Y
4. Collier TC
5. Fofana A
6. Cornel AJ
7. Lanzaro GC
(2015) Adaptive introgression in an African malaria mosquito coincident with the increased usage of insecticide-treated bed nets
PNAS 112:815–820.

https://doi.org/10.1073/pnas.1418892112
- PubMed
- Google Scholar
(2012) Adaptive introgression across species boundaries in Heliconius butterflies
PLOS Genetics 8:e1002752.

https://doi.org/10.1371/journal.pgen.1002752
- PubMed
- Google Scholar
Conference
1. Platt JC
(1999)
Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods

Advances in Large Margin Classifiers. pp. 61–74.
- Google Scholar
1. Prüfer K
2. Racimo F
3. Patterson N
4. Jay F
5. Sankararaman S
6. Sawyer S
7. Heinze A
8. Renaud G
9. Sudmant PH
10. de Filippo C
11. Li H
12. Mallick S
13. Dannemann M
14. Fu Q
15. Kircher M
16. Kuhlwilm M
17. Lachmann M
18. Meyer M
19. Ongyerth M
20. Siebauer M
21. Theunert C
22. Tandon A
23. Moorjani P
24. Pickrell J
25. Mullikin JC
26. Vohr SH
27. Green RE
28. Hellmann I
29. Johnson PL
30. Blanche H
31. Cann H
32. Kitzman JO
33. Shendure J
34. Eichler EE
35. Lein ES
36. Bakken TE
37. Golovanova LV
38. Doronichev VB
39. Shunkov MV
40. Derevianko AP
41. Viola B
42. Slatkin M
43. Reich D
44. Kelso J
45. Pääbo S
(2014) The complete genome sequence of a Neanderthal from the Altai Mountains
Nature 505:43–49.

https://doi.org/10.1038/nature12886
- PubMed
- Google Scholar
1. Prüfer K
2. de Filippo C
3. Grote S
4. Mafessoni F
5. Korlević P
6. Hajdinjak M
7. Vernot B
8. Skov L
9. Hsieh P
10. Peyrégne S
11. Reher D
12. Hopfe C
13. Nagel S
14. Maricic T
15. Fu Q
16. Theunert C
17. Rogers R
18. Skoglund P
19. Chintalapati M
20. Dannemann M
21. Nelson BJ
22. Key FM
23. Rudan P
24. Kućan Ž
25. Gušić I
26. Golovanova LV
27. Doronichev VB
28. Patterson N
29. Reich D
30. Eichler EE
31. Slatkin M
32. Schierup MH
33. Andrés AM
34. Kelso J
35. Meyer M
36. Pääbo S
(2017) A high-coverage Neandertal genome from Vindija Cave in Croatia
Science 358:655–658.

https://doi.org/10.1126/science.aao1887
- PubMed
- Google Scholar
(2015) Evidence for archaic adaptive introgression in humans
Nature Reviews. Genetics 16:359.

https://doi.org/10.1038/nrg3936
- PubMed
- Google Scholar
1. Racimo F
2. Gokhman D
3. Fumagalli M
4. Ko A
5. Hansen T
6. Moltke I
7. Albrechtsen A
8. Carmel L
9. Huerta-Sánchez E
10. Nielsen R
(2017a) Archaic Adaptive Introgression in TBX15/WARS2
Molecular Biology and Evolution 34:509–524.

https://doi.org/10.1093/molbev/msw283
- PubMed
- Google Scholar
(2017b) Signatures of Archaic Adaptive Introgression in Present-Day Human Populations
Molecular Biology and Evolution 34:296–317.

https://doi.org/10.1093/molbev/msw216
- PubMed
- Google Scholar
1. Ragsdale AP
2. Gravel S
(2019) Models of archaic admixture and recent history from two-locus statistics
PLOS Genetics 15:e1008204.

https://doi.org/10.1371/journal.pgen.1008204
- PubMed
- Google Scholar
1. Reich D
2. Green RE
3. Kircher M
4. Krause J
5. Patterson N
6. Durand EY
7. Viola B
8. Briggs AW
9. Stenzel U
10. Johnson PL
11. Maricic T
12. Good JM
13. Marques-Bonet T
14. Alkan C
15. Fu Q
16. Mallick S
17. Li H
18. Meyer M
19. Eichler EE
20. Stoneking M
21. Richards M
22. Talamo S
23. Shunkov MV
24. Derevianko AP
25. Hublin JJ
26. Kelso J
27. Slatkin M
28. Pääbo S
(2010) Genetic history of an archaic hominin group from Denisova Cave in Siberia
Nature 468:1053.

https://doi.org/10.1038/nature09710
- PubMed
- Google Scholar
1. Sanchez T
2. Cury J
3. Charpiat G
4. Jay F
(2021) Deep learning for population size history inference: design, comparison and combination with approximate bayesian computation
Molecular Ecology Resources 1:13224.

https://doi.org/10.1111/1755-0998.13224
- Google Scholar
1. Sankararaman S
2. Mallick S
3. Dannemann M
4. Prüfer K
5. Kelso J
6. Pääbo S
7. Patterson N
8. Reich D
(2014) The genomic landscape of Neanderthal ancestry in present-day humans
Nature 507:354–357.

https://doi.org/10.1038/nature12961
- PubMed
- Google Scholar
(2016) The Combined Landscape of Denisovan and Neanderthal Ancestry in Present-Day Humans
Current Biology : CB 26:1241–1247.

https://doi.org/10.1016/j.cub.2016.03.037
- PubMed
- Google Scholar
1. Schrider DR
2. Kern AD
(2018) Supervised Machine Learning for Population Genetics: A New Paradigm
Trends in Genetics : TIG 34:301–312.

https://doi.org/10.1016/j.tig.2017.12.005
- PubMed
- Google Scholar
1. Scott TM
2. Guo H
3. Eichler EE
4. Rosenfeld JA
5. Pang K
6. Liu Z
7. Lalani S
8. Bi W
9. Yang Y
10. Bacino CA
11. Streff H
12. Lewis AM
13. Koenig MK
14. Thiffault I
15. Bellomo A
16. Everman DB
17. Jones JR
18. Stevenson RE
19. Bernier R
20. Gilissen C
21. Pfundt R
22. Hiatt SM
23. Cooper GM
24. Holder JL
25. Scott DA
(2020) BAZ2B haploinsufficiency as a cause of developmental delay, intellectual disability, and autism spectrum disorder
Human Mutation 41:921–925.

https://doi.org/10.1002/humu.23992
- PubMed
- Google Scholar
1. Setter D
2. Mousset S
3. Cheng X
4. Nielsen R
5. DeGiorgio M
6. Hermisson J
(2020) VolcanoFinder: Genomic scans for adaptive introgression
PLOS Genetics 16:e1008867.

https://doi.org/10.1371/journal.pgen.1008867
- PubMed
- Google Scholar
1. Sheehan S
2. Song YS
(2016) Deep Learning for Population Genetic Inference
PLOS Computational Biology 12:e1004845.

https://doi.org/10.1371/journal.pcbi.1004845
- PubMed
- Google Scholar
Preprint
(2014) Deep inside convolutional networks: visualising image classification models and saliency maps
arXiv.

https://arxiv.org/abs/1312.6034
- Google Scholar
1. Speidel L
2. Forest M
3. Shi S
4. Myers SR
(2019) A method for genome-wide genealogy estimation for thousands of samples
Nature Genetics 51:1321–1329.

https://doi.org/10.1038/s41588-019-0484-x
- PubMed
- Google Scholar
Preprint
(2015) Striving for simplicity: the all convolutional net
arXiv.

https://arxiv.org/abs/1412.6806
- Google Scholar
(2018) Model-based detection and analysis of introgressed neanderthal ancestry in modern humans
Molecular Ecology 27:3873–3888.

https://doi.org/10.1111/mec.14565
- PubMed
- Google Scholar
1. Suarez-Gonzalez A
2. Hefer CA
3. Christe C
4. Corea O
5. Lexer C
6. Cronk QC
7. Douglas CJ
(2016) Genomic and functional approaches reveal a case of adaptive introgression from Populus balsamifera (balsam poplar) in P. trichocarpa (black cottonwood)
Molecular Ecology 25:2427–2442.

https://doi.org/10.1111/mec.13539
- PubMed
- Google Scholar
(2019) Estimating the Genome-wide Mutation Rate with Three-Way Identity by Descent
American Journal of Human Genetics 105:883–893.

https://doi.org/10.1016/j.ajhg.2019.09.012
- PubMed
- Google Scholar
1. Torada L
2. Lorenzon L
3. Beddis A
4. Isildak U
5. Pattini L
6. Mathieson S
7. Fumagalli M
(2019) ImaGene: a convolutional neural network to quantify natural selection from genomic data
BMC Bioinformatics 20:337.

https://doi.org/10.1186/s12859-019-2927-x
- PubMed
- Google Scholar
Preprint
1. Turner R
2. Hung J
3. Frank E
4. Saatci Y
5. Yosinski J
(2019) Metropolis-Hastings generative adversarial networks
arXiv.

https://arxiv.org/abs/1811.11357
- Google Scholar
1. Uchiyama Y
2. Yanagisawa K
3. Kunishima S
4. Shiina M
5. Ogawa Y
6. Nakashima M
7. Hirato J
8. Imagawa E
9. Fujita A
10. Hamanaka K
11. Miyatake S
12. Mitsuhashi S
13. Takata A
14. Miyake N
15. Ogata K
16. Handa H
17. Matsumoto N
18. Mizuguchi T
(2018) A novel CYCS mutation in the α-helix of the CYCS C-terminal domain causes non-syndromic thrombocytopenia
Clinical Genetics 94:548–553.

https://doi.org/10.1111/cge.13423
- PubMed
- Google Scholar
1. Veeramah KR
2. Johnstone L
3. Karafet TM
4. Wolf D
5. Sprissler R
6. Salogiannis J
7. Barth-Maron A
8. Greenberg ME
9. Stuhlmann T
10. Weinert S
11. Jentsch TJ
12. Pazzi M
13. Restifo LL
14. Talwar D
15. Erickson RP
16. Hammer MF
(2013) Exome sequencing reveals new causal mutations in children with epileptic encephalopathies
Epilepsia 54:1270–1281.

https://doi.org/10.1111/epi.12201
- PubMed
- Google Scholar
1. Vernot B
2. Tucci S
3. Kelso J
4. Schraiber JG
5. Wolf AB
6. Gittelman RM
7. Dannemann M
8. Grote S
9. McCoy RC
10. Norton H
11. Scheinfeldt LB
12. Merriwether DA
13. Koki G
14. Friedlaender JS
15. Wakefield J
16. Pääbo S
17. Akey JM
(2016) Excavating Neandertal and Denisovan DNA from the genomes of Melanesian individuals
Science 352:235–239.

https://doi.org/10.1126/science.aad9416
- PubMed
- Google Scholar
1. Vernot B
2. Akey JM
(2014) Resurrecting surviving Neandertal lineages from modern human genomes
Science 343:1017–1021.

https://doi.org/10.1126/science.1245938
- PubMed
- Google Scholar
1. Villanea FA
2. Schraiber JG
(2019) Multiple episodes of interbreeding between Neanderthal and modern humans
Nature Ecology & Evolution 3:39.

https://doi.org/10.1038/s41559-018-0735-8
- PubMed
- Google Scholar
Preprint
1. Wang Z
2. Wang J
3. Kourakos M
4. Hoang N
5. Lee HH
6. Mathieson I
7. Mathieson S
(2020) Automatic inference of demographic parameters using generative adversarial networks
bioRxiv.

https://doi.org/10.1101/2020.08.05.237834
- Google Scholar
(2006) Adaptive introgression of herbivore resistance traits in the weedy sunflower Helianthus annuus
The American Naturalist 167:794–807.

https://doi.org/10.1086/504606
- PubMed
- Google Scholar
1. Yang XH
2. Huang S
(1999) PFM1 (PRDM4), a new member of the PR-domain family, maps to a tumor suppressor locus on human chromosome 12q23-q24.1
Genomics 61:319–325.

https://doi.org/10.1006/geno.1999.5967
- PubMed
- Google Scholar
Conference
(2019)
INVASE: instance-wise variable selection using neural networks

International Conference on Learning Representations. pp. 1–24.
- Google Scholar
1. Zammit NW
2. Siggs OM
3. Gray PE
4. Horikawa K
5. Langley DB
6. Walters SN
7. Daley SR
8. Loetsch C
9. Warren J
10. Yap JY
11. Cultrone D
12. Russell A
13. Malle EK
14. Villanueva JE
15. Cowley MJ
16. Gayevskiy V
17. Dinger ME
18. Brink R
19. Zahra D
20. Chaudhri G
21. Karupiah G
22. Whittle B
23. Roots C
24. Bertram E
25. Yamada M
26. Jeelall Y
27. Enders A
28. Clifton BE
29. Mabbitt PD
30. Jackson CJ
31. Watson SR
32. Jenne CN
33. Lanier LL
34. Wiltshire T
35. Spitzer MH
36. Nolan GP
37. Schmitz F
38. Aderem A
39. Porebski BT
40. Buckle AM
41. Abbott DW
42. Ziegler JB
43. Craig ME
44. Benitez-Aguirre P
45. Teo J
46. Tangye SG
47. King C
48. Wong M
49. Cox MP
50. Phung W
51. Tang J
52. Sandoval W
53. Wertz IE
54. Christ D
55. Goodnow CC
56. Grey ST
(2019) Denisovan, modern human and mouse TNFAIP3 alleles tune A20 phosphorylation and immunity
Nature Immunology 20:1299–1310.

https://doi.org/10.1038/s41590-019-0492-0
- PubMed
- Google Scholar
Software
1. Zarr Development Team
(2020) Zarr, version 2.4.0
Zarr.

https://zarr.readthedocs.io/en/stable/
(2020) The Impact of Recessive Deleterious Variation on Signals of Adaptive Introgression in Human Populations
Genetics 215:799–812.

https://doi.org/10.1534/genetics.120.303081
- PubMed
- Google Scholar

Article and author information

Author details

Graham Gower

Lundbeck GeoGenetics Centre, Globe Institute, Faculty of Health and Medical Sciences, University of Copenhagen, Copenhagen, Denmark

Contribution
Software, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

For correspondence
graham.gower@gmail.com

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-6197-3872
Pablo Iáñez Picazo

Lundbeck GeoGenetics Centre, Globe Institute, Faculty of Health and Medical Sciences, University of Copenhagen, Copenhagen, Denmark

Contribution
Software, Formal analysis, Visualization, Methodology, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-7174-3264
Matteo Fumagalli

Department of Life Sciences, Silwood Park Campus, Imperial College London, London, United Kingdom

Contribution
Conceptualization, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4084-2953
Fernando Racimo

Lundbeck GeoGenetics Centre, Globe Institute, Faculty of Health and Medical Sciences, University of Copenhagen, Copenhagen, Denmark

Contribution
Conceptualization, Supervision, Funding acquisition, Methodology, Writing - original draft, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-5025-2607

Funding

Villum Fonden (00025300)

Fernando Racimo

Leverhulme Trust (RPG-2018-208)

Matteo Fumagalli

Lundbeckfonden (R302-2018-21555)

Fernando Racimo

Novo Nordisk Fonden (NNF18SA0035006)

Fernando Racimo

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

Acknowledgements

We thank Andrew Kern, Martin Sikora, Flora Jay and Anders Albrechtsen, as well as three anonymous reviewers, and the members of the Racimo group and the PopSim consortium, for helpful advice and discussions. We also thank Murray Cox and Georgi Hudjashov for facilitating access to the IGDP data. Funding FR and GG were supported by a Villum Fonden Young Investigator award to FR (project no. 00025300). MF was supported by a Leverhulme Research Project grant (RPG-2018–208) and the Imperial College European Partners Fund. FR was also supported by a Lundbeckfonden grant (R302-2018-2155) and a Novo Nordisk Fonden grant (NNF18SA0035006) to the GeoGenetics Centre. The funding sources were not involved in study design, data collection and interpretation, or the decision to submit the work for publication.

Copyright

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

Metrics

5,111

views
546

downloads
68

citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

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

Downloads (link to download the article as PDF)

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

Mendeley

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

Graham Gower
Pablo Iáñez Picazo
Matteo Fumagalli
Fernando Racimo

(2021)

Detecting adaptive introgression in human evolution using convolutional neural networks

eLife 10:e64669.

https://doi.org/10.7554/eLife.64669

Categories and tags

Research organism

Human

Share this article

Cite this article

A schematic overview of how genomatnn detects adaptive introgression.

CNN performance on validation simulations for Demographic Model A.

Saliency maps, showing the CNN’s attention across the input matrices for each simulated scenario, calculated for the CNN trained on Demographic Model A, filtered for beneficial allele frequency >0.25.

Comparison of Manhattan plots using beta-calibrated output probabilities for different class ratios.

Top ranking gene candidates corresponding to Neanderthal AI in Europeans.

Application of the trained CNN to the Vindija and Altai Neanderthals, and 1000 genomes populations YRI and CEU.

Top ranking gene candidates corresponding to Denisovan AI in Melanesians.

Application of the trained CNN to the Altai Denisovan and Altai Neanderthal, 1000 genomes YRI populations, and IGDP Melanesians.

Top ranking gene candidates corresponding to Neanderthal AI in Europeans.

Top ranking gene candidates corresponding to Denisovan AI in Melanesians.

Loss and accuracy for CNNs after training for three epochs, as reported by Keras/Tensorflow, for the training and validation datasets.

Parameter values used for simulating Demographic Model A1.

Haplotype plot for the candidate region chr1:104500001–104600000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr2:109360001–109460000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr2:160160001–160280000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr3:114480001–114620000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr4:54240001–54340000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr5:39220001–39320000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr6:28180001–28320000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr8:143440001–143560000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr9:16700001–16820000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr12:85780001–85880000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr19:20220001–20380000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr19:33580001–33740000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr20:62100001–62280000 in the Neanderthal-into-European AI scan.

Haplotype plot for the candidate region chr21:25840001–25940000 in the Neanderthal-into-European AI scan.

Genotype plot for the candidate region chr2:129960001–130060000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr3:3740001–3840000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr4:41980001–42080000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr5:420001–520000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr6:74640001–74740000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr6:81960001–82060000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr6:137920001–138120000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr7:25100001–25200000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr7:38020001–38120000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr7:121160001–121260000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr8:3040001–3140000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr12:84640001–84740000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr12:108240001–108340000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr12:114020001–114280000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr14:61860001–61960000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr14:63120001–63220000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr14:96700001–96820000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr15:55260001–55400000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr16:62600001–62700000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr16:78360001–78460000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr18:22060001–22160000 in the Denisovan-into-Melanesian AI scan.

Genotype plot for the candidate region chr22:19040001–19140000 in the Denisovan-into-Melanesian AI scan.

Author details

Graham Gower

Contribution

For correspondence

Competing interests

Pablo Iáñez Picazo

Contribution

Competing interests

Matteo Fumagalli

Contribution

Competing interests

Fernando Racimo

Contribution

Competing interests

Downloads (link to download the article as PDF)

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

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

Categories and tags

Research organism

Further reading