Y-mer: A k-mer based method for determining human Y chromosome haplogroups from ultra-low sequencing depth data

doi:10.21203/rs.3.rs-5042960/v1

Download PDF

Method Article

Y-mer: A k-mer based method for determining human Y chromosome haplogroups from ultra-low sequencing depth data

https://doi.org/10.21203/rs.3.rs-5042960/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Determining genetic ancestry of an individual is challenging from poorly preserved or mixed samples that permit only ultra-low sequence (ulcWGS) depth < 0.1x to be gained at target loci. Leveraging the recent advances in telomere-to-telomere sequencing of the whole genomes with long reads we show first in a simplified example how short DNA string (k-mer) copy numbers at two different types of repeat arrays correlate with basal chromosome Y (chrY) haplogroups (HG-s). We develop a new k-mer based method Y- mer and show how information from hundreds of thousands of k-mers in distance-based models enables accurate inference of chrY haplogroup from WGS sequence at depth less than 0.01x without additional PCR or capture. We test the performance of Y-mer on ancient DNA and prenatal screening data showing its potential for genetic ancestry inference for cell free, forensic and ancient DNA research from short read WGS data.

aDNA

ulcWGS

NIPS

NIPT

k-mer

human

chrY

haplogroups

prediction

Human Y chromosome was due to its high complexity the last chromosome of the human genome to be completely sequenced (Rhie et al., 2023). Approximately half of its length is composed of euchromatin, regions recombining (PAR1 and PAR2), and regions not recombining with human X chromosome, while the other half, including its centromere and the largest heterochromatin region (Yq12) in the entire human genome, remained uncharted in the human reference genome until the recent completion of the Y chromosome assembly (Rhie et al., 2023). Skov et al. (Skov, The Danish Pan Genome Consortium and Schierup, 2017) study of 62 Danish Y chromosomes sequenced to high coverage with a wide range of library insert sizes and a study of Y chromosome assemblies of 7 major Y chromosome haplogroups by Esteller-Cucala (Esteller-Cucala et al., 2023) highlighted major gaps in the reference (GRCh38) sequence and uncovered high degree of variation in the centromeric region and structural variants restricted to certain haplogroups. Furthermore, two recent studies that used high coverage long read sequence data to determine the sequence of Y chromosome from telomere to telomere, revealed highly variable numbers of amplicon genes and up to two-fold difference in the total length of the Y chromosome (Hallast et al., 2023; Rhie et al., 2023). The main source of Y chromosome length variation is the Yq12 region, which showed a high degree of length variation not only across but also within the two most represented haplogroups, E1 and O2 in a sample of 43 Y chromosomes that have been completely sequenced (Hallast et al., 2023). Certain repeat regions, e.g. the TSPY array, showed, however, length variation that was constrained by phylogeny. To what extent differences in the size of the individual constituent blocks of Y chromosome sequence are phylogenetically informative remains to be systematically studied with larger sample sizes.

Y chromosome haplogroups are routinely determined in ancient DNA studies as they are indicative of male-specific gene flow and important for clarifying specific relationships among individuals whose relatedness has been identified from autosomal data. They are also relevant for individual identification in forensic practice. However, Y chromosome haplogroup determination from shotgun sequence data is challenging below 0.1× and 0.01× coverage because most branches in the Y phylogeny are defined by less than 10 and only a few with > 100 SNVs, respectively. The increase of sequence coverage via capture or additional shotgun sequencing may be impractical in cases of poorly preserved samples due to high fragmentation, damage or low proportion of endogenous DNA. Here, in this study, we explore the potential of predicting Y chromosome haplogroups from ultra-low coverage sequence data with a new k-mer based method Y-mer, which uses the counts of k-mers unique to the complete sequence of human Y chromosome. We will first examine the correlations between hundreds of thousands of k-mers and Y chromosome haplogroups determined from high quality SNV data in the 1000 Genomes Project males (1000 Genomes Project paper 2015) and Estonian Biobank high coverage sequences (Mitt et al., 2017). We develop lower and higher resolution models of haplogroup prediction with Y-mer in training sets of regional and global samples and determine the prediction accuracy in relevant validation sets. We determine the lower boundary of sequence coverage for the approach by down sampling the coverage of Y chromosomes of known haplogroups and, finally, apply the approach on ultra-low coverage ancient DNA and non-invasive prenatal screening (NIPS) data where low coverage Y chromosome reads of the male fetus are presented together with majority of reads derived from the mother.

We test whether and how accurately human Y chromosome haplogroups can be predicted with Y-mer from ultra-low coverage sequence data using counts of k-mers unique to human Y chromosome. As nearly one half of the Y chromosome sequence of the HG002 assembly (Rhie et al., 2023) is uniquely mappable with short male specific k-mers, corresponding to ~ 30Mbp, which is roughly 1800 times of the length of human mitochondrial genome, it is expected that such approach may have a major potential for the study of genetic ancestry of poorly preserved or mixed samples.

Simplified model of Y chromosome haplogroup clustering from three repeat arrays in chrY data from a single population

Firstly, as a proof of principle, we tested the feasibility of the prediction of Y chromosome haplogroups from the counts of k-mers using high coverage Y chromosome sequence data available for 198 Estonian male individuals in which Y chromosome haplogroups had been determined confidently and at high resolution from SNVs (Table S1). We estimated the counts of three k-mers from chrY repeat arrays in each individual, normalized by the counts of a k-mer selected from the LINE element (Table S2). Figure 1 presents the distribution of Y chromosome haplogroups in the 198 individuals by the normalized counts of two k-mers (plots based on three k-mers shown in Figure S1 and 3-dimensional plot in https://bioinfo.ut.ee/randomtandem/EGV.html). We find that the clustering of the male samples by the k-mer ratios clearly separates the basal haplogroups N3a, I, R1a, R1b, E in the Estonian sample. The clustering by haplogroups can be observed in each pairwise plot, with the usage of each additional repeat adding more resolution. For example, in the HET-CEN plot (Figure S1A), haplogroups N and R can be clearly separated from each other while it is impossible to differentiate I from N3a. The latter two haplogroups can be distinguished by the counts of the VNTR k-mer (Figure S1B,C). High copy number of the chosen CEN k-mer further separates R1a and R1b lineages from other haplogroups (Figure S1A,C) suggestive of centromere length shortening in the R1 lineage between 20–50 Kya after its formation (Karmin et al., 2015). However, analysis of the k-mer counts of these three repeat arrays in 1000G individuals with broader range of Y chromosome haplogroups and their sub-clades sampled in a geographically wide range of populations across the world revealed more complex patterns of repeat counts by haplogroups with lesser resolution (data not shown).

Complex distance-based models for predicting Y chromosome haplogroups in multi-population data

To explore whether with larger numbers of chrY k-mers it is possible to robustly predict chrY haplogroups from multiple populations across the globe we generated and tested models based on hundreds of thousands to millions of k-mers. We created models using five different sets of high coverage genomes as sources of chrY data (Fig. 2), the first with only a single chrY (HG002, Table S3) T2T long-read assembly, the second with 21 T2T Y chromosome assemblies from different haplogroups (Hallast et al. 2023, Table S4), and, finally, three sets of short-read data for 110–222 Y chromosomes from 1000 Genomes Project (The 1000 Genomes Project Consortium et al., 2015) and EGC data (Mitt et al., 2017). Using glistmaker of GenomeTester4 (Kaplinski, Lepamets and Remm, 2015) we extracted lists of k-mers specific to chrY, including up to 14 million 25-mers (see Material and Methods for details), from which sets of 10,000 to 100,000 k-mers were chosen for downstream analyses per each haplogroup used in the given model (Fig. 1). We trained haplogroup predictions on three different geographic ranges, 110W set representing 11 worldwide haplogroups, 213E representing 22 haplogroups common across Europe and 222NE set representative of haplogroups common specifically in Northeast Europe (Table S6,S7). In each model the haplogroup prediction was made with distance-based clustering by assigning individual to the haplogroup with the smallest distance considering all k-mers in the model (Appendix Model).

Model validation

We tested the performance of haplogroup prediction models in a validation set of 110 individuals (V110) selected from the 1000 Genome Project, HGDP and EGV data, including only individuals that had not been used in model training (Table S9). The haplogroup of each individual of the validation set had been determined from SNV data. The validation set included 10 individuals from each of the 11 basal haplogroups (AB, C, E, G, H, IJ, LT, N, O, Q, and R) chosen to represent human Y chromosome diversity at the global scale. Each individual in the validation set was downsampled to 10 different coverage values, within the range of 1x to 0.00001×, to test the limits of the approach.

Firstly, to clarify how many different Y chromosome assemblies are required as sources of k-mer lists to accurately determine Y chromosome haplogroups from ultra-low coverage data we compared the performance of M1W, M21W, and M110W models (Fig. 3, Table S10) on the validation set. We found that M21W and M110W models both offered high accuracy (> 0.95) from 1x down to 0.001x coverage while using only a single Y chromosome assembly as a source performed poorly at all coverages, with a 0.65 accuracy at maximum for 1× coverage. The M110W model showed higher accuracy than M21W between 0.0005× and 0.005× coverage while at higher coverage ranges the difference was not notable.

Next, to determine how many k-mers are required to accurately predict chrY haplogroups, we compared the performance of models based on 10K to 100K k-mers selected from each haplogroup. In each tested model we used the 21Y, the long-read assembly of 21 Y chromosomes (Hallast et al. 2023), as the source for k-mer lists. We found that all models performed highly accurately (> 0.95) at coverage equal or higher than 0.0005x (Fig. 4). Models based on 10K k-mers were significantly less accurate at coverage range below 0.0005x. While at the lowest coverage range, 0.00001x, the 20K model also performed poorly, we observed no major differences in the performance of models based on 30K or more k-mers on the coverage range tested.

We tested the effect of contamination on the accuracy of haplogroup prediction with Y-mer on downsampled replicas of a Finnish individual from haplogroup N3. To this test sample, we added reads from another individual, from haplogroup R1a, while keeping the total number of reads constant at coverage 0.01x. We determined with Y-mer the haplogroups of such pools of individual reads, showing that Y-mer haplogroup prediction remains highly robust with contamination rates up to 30% (Fig. 5).

Model testing on ancient DNA data

Our tests of haplogroup prediction accuracy on validation set V110 included individuals from the same modern populations from which the training sets were formed. To explore the behavior of our models on datasets with haplogroup composition different from the modern training sets we turned to aDNA data. We first tested the haplogroup prediction models on ancient genomic data of 91 male individuals from Eurasian Steppe Belt, dated to 1500–4500 years BP and whole genome sequenced to depth range of 0.029-8.7x (Damgaard et al., 2018). Beyond the temporal difference between our modern references and the Damgaard data, none of the modern Steppe Belt or Central Asian populations were included in the data that we used for model training. For 47 individuals basal Y chromosome haplogroups were reported by Damgaard et al. 2018 while a sizeable proportion of 44 male samples had no haplogroup assignment, likely due to their low coverage. To generate a list of expected haplogroups for model testing we called 256,463 binary haplogroup informative SNVs and determined the likeliest chrY haplogroups of all the 91 Steppe Belt individuals in Damgaard et al. data (Table S11) following the approach described in Hui et al. 2024 (Hui et al., 2024). All haplogroup assignments we made for the 44 individuals with previous assignments matched at the base haplogroup level those made by Damgaard et al. 2018 (Damgaard et al., 2018). A number of haplogroups detected in this ancient data set, such as C3, I3, N5, O6a, Q1c, Q1g, R1b13, and R1b16, are either extremely rare in modern data sets and/or not included in our training sets.

To test the performance of Y-mer on the phylogenetically diverse ancient Steppe Belt data, we predicted the chrY haplogroups for all the 91 individuals using 6 different k-mer based models (Table S12). At the most basal level of 11 haplogroups predicted from k-mer lists derived from 21 Y chromosome T2T assemblies (model M21W), we observed high rate (94%) of matches between the haplogroup calls made from SNV and k-mer data (Table 1). In models trained to predict haplogroups at higher resolution (M21E, M222E, M21NE, M222NE) we observed accuracy range 72–86% with more mismatches, in particular in haplogroups not used in model training, such as C and O in the M21NE and M222NE models. In the case of model M222NE, all 78 individuals whose SNV-determined haplogroup was included in the model were predicted correctly. In the 13 cases where the haplogroup expected from SNV data was not used in our model training, the model prediction was in 6 cases a phylogenetically closely related lineage (such as I2 instead of I3 and R1b3 instead of R1b16) and in 7 cases haplogroup LT, which was phylogenetically more distant than the closest haplogroup in the training set. Instances where haplogroup predictions were incorrect included, notably, haplogroup LT assignments made with high confidence (p < 0.05), whereas only a small subset of wrong predictions were due to low confidence calls. As filtering for the p-values did not appear to increase the accuracy of haplogroup predictions we considered the consistency of haplogroup predictions made by different models. We found that for 45 individuals, the haplogroup assignments were consistent by all 6 models we examined and these predictions were always correct. Notably, the correct predictions include the individual with the lowest coverage. Using p < 0.05 threshold for haplogroup predictions improved but only modestly the accuracy (Table 1).

Table 1

Accuracy of haplogroup assignment models on three ancient DNA data sets.
		Gretzinger		Damgaard		Saag
	p-value	no thr	< 0.05	no thr	< 0.05	no thr	< 0.05
Haplogroup assignment models
M21W	N	252	180	91	81	23	17
	Accuracy	0.571	0.700	0.945	0.963	0.826	0.941
M110W	N	252	225	91	82	23	21
	Accuracy	0.900	0.970	0.835	0.866	0.913	0.952
M21E	N	250	225	70	57	23	17
	Accuracy	0.784	0.822	0.771	0.860	0.826	0.941
M213E	N	250	248	70	68	23	23
	Accuracy	0.992	0.992	0.929	0.941	1	1
M21NE	N	123	78	26	24	15	14
	Accuracy	0.537	0.692	1	1	0.333	0.357
M222NE	N	123	70	26	26	15	8
	Accuracy	0.423	0.529	1	1	0.600	0.625
Sub-clade assignment models
M80R1	N	75	48			13	11
	Accuracy	0.853	0.979			0.846	0.909
M43I1	N	15	10
	Accuracy	0.800	1

Note, N – number of individuals considered in tests, ‘no thr’ – no p-value threshold used; ‘<0.05’ - keeping only haplogroup assignments supported by a p-value < 0.05. Gretzinger data (Gretzinger et al., 2022); Damgaard data (Damgaard et al., 2018); Saag data (Saag et al., 2017, 2019).

To test the accuracy of the performance of higher resolution models M21E, M213E, M21NE and M222NE on data with haplogroup composition more similar to our training set we turned to two additional ancient DNA data sets, one from Early Medieval Northwest Europe (Gretzinger et al., 2022) and one from Bronze Age and Iron Age Estonia (Saag et al., 2017, 2019). These data sets include individuals separated from our training sets by time while still having similar chrY haplogroup composition to them at subclade level (Table S13-14). We find limited or no improvement compared to the Damgaard et al. data in the performance of M21E, M21NE, and M222NE models (Table 1). However, the M213E model yielded high (> 0.98) accuracy of basal haplogroup assignments in both Gretzinger et al. and Saag et al. data sets. As the M213E model differs from M222NE by its lack of resolution between recently diverged R1b and I1 sub-clades, it is notable that the accuracy of haplogroup predictions made after applying M80R1 and M43I1 models separately on individuals assigned to R1 and I1 clades, respectively, was considerably higher (> 0.8) than in case of the M222NE model (< 0.53), in which these sub-clade assignments were made in context of higher level haplogroup assignments of the global chrY phylogeny.

Model testing on non-invasive prenatal screening (NIPS) data

Another source of low coverage data from which Y chromosome prediction with conventional SNV-based methods would be extremely challenging are non-invasive prenatal screening (NIPS) data generated from circulating free DNA (cfDNA) samples of maternal blood widely used to identify chromosomal structural abnormalities at various stages of fetal development. Using NIPS data of pregnant individuals from China (Xu et al., 2018) and Estonia (Žilina et al., 2019) and focusing on male fetus cases, where the fetal sex had been independently determined, we show that Y-mer is able to predict with high confidence fetal haplogroups from sequencing depths 0.006-0.12x in Estonian samples (Table S15) and 0.0009-0.009x in Chinese samples (Table S16). While we do not have reference data for individual NIPS samples to estimate the accuracy of Y-mer, the haplogroup predictions in both Chinese and Estonian cohorts appear to match for most haplogroups their frequency expectations for China and Estonia made from independently sampled cohorts (Table S17). In case of Chinese NIPS haplogroup predictions we observe overall good match between predicted and observed haplogroup frequencies with somewhat higher than expected frequency of C and lower than expected O3 (O2a2 by ISOGG nomenclature) frequencies, which could reflect differences in recruitment in the NIPS and 1000 Genomes Project data of individuals from Southern and Northern China where these haplogroups show opposing frequency gradients (Li et al., 2023).

Determining individual ancestry from traces of DNA derived from forensic context or from poorly preserved ancient human remains from archaeological context is typically challenging due to low number of host DNA molecules in the sample, their fragmentation and post-mortem damage. Because of the lack of recombination in Y chromosome, its variation accumulates over time along a simple phylogenetic tree. The main branches of this tree can be robustly inferred even from low coverage data because of the hierarchic redundance of the accumulated variation. Furthermore, higher Y chromosome regional differentiation observed among present-day populations, compared to mitochondrial and autosomal diversity (Karmin et al., 2015), makes Y chromosome haplogroup prediction attractive for genetic ancestry inference. In this study we have shown that human Y chromosome haplogroups can be predicted accurately and efficiently from ultralow coverage sequence data using methods that determine the relative abundance of male-specific k-mers. With a simplified example, illustrated on present-day Estonian Biobank chrY data generated at high sequence depth and quality (Mitt et al., 2017), we showed that the three most common basal haplogroups in Estonia, N3, R1, and I1, can be differentiated from each other on the basis of the k-mer abundance of just two chrY k-mers from repetitive regions DYZ3 and DYZ19. The observed haplogroup-specific differences in k-mer abundance could reflect structural changes in the phylogenetic ancestry of the given haplogroups, e.g. shortening of the centromeric DYZ3 region in the ancestry of haplogroup R1 and shortening of the DYZ19 region in the I1 ancestry. Differentiation at sub-clade level and separation of minor haplogroups J2, I2 and E2 was, however, not possible with these two k-mers only. Furthermore, haplogroup prediction models based only on a small number of repeats are likely to be vulnerable to parameters such as sequencing depth, which motivated us to explore models based on a larger number of k-mers.

The accuracy, resolution, and robustness of haplogroup prediction, whether using the SNVs or the Y-mer approach described here, would depend, besides the quality and sequence coverage of the sample, on the size and diversity of the reference panel as well as the number of informative variants being used. Our tests with global and European training and validation sets showed that models using just a single T2T Y chromosome reference genome as the source of k-mer extraction performed poorly, with lower than 80% accuracy across the tested coverage range of down-sampled validation data. Models using k-mers extracted from 21, 110, and 213 different Y chromosome sources performed better, with accuracy higher than 95%, in coverage range > 0.001x (Fig. 3). As we did not observe major differences in the performance of M21W versus M213W and M21E versus M213E models, we can conclude that a phylogenetically diverse panel of more than 20 Y chromosomes suffices as a k-mer source for basal haplogroup determination.

We showed that haplogroup predictions with Y-mer are robust to contamination in a model (M21E) trained with individuals from phylogenetically distinct haplogroups (Fig. 5). However, models trained with regionally specific sets of haplogroups (M21NE), including sub-clades of I1 and R1b that have separated only within the last 5,000 years (Karmin et al. 2015), performed less well even without the presence of contamination (Table 1). This drop of performance is likely caused by increasingly higher proportions of k-mers that overlap the k-mer lists extracted for the given sub-clades. This result highlights the need to use additional filters that remove k-mer overlaps in future developments of detailed sub-clade prediction models, or, with our current approach, the need to use both diverse and balanced training sets. Even though we did not see improvement between M21W and M110W models (Fig. 3) it is possible that with models that require higher haplogroup resolution the number of Y chromosome sources from which k-mers are retrieved will also need to be adapted.

Our analyses revealed that for robust haplogroup prediction the number of k-mers selected to distinguish each individual haplogroup included in the model has to be sufficiently high (> 10,000). Our comparisons of models with increasingly higher numbers of k-mers, however, showed no detectable increase in accuracy (Fig. 4) in models with more than 20,000 haplogroup-specific k-mers, suggesting that for computation efficiency models with 20,000–50,000 k-mers are likely to represent the most optimal solution.

When applying Y-mer on validation sets different in their haplogroup composition from the training set we observed higher rates of mismatches, particularly in association with rare sub-clades. Predictions supported by models trained at different haplogroup levels appeared to be mostly correct suggesting that applying multiple models on the same data can be helpful in distinguishing predictions that are robustly supported by multiple models from those that have lower confidence and are supported only by individual predictions. The necessity to apply multiple models on the same data was further illustrated in case of the ancient DNA data from the Steppe Belt, where our basal haplogroup prediction model M21W made haplogroup calls at high accuracy (~ 95%), while predictions with region-specific haplogroup compositions, adapted on Europe or more specifically in case of some models on Northeast European data, showed higher number of mismatches, particularly in case of haplogroups that were not included in the model. While we could see that haplogroup predictions that were supported by multiple models appeared mostly to be correct, this case study highlighted the need also for caution in choosing the models with appropriate haplogroup composition and the training sets for the future use of Y-mer tool in ancient DNA studies. Preliminary insights of haplogroup composition of an ancient cohort through SNV analyses of higher coverage individual samples will be advisable to inform such models. Where such high coverage data is obtainable it can substantially increase the accuracy of Y-mer in predicting haplogroups from lower coverage range of data. Human chrY 0.001x sequencing depth appeared to be sufficient for accurate haplogroup prediction in our tests where the validation set was down-sampled to lower coverage (Fig. 3).

Analyses of two ancient DNA data sets from Europe (Saag et al., 2017, 2019; Gretzinger et al., 2022) revealed high accuracy of base haplogroup prediction with M213E model which does not differentiate between recently diverged sub-clades of R1b and I1. The M222NE model that was designed to predict these sub-clades with a training set that combines them with the global range of haplogroup diversity performed less well (accuracy < 0.53) than the case where either the M110W or M213E model’s haplogroup R1 or I1 predictions were further resolved with subclade-specific models, which showed high (> 0.95) accuracy for high confidence (p < 0.05) calls. These results suggest that a two-stage strategy in haplogroup prediction, whereby the base haplogroup is called first and the sub-clade determination is performed separately may be preferable over models in that combine different levels of haplogroup diversity.

Our analyses of Chinese and Estonian NIPS data further confirmed the good performance of Y-mer at low coverage range, 0.001-0.12x, of the male fetus Y chromosome data as we obtained haplogroup frequency profiles similar to those expected from relevant reference data. While these results show that the most common haplogroups in Europe and Asia can be predicted with sufficient accuracy, the drop of accuracy we observe with models that entail haplogroup distinction at sub-clade level further emphasize the limitation of our currently described approach for purposes that require both high confidence and resolution, such as the determination of genetic relatedness. However, in cases where genetic relatedness has been determined independently, e.g. via identity-by-descent or identity-by-state methods (Monroy Kuhn, Jakobsson and Günther, 2018; Popli, Peyrégne and Peter, 2023), Y-mer analyses can be useful for testing (ruling out) the plausibility of patrilineal relatedness, even when acknowledging that a match at a generic Y chromosome haplogroup level cannot constitute a proof of patrilineal relationship.

In conclusion, we present a new k-mer based tool, Y-mer, for predicting Y chromosome haplogroups. We show that Y-mer is able accurately to predict basal chrY haplogroups from ultralow (> 0.001x) coverage data. As such, it is an approach that can be useful in situations where basic, low resolution, information about individual ancestry is required while higher coverage sequencing of the samples is either not possible or practical, e.g. due to costs or inavailability of sufficient quantities of the sample. For ancient DNA studies or forensic case analyses this approach can potentially make more individual samples available for Y chromosomal ancestry analyses, which can be more informative when high coverage/quality data is already available for a subset of the samples. We show that Y-mer performs more accurately when its models are trained on data that match the haplogroup composition of the target group, which highlights the needs for tailored approach in cases where detailed sub-haplogroup level distinctions are required. Besides its possible uses for Y chromosome data, the k-mer based approach described here is potentially extendable also on ancestry analyses of the autosomal genome. Considering the high rate of genetic variation detected in centromeric regions assembled with long read sequences and the low rate of recombination in the pericentromeric regions (Logsdon et al., 2024), the study of k-mers from (peri-)centromeric haplotypes may offer new prospects for autosomal ancestry scanning from low coverage data, similarly, though not identically, to Y-chromosome analyses described here, considering the differences in inheritance of autosomal and Y chromosome DNA. Alternatively, genome-wide scans of population-specific peaks of autosomal k-mer abundance could be screened and used in ancestry mapping of low coverage data. Development of such tools would require larger ancestry diverse reference panels such as graph-based pangenomes that are currently being developed and likely become available in the near future.

Human Y chromosome nomenclatures

Most previous population genetic studies have used for ancestry analyses and Y chromosome haplogroup determination only the X-degenerate regions of Y chromosome that capture less than one fifth of the length of Y chromosome and are deemed to be amenable for short-read sequence mapping (Francalacci et al., 2013; Mendez et al., 2013; Poznik et al., 2013; Wei et al., 2013; Hallast et al., 2015; Karmin et al., 2015; The 1000 Genomes Project Consortium et al., 2016). Human Y chromosome haplogroups are defined by the combined presence of unique sets of allele variants, typically SNVs, that co-occur in individuals who share patrilineal ancestry with the exclusion of other individuals examined. Numerous attempts have been made to update the Y chromosome haplogroup nomenclature since its establishment on the basis of 245 markers in a globally representative sample (Consortium, 2002). This alphanumeric nomenclature system continued to be updated from 2005 to 2020 by a large group of citizen scientists (https://isogg.org/tree/). As a result of ever-increasing volumes of sequence data some updated ISOGG haplogroup labels exceed 20 characters while many sub-clades of the Y phylogeny remain poorly labeled. To find shorter and more stable alternative solutions van Oven et al (2013) proposed a minimal reference tree and Karmin et al., 2015 a shorter time-depth constrained haplogroup labeling system. The latter system will be used when referring to haplogroup names in this study. In cases where data sources have used parallel systems, such as ISOGG-based labeling, these will be explicitly referred to for clarity.

Selecting chrY haplogroups to model (Tables S6)

In order to explore the potential of using k-mer based metrics to capture Y chromosome variation at different evolutionary time depths we dissected the single nucleotide variants (SNV) based Y chromosome phylogeny of the 1000 Genomes Project data (The 1000 Genomes Project Consortium et al., 2016) at different levels of depth (Table S6). Firstly, at the broader global level (level 1, “WORLD”, or “W"), we distinguish 11 branches older than 20 thousand years that are representative of variation across the world. Next, at level 2 we focus on 22 branches that are informative of regional differences across Europe (“EUROPE”, “E“). At level 3 our focus on 23 branches includes sub-clades that are younger than 5 thousand years and these clades are particularly informative for the study of Y chromosome variation in Northeast-Europe (“NE_EUROPE”, “NE“). We assembled three lists of individuals, 110Y, 213Y, and 222Y representing these three layers of phylogeny, respectively, for k-mer based model training and testing including WGS data from 1000 Genomes Project (The 1000 Genomes Project Consortium et al., 2015), EGV (Mitt et al., 2017), and HGDP (Bergström et al., 2020).

Selecting chrY specific k-mers for the haplogroup prediction models

Four different k-mer selection approaches were tested in this study. In the simplified model we chose a single k-mer from each of the three chrY repeat arrays (Table S1) from HG002 assembly (Rhie et al., 2023). The chrY k-mer counts were normalized by k-mer counts of a LINE element, dispersed widely across the genome. The second selection (1Y, HG002) has all k-mers in DYZ1, DYZ2, DYZ3 and DYZ19 regions not presented in other chromosomes in the CHM13v2.0 assembly in total of 618,393 k-mers (Table S2).

Third a selection of 21Y chromosomes of different haplogroup assemblies from Hallast et al. (2023) was used for k-mer selection. All chrY sequences of the 21 chrY assemblies were merged from which a list of all possible 25-mers was extracted. From this list all 25-mers that were found to occur 6 or more times in a k-mer list of female WGS-s (15 from 1000G and 15 from EGC, Table S8) were excluded. The comparisons of k-mer lists were performed with glistcompare function of the GenomeTester4 package (Kaplinski, Lepamets and Remm, 2015). In the end, 12,003,307 k-mers, were retained for further analyses (Tables S5).

Last, fourth selection (110W, 213E, 222NE) is based on unassembled short-read sequencing data from 110, 213 and 222 individuals (Table S7) from across the world including 1000 Genomes Project data (The 1000 Genomes Project Consortium et al., 2015) and Estonian Genome Centre high coverage whole genome sequence data (Mitt et al., 2017). For each model, the union list of different haplogroup intersection lists was merged after the exclusion of k-mers found in the female list (Table S8). In total, approximately 14 million k-mers were retained for further analyses.

k-mer manipulations

GenomeTester4 allows the user to create k-mer lists from fasta and fastq files with its glistmaker function. The k-mer lists can be compared with glistcompare function and the frequency of each k-mer can be obtained with the glistquery function.

Selecting haplogroup-specific k-mers for models applied on three training sets

In each haplogroup prediction model, based on different chrY data sources and k-mer lists (Fig. 2), we used as a training set one of three lists of individuals: WORLD training set, corresponding to 110Y set of Y chromosomes from chrY data sources; EUROPE training set corresponding to 213Y and NE_EUROPE training set corresponding to 222Y. For each model, we queried Y chromosome specific k-mer frequencies in all individuals (Table S7) in a given training set, normalized with chrY median sequencing depth per individual using 5,071,089 k-mers extracted from GRCh38 (Sauk et al., 2018).

As the next step, we determined for each haplogroup distinctive k-mers that were either more or less frequent in individuals of the given haplogroup compared to the pool of individuals from all other haplogroups for that model. Mann-Whitney test, implemented in the MWS.R script, was used to determine 10,000 k-mers with the highest and 10,000 k-mers with the lowest test statistic values for the given haplogroup. We excluded k-mers used in sequence depth estimation and explored models with up to 100,000 k-mers per haplogroup (Table S5).

Validation set

To test the accuracy of the haplogroup prediction models we generated a validation set V110 composed of 110 individuals chosen from high and low coverage data sources. Only individuals that had not been used in the training set were considered for the inclusion in V110. The V110 set (Table S9) mirrors in its haplogroup composition the global training set 110W (Table S6) – each of the eleven basal haplogroups are represented by 10 individuals. The V110 set is predominantly composed of individuals of the 1000 Genomes Project data (The 1000 Genomes Project Consortium et al., 2015). As AB HG-s samples have limited representation, we added AB individuals from Human Genome Diversity Panel (HGDP) data (Bergström et al., 2020). Similarly, to increase the representation of haplogroup N variation, individuals from the EGC dataset (Mitt et al., 2017) were added.

Ethics approval and consent to participate

Human subjects: Variants from the individuals used to create the model are part of the gnomAD database (https://gnomad.broadinstitute.org/). aDNA and NIPS samples were originally published elsewhere.

Availability of Data and Materials

The workflow outlining the steps involving the R scripts (MWS.R, MODEL.R, DILUTER.R, and PREDICTER.R) is presented in Figure 2. These R scripts, along with relevant documentation on their input and output formats, are available at https://github.com/bioinfo-ut/Y-mer/. Running MWS.R and MODEL.R requires 35GB of storage space per individual. Ready-to-use model files, including relevant k-mer dictionaries, PREDICTER.R, and instructions, are available at https://bioinfo.ut.ee/randomtandem/mudelid/.

Funding

This work was funded by the Fonds Wetenschappelijk Onderzoek Grants G0A4521N (T.K.), G050822N (T.K.), and S003422N (T.K.); the Estonian Research Council Grant PRG1076 (A.S., K.K.); the Horizon Europe NESTOR project Grant 101120075 (A.S., K.K.); and partly by the Estonian Business and Innovation Agency Grant RE.5.04.23-0214 (SCANS) (T.P., L.K., M.R.), as well as by Estonian Research Council Grant TEM-TA35, 2021-2027.1.01.24-0627 (T.P., L.K., M.R.).

Authors' Contributions

T.P. and T.K. initiated the project. T.P., T.K., L.K., and M.R. designed the project until the final concept. K.K. prepared the Estonian NIPS samples. M.M. wrote the R scripts and built the statistical framework. K.M. named the project. T.P. and K.M. managed and processed the data. T.K., T.P., and M.M. wrote the first draft of the manuscrip. All authors contributed to editing the manuscript. T.K. and M.R. coordinated and equally funded the study, sharing last authorship.

Acknowledgements

We thank Lehti Saag for sharing various Estonian aDNA samples, Richard Villems for discussions on human genomes, Mart Kals for valuable advice on data presentation, and Reidar Andreson for setting up the web-based tool for Y-mer. Data analysis was performed using the facilities of the High-Performance Computing Center at the University of Tartu.

Bergström, A. et al. (2020) ‘Insights into human genetic variation and population history from 929 diverse genomes’, Science, 367(6484), p. eaay5012. Available at: https://doi.org/10.1126/science.aay5012.
Consortium, T.Y.C. (2002) ‘A Nomenclature System for the Tree of Human Y-Chromosomal Binary Haplogroups’, Genome Research, 12(2), pp. 339–348. Available at: https://doi.org/10.1101/gr.217602.
Damgaard, P.D.B. et al. (2018) ‘137 ancient human genomes from across the Eurasian steppes’, Nature, 557(7705), pp. 369–374. Available at: https://doi.org/10.1038/s41586-018-0094-2.
Esteller-Cucala, P. et al. (2023) ‘Y chromosome sequence and epigenomic reconstruction across human populations’, Communications Biology, 6(1), p. 623. Available at: https://doi.org/10.1038/s42003-023-05004-9.
Francalacci, P. et al. (2013) ‘Low-Pass DNA Sequencing of 1200 Sardinians Reconstructs European Y-Chromosome Phylogeny’, Science, 341(6145), pp. 565–569. Available at: https://doi.org/10.1126/science.1237947.
Gretzinger, J. et al. (2022) ‘The Anglo-Saxon migration and the formation of the early English gene pool’, Nature, 610(7930), pp. 112–119. Available at: https://doi.org/10.1038/s41586-022-05247-2.
Hallast, P. et al. (2015) ‘The Y-Chromosome Tree Bursts into Leaf: 13,000 High-Confidence SNPs Covering the Majority of Known Clades’, Molecular Biology and Evolution, 32(3), pp. 661–673. Available at: https://doi.org/10.1093/molbev/msu327.
Hallast, P. et al. (2023) ‘Assembly of 43 human Y chromosomes reveals extensive complexity and variation’, Nature, 621(7978), pp. 355–364. Available at: https://doi.org/10.1038/s41586-023-06425-6.
Hui, R. et al. (2024) ‘Genetic history of Cambridgeshire before and after the Black Death’, Science Advances, 10(3), p. eadi5903. Available at: https://doi.org/10.1126/sciadv.adi5903.
Kaplinski, L., Lepamets, M. and Remm, M. (2015) ‘GenomeTester4: a toolkit for performing basic set operations - union, intersection and complement on k-mer lists’, GigaScience, 4(1), p. 58. Available at: https://doi.org/10.1186/s13742-015-0097-y.
Karmin, M. et al. (2015) ‘A recent bottleneck of Y chromosome diversity coincides with a global change in culture’, Genome Research, 25(4), pp. 459–466. Available at: https://doi.org/10.1101/gr.186684.114.
Li, J. et al. (2023) ‘Comprehensive insights into the genetic background of Chinese populations using Y chromosome markers’, Royal Society Open Science, 10(9), p. 230814. Available at: https://doi.org/10.1098/rsos.230814.
Logsdon, G.A. et al. (2024) ‘The variation and evolution of complete human centromeres’, Nature, 629(8010), pp. 136–145. Available at: https://doi.org/10.1038/s41586-024-07278-3.
Mendez, F.L. et al. (2013) ‘An African American Paternal Lineage Adds an Extremely Ancient Root to the Human Y Chromosome Phylogenetic Tree’, The American Journal of Human Genetics, 92(3), pp. 454–459. Available at: https://doi.org/10.1016/j.ajhg.2013.02.002.
Mitt, M. et al. (2017) ‘Improved imputation accuracy of rare and low-frequency variants using population-specific high-coverage WGS-based imputation reference panel’, European journal of human genetics: EJHG, 25(7), pp. 869–876. Available at: https://doi.org/10.1038/ejhg.2017.51.
Monroy Kuhn, J.M., Jakobsson, M. and Günther, T. (2018) ‘Estimating genetic kin relationships in prehistoric populations’, PloS One, 13(4), p. e0195491. Available at: https://doi.org/10.1371/journal.pone.0195491.
Popli, D., Peyrégne, S. and Peter, B.M. (2023) ‘KIN: a method to infer relatedness from low-coverage ancient DNA’, Genome Biology, 24(1), p. 10. Available at: https://doi.org/10.1186/s13059-023-02847-7.
Poznik, G.D. et al. (2013) ‘Sequencing Y Chromosomes Resolves Discrepancy in Time to Common Ancestor of Males Versus Females’, Science, 341(6145), pp. 562–565. Available at: https://doi.org/10.1126/science.1237619.
Rhie, A. et al. (2023) ‘The complete sequence of a human Y chromosome’, Nature, 621(7978), pp. 344–354. Available at: https://doi.org/10.1038/s41586-023-06457-y.
Saag, Lehti et al. (2017) ‘Extensive Farming in Estonia Started through a Sex-Biased Migration from the Steppe’, Current Biology, 27(14), pp. 2185-2193.e6. Available at: https://doi.org/10.1016/j.cub.2017.06.022.
Saag, Lehti et al. (2019) ‘The Arrival of Siberian Ancestry Connecting the Eastern Baltic to Uralic Speakers further East’, Current biology: CB, 29(10), pp. 1701-1711.e16. Available at: https://doi.org/10.1016/j.cub.2019.04.026.
Sauk, M. et al. (2018) ‘NIPTmer: rapid k-mer-based software package for detection of fetal aneuploidies’, Scientific Reports, 8(1), p. 5616. Available at: https://doi.org/10.1038/s41598-018-23589-8.
Skov, L., The Danish Pan Genome Consortium and Schierup, M.H. (2017) ‘Analysis of 62 hybrid assembled human Y chromosomes exposes rapid structural changes and high rates of gene conversion’, PLOS Genetics. Edited by K.D. Makova, 13(8), p. e1006834. Available at: https://doi.org/10.1371/journal.pgen.1006834.
The 1000 Genomes Project Consortium et al. (2015) ‘A global reference for human genetic variation’, Nature, 526(7571), pp. 68–74. Available at: https://doi.org/10.1038/nature15393.
The 1000 Genomes Project Consortium et al. (2016) ‘Punctuated bursts in human male demography inferred from 1,244 worldwide Y-chromosome sequences’, Nature Genetics, 48(6), pp. 593–599. Available at: https://doi.org/10.1038/ng.3559.
Wei, W. et al. (2013) ‘A comparison of Y-chromosomal lineage dating using either resequencing or Y-SNP plus Y-STR genotyping’, Forensic Science International: Genetics, 7(6), pp. 568–572. Available at: https://doi.org/10.1016/j.fsigen.2013.03.014.
Xu, H. et al. (2018) ‘Informative priors on fetal fraction increase power of the noninvasive prenatal screen’, Genetics in Medicine, 20(8), pp. 817–824. Available at: https://doi.org/10.1038/gim.2017.186.
Žilina, O. et al. (2019) ‘Creating basis for introducing non‐invasive prenatal testing in the Estonian public health setting’, Prenatal Diagnosis, 39(13), pp. 1262–1268. Available at: https://doi.org/10.1002/pd.5578.

No competing interests reported.

Download PDF

Reviewers agreed at journal
12 Nov, 2024
Reviewers agreed at journal
11 Nov, 2024
Reviewers invited by journal
08 Nov, 2024
Editor assigned by journal
12 Sep, 2024
Submission checks completed at journal
06 Sep, 2024
First submitted to journal
06 Sep, 2024

You are reading this latest preprint version

Y-mer: A k-mer based method for determining human Y chromosome haplogroups from ultra-low sequencing depth data

Status:

Version 1

Abstract

Figures

INTRODUCTION

RESULTS

Complex distance-based models for predicting Y chromosome haplogroups in multi-population data

Model validation

Model testing on ancient DNA data

Model testing on non-invasive prenatal screening (NIPS) data

DISCUSSION

METHODS

Human Y chromosome nomenclatures

Selecting chrY haplogroups to model (Tables S6)

Selecting chrY specific k-mers for the haplogroup prediction models

k-mer manipulations

Selecting haplogroup-specific k-mers for models applied on three training sets

Validation set

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1