1. Introduction

Sequencing of the Neanderthal genomes [1] revealed that Neanderthals contributed to around 2% of the genomes of present-day humans with non-African ancestry [2,3]. The discovery and subsequent sequencing of another archaic human group [4,5], called the Denisovans, revealed that they too contributed to present-day human genome pool, with around 0.1% to Eurasian populations and 2-3% to people from Papua New Guinea [6]. Both Neanderthal and Denisovan DNA have impacted many human traits, including health, particularly the immune system [7-9].

Since most tools to detect archaic introgression are based on finding homologous sequences between sequenced archaic genomes and modern human samples, every time a new archaic genome is sequenced, introgression maps need to be updated [10]. Moreover, since the limited number of sequenced archaic human genomes might not represent all genetic variation from the introgressed population, the results from those methods are biased towards the sequenced archaic genomes to date.

To overcome this limitation, hmmix—pronounced h-m-m-mix—determines the archaic contribution to present-day human genomes without the need for archaic human reference genomes [11].

2. How the Model Works

The method was designed for finding introgressed segments from a source population (referred to as Archaic, Figure 1) in a test/target genome (referred to as Ingroup, Figure 1). This is particularly useful when the source population is not sampled either because it is extinct, when DNA has yet to be recovered from it or the samples available are substantially diverged from the introgressed population, as is the case for Denisovan introgression. In this particular case, the sampled Denisovan genome diverged hundreds of thousands of years from the introgressed populations [12]. The method’s trade-off is to assume that we have access to a set of individuals closely related to the ingroup without introgression from the source population (referred to as Outgroup, Figure 1), that share a great proportion of variation with the ingroup. Thus, hmmix requires only genetic data from the ingroup and the outgroup populations, but no archaic reference sequence is needed (Figure 1). Even though hmmix terminology - archaic, ingroup and outgroup - refers to the archaic introgression scenario, hmmix can be used more generally to detect introgression from any donor/source population into a test/target population only if there exists an outgroup/unadmixed population.

Figure 1 Overview of the method. Diagram of two gene trees for an introgressed segment (blue) and a non-introgressed segment (orange) for an ingroup individual and five outgroup individuals. When removing variants in the ingroup genome (orange and blue rectangle below) found also in the outgroup (grey) only variants private to the ingroup persist (dark orange). The difference in private variant density is how the model distinguishes ingroup segments from archaic introgressed segments.

Another advantage of the method, compared to other reference-free approaches, is that it is run independently for each ingroup genome, setting it apart from LD-based methods such as S* [13] and S’ [14]—Although we note that later versions of S* can be run on single genomes [15,16]. This allows the method to discover archaic segments independently of their frequency in the ingroup population.

To identify introgressed segments with hmmix, the method first removes the variation shared between the ingroup and the outgroup (Figure 1). Depending on the demographic history and sample size of the outgroup, this step can eliminate most of the variants that originated in their common ancestor. For instance, if the outgroup population is structured, the likelihood of identifying shared variants increases, especially when a great number of outgroup genomes are sampled. We have created a github repository (https://github.com/LauritsSkov/Determine_outgroup_size/) for the purpose of determining how much shared variation is removed depending on demography and sample size based on simulations.

The remaining SNPs in the ingroup individual are considered private SNPs. Their density, measured in genomic windows, is proportional to the minimum coalescence time with the outgroup: archaic segments show a higher SNP density than non-archaic segments because of their deeper coalescent times (T_archaic > T_ingroup, Figure 1). The power and accuracy of hmmix is high if the coalescence time between introgressed archaic segments and non-introgressed segments are very different (T_archaic >> T_ingroup).

In the context of Neanderthal introgression, the split time between present-day Sub-saharan Africans and non-Africans is estimated at around 50,000–70,000 years ago [6,17-19] while the split between modern humans and Neanderthals/Denisovans is thought to be 580,000 years ago [10]. Thus, the SNP density in the archaic state is roughly twelve times higher than the human ingroup state, which makes hmmix a good tool to infer archaic introgression in this scenario [11].

While the outgroup in hmmix is typically defined as a population free of archaic introgression—often Sub-Saharan Africans—the method is flexible. By adjusting the choice of outgroup, users can target population-specific admixture events. For instance, Denisovan introgression in Papuans can be investigated by adding any population unaffected by that specific Denisovan gene flow to the outgroup, such as European, or Asian genomes [11]. In this setup, hmmix highlights archaic components private to Papuans relative to the chosen outgroup. However, not all archaic content found with this method corresponds to the Papuan-specific Denisovan admixture. Any Neanderthal ancestry private to Papuans, when absent from the outgroup either due to genetic drift or other factors, would also be included in the inferred archaic component found with this method.

Initially, the emission parameters and transition parameters between states are not known and have to be learned from the data using hmmix’s EM optimization (Baum–Welch training). Once the parameters are optimized, the user can identify the location of archaic segments with the decoding methods provided by hmmix. The model also takes local mutation rate and genomic callability into account in its parametrization. Finally, the classification of archaic segments can be done post-processing by comparing the segments to sequenced archaic genomes.

Another feature of the hmmix framework is that the trained HMM parameters have biological interpretations: The HMM emission parameters are proportional to the coalescence times with the outgroup (T_archaic and T_ingroup) and the HMM transition parameters relate to the admixture proportion and admixture time [11].

3. How to Run the Method and Online Resources

An in-depth tutorial on how to use hmmix is hosted on github: https://github.com/LauritsSkov/Introgression-detection/.

This tutorial covers in depth the following points and it uses an individual from the publicly available 1000 genomes project as an example. The output at each step is shown and available so the user can verify that hmmix is working correctly.

4. Limitations and Recommendations

Here, we discuss issues arising from the underlying demography of the target populations. We assume that variants can be whole genome genotyped accurately. Any technical issue in this regard, would limit the performance of the method. In general, sequencing errors would decrease the power of the tool, making SNP density more similar between archaic and non-archaic states. Additionally, sparse data such as low coverage data (e.g, aDNA) or array genotyping treated without proper care leads to inaccurate segment calls with fragmented segments and increased false negative rate. For the aDNA case, we recommend using tools specifically conceived for this purpose [30].

hmmix makes multiple assumptions about the demography of the analyzed groups (Figure 2a). If these assumptions are not met, hmmix is not advised to be used. First, hmmix is not able to detect archaic introgression present in both the ingroup and outgroup population (Figure 2b). This is because these segments are removed when shared variation with the outgroup is filtered out from the ingroup genome. If you suspect that kind of introgression in your system, it is advised to use an approach that does not require genomes from the unadmixed outgroup such as ARGweaver-D [31] or ArchIE [32], both requiring instead prior knowledge about the underlying demographic model. Alternatively, you can use other models that do not rely on an unadmixed outgroup but do use archaic reference genomes instead, such as IBDmix [33].

Figure 2 How demography and outgroup size changes coalescence times for archaic and ingroup states. In all panels, the msprime [34] simulated demography and the resulting distribution of true coalescence times for the ingroup state (orange) and archaic state (blue) are shown. Note this demography does NOT correspond to neanderthal/human introgression but is chosen for illustrative purposes. The y-axis is shared for the demography illustration and the coalescence time distribution plot. All panels are truncated at 600 thousand years ago. If the coalescence time distributions are very distinct, hmmix classifies introgressed segments more accurately. Each panel shows 10,000 simulated segments (50kb long) without recombination using msprime. The split time between ingroup and outgroup is 50,000 years ago and the split between archaic and outgroup is 150,000 years ago. The effective population size of all population size is shown for each population at the bottom. The dotted horizontal lines indicate the estimated emission parameter for the ingroup (orange) and archaic (blue) state by hmmix after Baum-Welsch training. The solid lines represent the mean true simulated coalescence time for each state. (a) Representation of the base-line scenario used for comparison with other scenarios. (b) introgression into the common ancestor of the ingroup and outgroup is added. (c) Geneflow between the ingroup and outgroup is added in the model, leading to many introgressed archaic segments also being present in the outgroup. (d) The effective population size of ingroup and outgroup are increased to 20k while the archaic population is kept at 2k. This increases the coalescence times especially for the archaic/outgroup state which is 856 thousand years ago. (e) The number of individuals available from the outgroup is 1, which leads to a greater overlap of coalescence times for archaic and non-archaic segments. (f) No gene flow from archaic to ingroup is simulated here. hmmix still attempts to fit two states. Note there are two inferred states (two dotted horizontal lines) even though there is only one true state.

Second, it assumes that there is none or limited geneflow in the direction ingroup > outgroup. This kind of introgression transfers archaic variation into the outgroup, which leads to a similar scenario as the one above (Figure 2c). There is evidence that back migrations from non-African populations to Africa likely masked introgressed Neanderthal segments in non-Africans [33], but not to an extent that archaic introgressed segments are completely removed. This effect can be reduced by selecting outgroup individuals with limited evidence of admixture with the ingroup population. Geneflow from outgroup > ingroup is not a concern as this only lowers T_ingroup while keeping the same T_archaic, thus not affecting the inference of archaic segments in the ingroup.

Third, the model increases its power when the distributions of coalescence times with the outgroup are sufficiently different between introgressed and non-introgressed states. In scenarios where the archaic population splits off soon before the split of the ingroup/outgroup population (T_archaic ~ T_ingroup), hmmix has trouble finding the correct archaic segments. It is also not possible if we are unable to remove enough variation that happened in the common ancestral population between ingroup and outgroup. This could be the case when there is a very large ingroup and outgroup ancestral population (Figure 2d) or a low number of samples from the outgroup are available (Figure 2e).

The demographic history of the outgroup also affects how much common variation can be removed from the ingroup. Scenarios in which outgroup lineages are prevented from coalescing within the outgroup, increase the change of coalescence with the ingroup in the common ancestor of ingroup/outgroup, which translates to more variants shared between ingroup and outgroup. Rapid coalescence in the common ancestor of ingroup/outgroup increases the contrast between T_archaic and T_ingroup. Overall, it is ideal if the outgroup population is large and structured and the common ancestor of ingroup/outgroup has a small population size.

Lastly, the coalescence times between ingroup and outgroup do not happen at the same time point, as hmmix assumes, but instead are distributed across time—exponentially distributed in a neutral scenario -, making some coalescences much older than the mean (Figure 2f). In general, this is not a problem for hmmix and can only create a few false positives. However, in the absence of introgression, hmmix always fits two emission parameters, or as many as the user defines, corresponding to each hidden state (archaic human and modern human) from the continuous distribution of coalescences, i.e. overfitting the data (Figure 2f). Currently, there is no guideline to check if the ingroup has introgression or not as hmmix is not a formal test of admixture, rather it splits the genome into how many states the user defines and infers their emission parameters.

One way to test if these states are biologically meaningful is to perform computer simulations under multiple demographies and compare the emission and transition parameters obtained after the training stage. One should also investigate the length distribution of decoded segments - false positive segments which are due to deep coalescence within the common ancestor of the ingroup and outgroup tend to be short due to recombination.

Alternatively, other genome-wide tests can be performed in advance such as f3 or f4 statistics. We stress the importance of this step as one might otherwise mistakenly infer that there is introgression from an unsampled population. We refer to this review on methods for inferring introgression [35].

To allow users of hmmix to easily simulate data from a chosen demography we make a github repository available with the code to run the above examples and encourage users to play with the different parameters to test how the coalescence time distributions differ for a specific demography that is relevant for the tested population (https://github.com/LauritsSkov/SimulateDemography/).

5. Closing Remarks

hmmix is a powerful tool to detect introgression without a reference genome of the introgressing population on a demography similar to the modern and archaic human admixture scenario. hmmix is being maintained and new features are still being added. We recently updated the tool with two new modules, enabling it to perform hybrid decoding and sampling hidden state sequences from the posterior probabilities [20,21]. In conclusion, hmmix has been successfully used for 8 years and as more genomes from different species become available we anticipate that it will continue to be of use to find introgression in the future.

Declarations

References