HapCompass: a fast cycle basis algorithm for accurate haplotype assembly of sequence data

Abstract

Genome assembly methods produce haplotype phase ambiguous assemblies due to limitations in current sequencing technologies. Determining the haplotype phase of an individual is computationally challenging and experimentally expensive. However, haplotype phase information is crucial in many bioinformatics workflows such as genetic association studies and genomic imputation. Current computational methods of determining haplotype phase from sequence data--known as haplotype assembly--have difficulties producing accurate results for large (1000 genomes-type) data or operate on restricted optimizations that are unrealistic considering modern high-throughput sequencing technologies. We present a novel algorithm, HapCompass, for haplotype assembly of densely sequenced human genome data. The HapCompass algorithm operates on a graph where single nucleotide polymorphisms (SNPs) are nodes and edges are defined by sequence reads and viewed as supporting evidence of co-occurring SNP alleles in a haplotype. In our graph model, haplotype phasings correspond to spanning trees. We define the minimum weighted edge removal optimization on this graph and develop an algorithm based on cycle basis local optimizations for resolving conflicting evidence. We then estimate the amount of sequencing required to produce a complete haplotype assembly of a chromosome. Using these estimates together with metrics borrowed from genome assembly and haplotype phasing, we compare the accuracy of HapCompass, the Genome Analysis ToolKit, and HapCut for 1000 Genomes Project and simulated data. We show that HapCompass performs significantly better for a variety of data and metrics. HapCompass is freely available for download (www.brown.edu/Research/Istrail_Lab/).

FIG 1.

The SNP-fragment matrix M is shown on the left containing four fragments and four SNPs. Each pairwise phasing relationship defined by the fragments is represented on the edges of the compass graph on the right. A positive edge indicates there is more evidence in the fragments for the phasing while a negative edge indicates evidence for the opposite phasing.

FIG 2.

A compass graph G_C is shown on the left with two conflicting cycles. One edge removal (s₂, s₃) makes G_C happy by removing two conflicting cycles in one step. All spanning trees (ST) of the happy G_C correspond to the same phasing, but only two are shown in the lower right corner.

FIG 3.

In this simulation study, reads of size 100 bp were simulated on chromosome 22 of the 1000 genomes CEU individual NA12878. Mate pair lengths were sampled at random from one of four normal distribution with means [10 kb, 50 kb, 100 kb, 250 kb] and standard deviations [1 kb, 5 kb, 10 kb, 25 kb]. With these parameters for sequencing, we require about 10 million reads to connect most of the SNPs of G_C.

FIG 4.

We fit a linear least squares regression line to the FMPR measurement for algorithms Genome Analysis ToolKit (GATK), HapCut, and HapCompass on chromosome 22 of the 1000 genomes data for the NA12878 individual (top) for the 1000 genomes data only (Table 1) and (bottom) for the 1000 genomes data supplemented with 10 million simulated reads of length 100 and sequence base error rate of 0.05 (Table 2).