YHap: software for probabilistic assignment of Y haplogroups from population re-sequencing data

Y haplogroup analyses are an important component of genealogical reconstruction, population genetic analyses, medical genetics and forensics. These fields are increasingly moving towards use of low-coverage, high throughput sequencing. However, there is as yet no software available for using sequence data to assign Y haplogroup groups probabilistically, such that the posterior probability of assignment fully reflects the information present in the data, and borrows information across all samples sequenced from a population. YHap addresses this problem.


INTRODUCTION
The non-recombining portion of haploid chromosome Y is passed intact from father to son with a mutation rate several times greater than autosomes (Xue et al. 2009). As such, patterns of variation in Y are widely used to uncover historical patterns of human migration; are important in genealogical reconstruction and have application in forensic analyses.
The Y Chromosome Consortium (YCC) published revised Y-chromosome DNA haplogroup tree in 2008, consisting of approximately 600 markers, which can be used to characterize 20 major global haplogroups, labeled A-T, as well as sub-classification into a total of 311 haplogroups at the finest level of resolution. Different major haplogroups have been found at high frequencies in different geographical regions, for example the E clade in Africa, and the O clade in Eastern Asia. Particular fine-level haplogroups are found in multiple locations, such as R1a in Eastern Europe, South Asia and Central Asia, indicating migration of R1a from Eurasian Steppes to the new world. The C3 haplogroup, found at high frequency throughout Asia is commonly interpreted as genealogical remnants of the empire of Genghis Khan 2 (Zerjal et al. 2003).
Y haplogroup assignment has traditionally been carried out by targeted genotyping using a combination of short tandem repeat typing, multi-plex PCR and minisequencing (Sanchez et al. 2003, Zerjal et al. 2002, often using a hierarchical strategy in order to first refine the major haplogroup, and subsequently genotype markers within that haplogroup which illuminate finer levels of resolution. Such a procedure requires substantial amount of wet-lab analysis, requires stringent replication and quality control to eliminate errors which can arise due to the limited amount of information collected at each step. More recently, personal genetic companies have included specific Y chromosome markers on custom genotyping arrays (Turner et al. 2008). Nevertheless, the resolution available from genotyping arrays is limited by markers included on the chip.
Very high coverage high throughput sequencing has the potential to capture all single nucleotide and insertion/deletion variation, and as such provide near-perfect assignment of individuals to Y haplogroups. As high coverage sequencing of large population samples remains expensive, low coverage population sequencing, in which each individual is sequenced at less than 2x haploid coverage is an attractive alternative, but this will not capture all individual-level variation. However, given the sharing of haplogroups within an ethnically homogenous population, it should be possible to borrow information across individuals within a population in order to improve haplogroup assignment.
In this application note, we present the YHap tool, which has been primarily designed for assigning haplogroups to low-coverage population re-sequencing data. YHap borrows information across all samples to assign samples to haplogroups probabilistically, thus providing an accurate representation of the inference which can be made from the data collected. YHap is a complete solution and can also be applied to high-coverage sequence data, as well as data from genotyping arrays.

METHODS
We use the set of haplogroups and mutations defined in (Karafet et al. 2008). We map the forward and reverse primers described in this manuscript to identify the reported strand of the variation in the GRCh37 reference. After strand correction, we identify whether the mutant allele is the equal to the alternative or reference allele, so that we can subsequently work in reference/alternative allele space on the forward strand, consistent with conventional genotype calling schema. Next, we map each mutation to its position on the pre-defined Y phylogenetic tree T. Finally we create a haplogroup matrix H of size Nref*M where M is the total number of nodes in T (including leaf and internal nodes) and Nref is the number of pre-defined Y markers. Each entry Hij = {Hijg} is a probability distribution vector expressing the probability of a sampled individual from the clade below node j carries allele g (in this case either the reference or alternate allele).
At leaf nodes, this probability vector is either {0,1} or {1,0}, and at internal nodes, it is the proportion of descendant leaf nodes with reference or alternate alleles, respectively.
To assign a sequenced individual to a specific haplogroup, we obtained genotype likelihoods at each putative variant site (inclusive of all markers in H) from chrY VCF files of 1000 genome project. This results in a matrix S of size L*N where L is the number of sequenced samples, N is the number of putative variants, and Gij = {Gijg} is a vector of genotype likelihoods. We then generated an augmented H* matrix by adding in extra sites in G but not H with probability vector H*ij = {0.5,0.5}. The pipeline is similar for genotype data, except that the genotype likelihoods are taken to be either 1, if , or 0 otherwise. We can calculate the assignment of each individual using (1) Where and = Gijg We can then calculate the posterior probability of each haplogroup amongst a set of haplogroups, where prior haplogroup probability distribution P(H * .l) is set to the uniform distribution, By restricting the set of haplogroups considered in equation (2), YHap can be customized to either only assign to within the major haplogroups (A through to T), or all possible haplogroups at the finest level of classification.
While this model is sufficient for assigning Y haplogroups individually, it does not capture shared information between sequenced samples adequately, particularly for low coverage sequencing. Given that a population sample will share individuals from the same haplogroup, and while none of these individuals may have enough depth at informative Y haplogroup markers, there is enough information across the pooled reads from all samples from the same haplogroup. However, we do not know a-priori which samples can be pooled as coming from the same haplogroup.
In order to pool information between samples, we treat the allele probability distribution {H*ijg} at markers present in the sequence data but not present as haplogroup markers, as parameters in our model. We update these parameters using expectation maximization, in which we first calculate the posterior probability assigning each sample j to each haplogroup l using equation 2, and then update the {H*ilg} to reflect the average of genotypes assigned to haplogroup l at position i, weighted by this posterior probability of assignment. In this way, the model learns which alleles are characteristic of the pre-defined haplogroups, and is thus able to more accurately assign individuals which may not have good coverage at those sites, but which show similarity to other individuals across the Y chromosome. The probabilities P(H * .l) are also updated at each step to reflect the proportion of haplogroups assigned in the population.

RESULTS
We applied YHap to low-coverage sequencing data generated in the pilot phase of the 1000 genome consortium which were also part of the Hapmap project, including 19 YRI, 16 JPT, 21 CEU and 9 CHB samples (Abecasis et al. 2010). Major Y haplogroups have been previously assigned to these samples as part of the Hapmap project. The average sequencing depth of these samples is 1.67X as described in 1000 genome Y chromosome analysis report. Compared to haplogroups previously obtained from the Hapmap project (Altshuler et al. 2010), YHap showed perfect assignment accuracy (Table 1). We also used YHap on the Hapmap combined phase 1,2,3 Y genotype data and obtained the same assignments previously reported with this data.
In order to investigate the ability of YHap to assign finer-level haplogroups we compared YHap results obtained at complete resolution (i.e. considering all haplogroup leaf nodes on the pre-defined Y phylogenetic tree) on both Hapmap genotype data and also 1000 genomes low coverage sequence data (Suppl. Table 1). We see that there is complete concordance at the major haplogroup level, and there is increasing uncertainty in assignment as the resolution of assignment increases, particularly using dense genotype data. We also observe that accuracy remains high amongst those assignments which YHap assigns high confidence.
Finally, in order to investigate the relationshis between sequencing depth and assignment accuracy, we randomly downsampled the original bam files from 1000 genome to 0.6X. For simplicity, we chose JPT to run the test. We see that downsampling increases the uncertainty of assignment (Suppl. Table 2), but YHap accuracy remains high amongst those assignments which are made with high posterior probability. This demonstrates that as the underlying amount of information decreases, YHap is still able to extract inference and accurately represent the uncertainty of this inference.
The total complexity for the whole procedure is O(N 2 T), conventionally, when using defaut settings, it will take almost 10min to locate 10~20 individuals and approximately 200 Mb memory. Numbers indicate the percentage of samples which had an assignment with a given lower bound on certainty from 1000 genomes sequence data (horizontal) and Hapmap genotype data (vertical), and different levels. Level 1 indicates main haplogroup assignment (A-T), while levels 2-5 indicate increasing levels of precision. Numbers in brackets indicate the concordance rate between predictions made by sequence and genotype data within each bin.