Shape-IT: new rapid and accurate algorithm for haplotype inference

Abstract

Background: We have developed a new computational algorithm, Shape-IT, to infer haplotypes under the genetic model of coalescence with recombination developed by Stephens et al in Phase v2.1. It runs much faster than Phase v2.1 while exhibiting the same accuracy. The major algorithmic improvements rely on the use of binary trees to represent the sets of candidate haplotypes for each individual. These binary tree representations: (1) speed up the computations of posterior probabilities of the haplotypes by avoiding the redundant operations made in Phase v2.1, and (2) overcome the exponential aspect of the haplotypes inference problem by the smart exploration of the most plausible pathways (ie. haplotypes) in the binary trees.

Results: Our results show that Shape-IT is several orders of magnitude faster than Phase v2.1 while being as accurate. For instance, Shape-IT runs 50 times faster than Phase v2.1 to compute the haplotypes of 200 subjects on 6,000 segments of 50 SNPs extracted from a standard Illumina 300 K chip (13 days instead of 630 days). We also compared Shape-IT with other widely used software, Gerbil, PL-EM, Fastphase, 2SNP, and Ishape in various tests: Shape-IT and Phase v2.1 were the most accurate in all cases, followed by Ishape and Fastphase. As a matter of speed, Shape-IT was faster than Ishape and Fastphase for datasets smaller than 100 SNPs, but Fastphase became faster -but still less accurate- to infer haplotypes on larger SNP datasets.

Conclusion: Shape-IT deserves to be extensively used for regular haplotype inference but also in the context of the new high-throughput genotyping chips since it permits to fit the genetic model of Phase v2.1 on large datasets. This new algorithm based on tree representations could be used in other HMM-based haplotype inference software and may apply more largely to other fields using HMM.

Figure 1

Schematic representation of a sample of n genotypes. In this example, the space of possible haplotypes S_i for individual i contains 4 haplotype pairs with 8 distinct haplotypes. The possible phases between heterozygous markers are shown in bold.

Figure 2

Representation of the execution trellis of the hidden Markov model used to compute the probability of a haplotype. The haplotypes h₁,..., h_2n-2 denote the previously sampled haplotypes which are used to compute the probability of the observed haplotype h. The sets {o₁,..., o_s} and {q₁(k), ..., q_s(k)} correspond respectively to the observed state sequence of haplotype h and to the hidden state sequence of haplotype h_k. The transition probability a_j(k,l) corresponds to the probability of jumping from hidden state q_j(k) of haplotype h_k to hidden state q_j+1(l) of haplotype h_l, and the emission probability b_j(k) corresponds to the probability of observing o_j given the hidden state q_j(k). To compute the probability of observing the sequence h = {o₁, ..., o_s} in this HMM, one must sum up the probabilities of observing h over all (2n - 2)^s possible sequences of s hidden states which is done efficiently by the forward algorithm.

Figure 3

Different representations of the space of possible haplotypes pairs S_i. The left panel (A) shows the list representation commonly used by haplotype software such as Phase v2.1. The lower right panel (C) shows the representation used by Shape-IT. White and black circles indicate the phases between the heterozygous SNPs. On this example we use the same genotype G_i described in Figure 1. For iterations as performed by Phase v2.1 (A), the list requires the exploration of 20 nodes (4 haplotype pairs × 5 SNPs). With the complete tree representation (B) 10 nodes need to be explored, and with the incomplete tree representation as performed by Shape-IT (C), only 7 nodes need to be explored. The difference observed between (B) and (C) results from the pruning strategy which avoids the exploration of the nodes with probability ≤ 0.01.

Figure 4

Algorithm 1 to compute the FDSL distribution on the complete haplotype tree.

Figure 5

Algorithm 2 to compute the FDSL distribution on the incomplete haplotype tree.

Figure 6

Accuracy of the different values tested for the threshold f in Shape-IT (grey boxes) compared to Phase v2.1 (black line). This comparison was done on 300 datasets of 50 Tag SNPs called CEU Illumina 50.