Divisive K-Means clustering to partition the data

K-Means clustering is a well-known technique for uniquely assigning samples to one of K clusters based on the similarity of some real-valued attribute. Assignments are performed so as to minimize the within-cluster sum of squares, ${\Sigma }_j^{2N}\parallel H_j-{\mu }_{C_j}{\parallel }^2,$, where μ_c is the cluster mean for class c and C_j denotes the cluster membership of haplotype j. This objective function can be minimized by iteratively assigning class labels to each observation based on which cluster mean is closest and then recalculating the cluster means. When the distance measure is Euclidean, convergence is guaranteed²². We implemented this routine, allowing up to ten iterations (or stopping when class labels stabilize). Ten iterations are often not sufficient for convergence but we only require a rough partitioning. We used K-Means with K=2 to partition clusters of haplotypes recursively. Each cluster is split into two smaller clusters. The recursion is terminated when the cluster size is <M. Upon termination, we “top up” the cluster at the leaf of the bifurcating tree with haplotypes from the closest leaf. This last step ensures each haplotype is compared to M-1 other haplotypes.

Formally, we wish to find a vector of labels C = {C₁…C_2N} where C_i tells us the cluster which haplotype H_i belongs to. We refer to C_j as the primary cluster of H_j. We also store a dictionary D of sets, which we call secondary clusters. Every set in D has exactly M<<2N elements. It is possible for a haplotype to be in more than one secondary cluster, but haplotypes will only ever have one primary cluster (determined by C_j). We assume to that we have the functions available; K-Means(H) and euc(H_j, μ_c).The first performs the K-means clustering routine with K=2 on H and returns a vector of class labels, the second which calculates the Euclidean distance between H_j and μ_c. Our algorithm is described using pseudo-code in the Supplementary Note.

On each iteration of the SHAPEIT3 algorithm, we find the primary cluster labels C and build the dictionary of secondary clusters D. For each individual i we choose K conditioning haplotypes from H* = D[C_2i-1] D[C_2i] based on the minimum Hamming distances but excluding individuals estimated to be IBD2 within the region ⁸. We set the cluster size to M=4,000 for all the experiments described in this paper. We note that the final set of K conditioning haplotypes are chosen based on minimum Hamming distance, whereas the K-means routine uses Euclidean distance. Hamming distance is appropriate for pairs of binary vectors, we use this where possible since such distances can be calculated very rapidly via lookup tables. The K-Means routine needs to use a distance measure appropriate for real numbers, since the means of the clusters are unlikely to be 0 or 1 rather something on the interval [0,1], hence Euclidean distance was used. We suspect that this part of the algorithm can be improved via careful choice of the SNPs used to calculate Euclidean distances between haplotypes. For example, by focusing preferentially on rare SNPs²¹.

We further decrease the computational cost of our routine by thinning the SNPs that K-means is performed on. For each K-means clustering performed (each leaf of the recursion), Euclidean distances are only calculated on SNPs l={8k+o : k 0… L/8} where o {0,1,2,3,4,5,7} and is sampled randomly within each leaf. In words, we only use every eighth SNP (starting from a random offset) in our clustering routine. So in practice 4NL/8 differences (2 centroids, 2N haplotypes, L/8 SNPs) are calculated per iteration of K-means when clustering N samples, we include this description for completeness but remove the 8 in the denominator henceforth for clarity. The initial iteration of K-Means clustering requires 4NL calculations for calculating the distance between the N haplotypes and two clusters for L SNPs. As a rough approximation, we assume that K-Means divides the data into two clusters of equal size. Consequently, after the d^th recursion, we have 2^d clusters of size N/2^d. The recursion continues until a cluster of size <M is reached. This will take log₂(N/M) recursions. Hence, we can derive the computational complexity of the clustering routine quite simply:

Since L and M are constants our clustering routine has O(NlogN) complexity.

Modified MCMC sampling routine

We implement two modifications to the MCMC sampling routine of SHAPEIT2. Firstly, when a Hamming distance of zero is found between a haplotype and at least one conditioning haplotype, we do not perform HMM calculations within that window on the current iteration. In such cases, it is unlikely that a different haplotype will be generated in preference to this perfect match. In other words, when there is evidence that locally that the algorithm has converged, then we stop updating that individual. We refer to this as adaptive algorithm termination. Therefore, we simply carry forward the haplotype generated by the previous iteration. As shown in Supplementary Figure 2, this can be expected to happen quite frequently as N increases leading to considerable savings in computation time.

The second modification is an approximation to the MCMC sampling routine, which offers significant additional gains in speed when using multiple threads. In SHAPEIT2, the parallelization scheme updates the phase of T samples in parallel using T threads conditioning on the N-T other samples already processed. The problem with this approach is a given thread becomes idle as soon as it has finished processing its assigned sample before the other threads, and therefore has to wait for the other threads to finish. Because of this synchronization that has to be repeated many times, the parallelization efficiency of this scheme decreases with the number of threads. In SHAPEIT3, we introduce a new parallelization scheme based on another approximation of the MCMC process. This scheme takes advantage of the fact that SHAPEIT updates the phase of samples within windows of usually few Mb: it splits the N samples into W overlapping windows prior to any MCMC iteration and distributes the N × W jobs to be done into T distinct queues, one for each thread. Then, a thread performs the window-sample pairs in its queue independently of the other queue and therefore does not need any particular synchronization, which results in better usage of all CPU cores. Once all the jobs are done, the haplotype for each sample are updated and we move on the next MCMC iteration. Of note, this scheme implicitly assumes a MCMC scheme in which the phase of all samples is updated in parallel at the same time. In other words, we move from a Markov chain state to the next by updating all samples simultaneously which differs to the standard MCMC scheme in which we move from as state to the other by updating only one sample. We refer to this as completely parallel updating. We assessed the performance of this approach by running SHAPEIT3 without the new clustering routine on the 10,072 sample subset of the UK Biobank dataset and compared it to SHAPEIT2 run on the same dataset. The median switch error rate was 1.1% for both methods, indicating little difference in performance when using this new MCMC routine (see Table 1). Performance for increasing numbers of threads is shown in Supplementary Table 4.

Creating a validation data set for assessing accuracy

SHAPEIT3 was initially developed and tested on the UK-BiLEVE dataset, which contained a large number of apparent sibling pairs but no identifiable mother-father-child trios. Regions in which siblings inherit exactly one chromosome from the same parent, (IBD1) can be phased using simple Mendelian rules. So we developed a novel method for inferring the identity by descent (IBD) status of siblings, with the aim of constructing a set of ‘known’ haplotypes against which we could compare accuracy of different methods.

We identified 384 pairs of individuals with a kinship coefficient >0.35, indicating they were likely first order relatives (sibling-pairs or parent-child duos). Parents will share 100% of the chromosomes IBD1 with their children and siblings will share (on average) their genomes 25% IBD0, 50% IBD1 and 25% IBD2. No pairs of individuals had close to 0% IBD0 suggesting these samples were sibling-pairs, not parent-child²³. By detecting IBD1 regions we can accurately phase loci within those regions using simple Mendel rules. This will generate accurate haplotypes for validation purposes, albeit partially phased ones.

We infer regions of IBD1 sharing by a two-stage approach. Firstly, we applied a simple Hidden Markov Model (HMM) to the genotypes of each sib-pair (see Supplementary Note), with three unobserved states indicating IBD status. However, because this model does not account for LD (IBD sharing from more distant ancestors) it is vulnerable to inferring spurious stretches of IBD1. This issue has been noted in previous work on phasing siblings with missing parents²⁴. We circumvent it via simple post-hoc filtering of short (<10cM) IBD1/IBD2 segments. Supplementary Figure 6 shows the consequences of this filtering step for the IBD state for a pair of probable siblings across chromosome 20. The filtering step did indeed reduce the SE for all methods (Supplementary Figure 7) and also flattened the distribution of switch error rates with respect to the percentage of phase resolved sites for the corresponding individual (Supplementary Figure 8).

We note that the switch-error rates when using these validation data were higher than the trio phased validation haplotypes created in the UK Biobank experiments. This suggests our sibling-pair validation haplotypes were not as accurate as those produced via trio phasing. We can think of two likely causes. While we filter short stretches of (likely spurious) IBD1, such a spurious segment may also flank a longer (true) IBD1 segment introducing some errors on the edges of the segment. Secondly, genotyping errors will induce incorrect phase, and there is greater power to detect such errors in the trio setting. Nevertheless, this validation data set provides a reasonable relative comparison of methods and was of great utility when developing SHAPEIT3.

UK BiLEVE was funded by a Medical Research Council (MRC) strategic award to MDT, IPH, DS and LVW (MC_PC_12010). We would like to thank the UK BiLEVE study investigators, array design team, analysis group and steering group. This research has been conducted using the UK Biobank Resource. This article presents independent research funded partially by the National Institute for Health Research (NIHR). The views expressed are those of the author(s) and not necessarily those of the NHS, the NIHR or the Department of Health. J.M acknowledges support from the ERC (Grant no. 617306).