Whole-genome haplotype reconstruction using proximity-ligation and shotgun sequencing

Experimental strategy of HaploSeq

In HaploSeq, we first perform proximity-ligation sequencing using the previously established Hi-C protocol³⁴ (Fig. 1a). In this protocol, cross-linked DNA is first digested with a restriction enzyme and then ligated to form artificial fragments. This method therefore can capture DNA fragments from two distant genomic loci that looped together in three-dimensional space in vivo^33–35 (Fig. 1a). Indeed, after shotgun DNA-sequencing of the resulting DNA library, paired-end sequencing reads have ‘insert sizes’ that range from several hundred base pairs to tens of millions of base pairs (Fig. 1a). Thus, the short DNA fragments generated in a Hi-C experiment can yield small haplotype blocks, long fragments ultimately can link these small blocks together (Fig. 1b). With enough sequencing coverage, such an approach has the potential to link variants in discontinuous blocks and assemble every such block into a single haplotype.

One complicating factor is that proximity ligation can capture interactions both in cis within an individual allele and in trans between homologous and nonhomologous chromosomes. Although nonhomologous trans interactions between different chromosomes do not affect phasing, interactions in trans between homologous chromosomes (referred to as h-trans hereafter) might complicate haplotype reconstruction if h-trans interactions were as frequent as cis interactions. Therefore, we set out to determine the relative frequency of h-trans versus cis interactions in proximity-ligation sequencing data. To accomplish this, we used a hybrid mouse embryonic stem (ES) cell line derived from a cross between two inbred homozygous strains (Mus musculus castaneous (CAST) and 129S4/SvJae (J129)), for which the parental inbred whole-genome sequences (WGS) were publicly available (Online Methods). As a result, the maternal and paternal haplotypes within this cell line are known a priori, and the frequency of interactions between alleles can then be explicitly tested. We performed a Hi-C experiment and generated >620 million usable 75-bp paired-end reads from these hybrid ES cells, corresponding to 30× coverage of the genome (Online Methods).

To determine the extent of intrahaplotype (cis) versus interhaplotype (h-trans) interactions, we used the prior haplotype information to distinguish reads from CAST and J129 alleles. We first visually checked the pattern of interactions between every allele, finding that the CAST and J129 alleles for each chromosome were largely self-interacting and distinct (Fig. 2a). Such a pattern has been previously observed in Hi-C studies and is analogous to the long-established concept that chromosomes occupy distinct, self-associated territories, known as “chromosome territories,” within the interphase nucleus^34,35. However, previous Hi-C studies did not distinguish whether the two alleles for a given chromosome also occupy distinct, individual, chromosome territories^34,35. Our results demonstrate that each allele occupies a distinct chromosome territory, with h-trans interactions occurring at a frequency of <2% of the total reads originating and ending on the two homologous chromosomes (Fig. 2b). These results indicate that the vast majority of Hi-C reads truly originated from cis interactions, which is essential for using them for haplotyping, as described below. In addition, the probability of a DNA read being in h-trans versus in cis appears to increase as a function of the insert size between the read pairs (Fig. 2c and Supplementary Fig. 2). Because of this trend, we capped the maximum insert size of Hi-C reads at 30 Mb to reduce the overall number of h-trans interactions to ~0.6% (Fig. 2d). Currently, we cannot determine if these rare h-trans interactions are due to noise in the data or to biological phenomena, such as homologous pairing of chromosomes³⁶. Regardless, these observations indicate that h-trans interactions are rare, a prerequisite for HaploSeq analysis to succeed.

Accurate chromosome-spanning haplotypes in mouse ES cells

Rare h-trans interacting reads and phenomena such as sequencing errors at the variant locations can cause erroneous connections between homologous chromosomes and complicate the reconstruction of haplotypes. To overcome these problems, we incorporated HapCUT²⁴ software into HaploSeq analysis to probabilistically predict haplotypes. HapCUT constructs a graph in which nodes correspond to heterozygous variants and edges correspond to overlapping sequence fragments that may link the variants. This graph might contain spurious edges resulting from sequencing errors or h-trans interactions. HapCUT then uses a max-cut algorithm to predict parsimonious solutions that are maximally consistent with the haplotype information provided by the set of input sequencing reads (Fig. 3a). Because proximity ligation generates larger graphs than conventional genome sequencing or mate-pair sequencing, we modified HapCUT to balance computing time and number of iterations, so that the haplotypes can be predicted with reasonable speed and high accuracy (Online Methods and Supplementary Data).

To test the ability of HapCUT to generate haplotype blocks from proximity-ligation and sequencing data, we again used the CAST×J129 mouse ES cell Hi-C data. In this instance, we did not distinguish a priori to which allele a sequencing read belonged. Instead, we allowed HapCUT to reconstruct de novo haplotype blocks of the heterozygous variants. We used the metrics of completeness, resolution and accuracy to assess the performance of our HaploSeq analysis in haplotype reconstruction (Supplementary Fig. 3).

The completeness metric of haplotype phasing is defined as the size in base pairs of the haplotype blocks generated. In general, HapCUT will generate several haplotype blocks of various sizes for each chromosome depending on heterozygous variant connections. To assess completeness, we analyzed the span of the haplotype blocks generated for each chromosome. We observed that each chromosome contains one block with the most heterozygous variants phased (MVP). In addition to the MVP block, a minority of heterozygous variants may be assigned to smaller blocks owing to their inability to be connected with the MVP block. The MVP block in our case spanned >99.9% of the phasable base-pairs for each chromosome (Table 1), demonstrating that HaploSeq analysis using Hi-C data can generate complete, chromosome-spanning haplotypes.

Although completeness is defined as the base-pair span of the MVP block, resolution is defined as the fraction of phased heterozygous variants relative to the total variants spanned in the MVP block (Supplementary Fig. 3). These MVP blocks generated for each chromosome are of high resolution, as we could phase about 95% of the heterozygous variants on any given chromosome (Table 1). As 99.6% of variants are covered by at least one read, the inability to link the 5% of heterozygous variants is primarily due to the inability to link heterozygous variants to the MVP haplotype block. As a result, although the MVP block spans the majority of the chromosome, it has gaps that in total contain ~5% of the heterozygous variants.

To assess the accuracy of the heterozygous variants within the MVP block, we compared the predicted haplotypes generated de novo by HaploSeq analysis with the known haplotypes of the CAST and J129 alleles. We defined accuracy as the fraction of phased heterozygous variants that were correctly phased in the MVP block (Supplementary Fig. 3). Of the variants that were assigned to the MVP haplotype block, we observed >99.5% accuracy in distinguishing between the two known haplotypes (Table 1). Lastly, as we had previously demonstrated that the h-trans interaction probability increases with the genomic distance separating two sequencing reads (Fig. 2c), we incorporated the h-trans interaction probabilities into our HapCUT algorithm (Online Methods). These conditions did not sacrifice the completeness of the haplotypes we generated. Instead, we observed a further improvement in the accuracy of the variants in the MVP block with a modest reduction of the resolution of the variants phased (Supplementary Fig. 4a,b). In summary, these results demonstrate that HaploSeq analysis yields complete, high-resolution and accurate haplotypes for all autosomes.

Previous haplotyping efforts have often combined different shotgun sequencing methods to improve phasing. For instance, whole-genome sequencing has been combined with mate-pair sequencing²⁴. To see if this approach would also improve haplotyping with proximity-ligation data, we simulated 20× coverage DNA sequencing data for conventional paired-end shotgun DNA sequencing (i.e., WGS), mate-pair sequencing, fosmids and proximity ligation. As expected, combining WGS with mate pair or fosmid data resulted in fragmented, incomplete haplotype blocks (Fig. 3b,c). In contrast, performing HaploSeq analysis using proximity ligation in combination with WGS data did not increase the completeness of the haplotypes generated (Fig. 3b) but did marginally improve their resolution (Fig. 3c), suggesting that adding WGS to HaploSeq analysis may be a viable strategy in cases where the resolution of haplotypes must be maximized.

Performance of HaploSeq depends on variant density

A distinct feature of the CAST×J129 ES cell line is the high density of heterozygous variants present throughout the genome. On average, there is a heterozygous variant every 150 bases, which is 7–10 times more frequent than in humans^1,2 (Fig. 4a). As a first test of the feasibility of using HaploSeq to generate haplotypes in human cells, we sub-sampled heterozygous variants in the CAST×J129 data so that the variant density mimics that in human populations. We then tested how lower variant density affects the ability of HaploSeq to reconstruct haplotypes. Although lower variant density did result in fewer usable reads (Supplementary Fig. 5), we still observed complete haplotypes over each chromosome with only a marginal decrease in accuracy (from ~99.6% to ~99.2%) (Supplementary Fig. 6a). However, the MVP block generated using a variant density similar to that observed in the human genome had a lower resolution, with fewer variants phased relative to results obtained when the density of variants was higher. Approximately 32% of heterozygous variants were phased in the MVP block (Supplementary Fig. 6a), instead of 95% in the high-density case (Table 1). In summary, a low density of variants does not affect completeness or accuracy, but does substantially affect the resolution of chromosome-spanning haplotypes by HaploSeq analysis.

HaploSeq analysis of a human individual

To realistically assess the ability of our method to phase haplotypes in humans, we performed HaploSeq on the GM12878 lymphoblastoid cell line. The complete haplotype of this cell line has been inferred by the 1000 Genomes Project from whole-genome sequencing data of a parent-child trio¹⁴ using WGS. We generated >262 million usable, 100-bp paired-end reads corresponding to ~17× coverage (Online Methods). HaploSeq generated chromosome-spanning haplotypes in all chromosomes of the GM12878 cells (Supplementary Fig. 6b and Fig. 4b,c). Of note, previous methods attempting haplotype reconstruction in humans have been unable to reconstruct haplotypes spanning across the highly repetitive centromeric regions of metacentric chromosomes^4,23,25–28. Using HaploSeq, we generated haplotypes that accurately spanned the centromere in all metacentric chromosomes with the exception of chromosome 9, where an erroneous linkage causes switching of haplotype calls at the centromere. Chromosome 9 has both a large 15-Mbp, poorly mapped centromere region and relatively lower usable coverage (13.7×). We hypothesized that additional coverage might offer us a better chance in accurately spanning the centromere. Therefore, we combined our Hi-C data with previously generated Hi-C and tethered chromosome confirmation capture (TCC) data. TCC is a Hi-C variant using solid support ligation³⁵ that generates similar data as a Hi-C experiment with slightly better ability to capture long-range chromatin interactions (Supplementary Fig. 7)³⁵. Using this combined data set, we increased the coverage of chromosome 9 to ~15×, which allowed accurate phasing of the entire chromosome. In summary, we generated complete chromosome- spanning haplotypes for all human chromosomes including chromosome X, albeit at reduced resolution of ~22% (Fig. 4b and Supplementary Fig. 6b), from just 17× genome sequencing using HaploSeq analysis.

Combining HaploSeq and local conditional phasing

Although we generated chromosome-spanning haplotypes using HaploSeq, we were unable to achieve a high resolution of variants phased owing to the low variant density in a human population. We reasoned that the gaps in the MVP block containing unphased variants (Fig. 4b) could be probabilistically linked to the MVP block using linkage disequilibrium patterns derived from population-scale sequencing data. For this purpose, we used the HaploSeq-generated, chromosome-spanning haplotype as a ‘seed haplotype’ to guide the local phasing using the Beagle (v4.0)³⁷ software and sequencing data from the 1000 Genomes Project¹⁵.

To initially assess the effectiveness of this approach, we simulated chromosome-spanning seed haplotypes in the GM12878 genome with different percentages of variants phased in the MVP block (Online Methods). Our simulation results suggest that we can accurately infer local phasing even at low-resolution seed haplotype inputs (3% error at 10% seed haplotype resolution; Fig. 5a). Owing to complex population structures, occasional mismatches occurred between phase predictions from local haplotypes predicted by Beagle and the HaploSeq-generated seed haplotype. To correct these mismatches, we filtered heterozygous variants with <100% agreement with the seed haplotype in a local neighborhood window surrounding the heterozygous variant. This filtering reduced the error rate to ~0.7% regardless of seed haplotype resolution (Fig. 5a). Consequently, the fraction of heterozygous variants for which we can infer local phasing increased with greater seed haplotype resolution (Fig. 5a). By contrast, altering the neighborhood window size did not substantially increase accuracy (Supplementary Fig. 8).

Encouraged by these results, we used the MVP chromosome-spanning haplotypes generated from HaploSeq analysis as seed haplotypes and performed local conditional phasing. Overall, we generated chromosome-spanning haplotypes with ~81% resolution at an average accuracy of ~98% (Table 2). Therefore, by coupling HaploSeq analysis and local conditional phasing, we achieved high-resolution and accurate chromosome-spanning haplotypes in humans.

Sequencing requirements for obtaining haplotypes by HaploSeq

From our local conditional phasing analysis, we deduced that a seed haplotype with ~20–30% resolution is sufficient to obtain accurate and high-resolution, chromosome-spanning haplotypes. A subsequent question therefore is, what are the minimal experimental requirements to achieve chromosome-spanning seed haplotypes with ~20–30% resolution? To investigate this, we simulated proximity-ligation sequencing data with varying read lengths and sequencing coverage. Based on our simulation, achieving chromosome-spanning haplotypes depends on obtaining a usable sequencing coverage of ~15× for most of the read lengths tested (Fig. 5b). However, chromosome-spanning seed haplotypes alone are not enough for achieving high-resolution haplotypes through local conditional phasing. In particular, the resulting sparse seed haplotype graph may limit the ability to generate final high-resolution haplotypes. To increase the resolution of the seed haplotype once complete seed haplotypes are obtained, one must increase coverage, either through higher sequencing depth or longer read lengths (Fig. 5c). We observed that 50- to 100-bp paired-end reads balanced completeness and resolution, and achieved the desired fraction of ~20–30% resolution at ~25–30× usable coverage. Although our simulations cannot address whether accuracy of predicted haplotypes depends on sequencing depth, we observed a relationship between haplotype accuracy and sequencing depth in the actual data. For instance, low coverage chromosomes such as 17 and 19 had a relatively lower accuracy. In particular, lower coverage might cause many variants to be linked with fewer edges, which in turn can propagate high-error structures to the entirety of chromosome-spanning haplotypes (Supplementary Fig. 9). Higher coverage can potentially increase accuracy even further, as observed in low-density CAST×J129 system. Therefore, 25–30× usable coverage with 100-bp paired-end reads is sufficient to achieve chromosome-spanning haplotypes with ~20–30% resolution and allow accurate local conditional phasing using HaploSeq analysis.

Cell culture and experimental methods

The F1 Mus musculus castaneus × S129/SvJae mouse embryonic stem cell line (F123 line) was a gift from the laboratory of E. Heard and has been described previously⁴⁴. These cells were grown in KnockOut Serum Replacement containing mouse ES cell media: DMEM 85%, 15% KnockOut Serum Replacement (Invitrogen), penicillin/streptomycin, 1× non-essential amino acids (Gibco), 1× GlutaMax, 1000 U/ml LIF (Millipore), 0.4 mM β-mercaptoethanol. F123 mouse ES cells were initially cultured on 0.1% gelatin-coated plates with mitomycin-C treated mouse embryonic fibroblasts (Millipore). Cells were passaged twice on 0.1% gelatin-coated feeder-free plates before harvesting. GM12878 cells (Coriell) were cultured in suspension in 85% RPMI media supplemental with 15% FBS and 1× penicillin/streptomycin.

Cells were harvested either in suspension (GM12878) or after trypsin treatment (F123 mouse ES cells). Formaldehyde fixation and Hi-C experiments were done as previously described³⁴. Of note, Hi-C data sets were used in our experimental method, but in principle, other proximity ligation methods such as Tethered Chromosome Conformation Capture (TCC) could also be used. We chose Hi-C as the primary proximity ligation method used in HaploSeq because it is a simpler method than TCC. As a result, this implies that the libraries produced by Hi-C are sufficient for successful haplotyping using HaploSeq despite the somewhat greater noise in the data, an important factor in the utility of the method for future applications.

Genotyping

Variant calls and genotypes for GM12878 were downloaded from ref. 43 and these were used for haplotype reconstruction. To clarify, “GM12878” is the name of the lymphoblastoid cell line whereas “NA12878” is the identifier for the individual from whom this cell line was derived. We have chosen to use GM12878 throughout this manuscript for the sake of consistency and clarity. Phasing Information for GM12878 was downloaded from 1000 Genomes Project¹⁴. The phasing of GM12878 by the 1000 Genomes Project covers only a fraction of heterozygous variants genotyped (by ref. 43) in this individual’s genome. Among these, only ~82% of the variants have phasing information owing to heterozygosity in both parents.

For generating genotype calls for the hybrid CAST×J129 cells, we downloaded parental genome sequencing data from publicly available databases. For Mus musculus castaneus, we downloaded the genome sequence from the European Nucleotide Archive (accession number ERP000042). S129/SvJae genome sequencing data was downloaded from the Sequence Read Archive (accession number SRX037820). Reads were aligned to the mm9 genome using Novoalign (http://www.novocraft.com) and using samtools⁴⁵, and we filtered out unmapped reads and PCR duplicates. The final aligned data sets were processed using the Genome Analysis Toolkit (GATK)⁴⁶. Specifically, we performed indel realignment and variant recalibration. The GATK Unified Genotyper was used to make single-nucleotide polymorphism (SNP) and indel calls. We filtered out variants that did not meet the GATK quality filters or that were called as heterozygous variants, as the genome sequencing was performed in homozygous parental inbred mice. The genotype calls in the parents were used both to determine the extent of interactions in cis versus h-trans to learn the phasing of hybrid CAST×J129 cells a priori to haplotype reconstruction.

Hi-C read alignment

For Hi-C read alignment, we aligned Hi-C reads to the mm9 (mouse) or the hg18 (human) genome. In each case, we masked any bases in the genome that were genotyped as SNPs in either Mus musculus castaneus or S129/SvJae (for mouse) or GM12878 (for humans). These bases were masked to “N” in order to reduce reference bias mapping artifacts. Hi-C reads were aligned iteratively as single-end reads using Novoalign. Specifically, for iterative alignment, we first aligned the entire sequencing read to either the mouse or human genome. Unmapped reads were then trimmed by 5 bp and realigned. This process was repeated until the read successfully aligned to the genome or until the trimmed read was less than 25 bp long. Iterative alignment is useful for Hi-C data because certain reads will span a proximity-ligation junction and fail to successfully align to the genome due to gaps and mismatches. Iteratively trimming unmapped reads has the potential to allow these reads to align successfully to the genome when the trimming removes the part of the read that spans the ligation junction. After iterative alignment of reads as single ends is complete, the reads are manually paired using in-house scripts. Unmapped and PCR duplicate reads are removed. The aligned data sets are then finally subjected to GATK indel realignment and variant recalibration.

Analysis of interaction frequencies between homologous chromosomes

When we aligned the Hi-C data, a paired-end read could either have both ends mapped to the same chromosome (intra-chromosomal) or mapped to different chromosomes (inter-chromosomal). However, the initial mapping of the Hi-C data utilized a haploid reference genome and did not distinguish to which of the two homologous copies of a chromosome an individual sequencing read maps. As a result, read pairs that initially map as intrachromosomal should be broken down into reads that occur on the same homologous chromosome (which are truly in cis) and reads that map between the two homologous pairs (which we define in this manuscript as h-trans).

To determine the extent of reads that are in cis versus h-trans, we first distinguished to which allele an individual read mapped. This is done simply by identifying reads that overlap with variant locations in the genome and then determining which allele the sequenced base at the variant location corresponds to. Once this information is obtained, we can determine the frequency with which regions interact in cis versus h-trans (as shown in Fig. 2b). For Figure 2c and Supplementary Figure 2, we checked interacting variants (through reads) at any given insert size.

Usable coverage as defined by intra- and interchromosomal reads

For phasing using HapCUT, we utilize both intrachromosomal and interchromosomal reads. For interchromosomal reads, we consider each interchromosomal read pair as two single-end reads, as the paired information for such reads is not useful for phasing. In contrast, all intrachromosomal reads are considered for phasing. The probability of a single read to harbor more than one variant is small, especially in humans where the variant density is relatively low. This, in combination with the fact that only the paired intrachromosomal reads will have large insert sizes, means that the vast majority of reads that contribute to the success of haplotype phasing are the intra-chromosomal reads. Therefore, we define the “usable coverage” as the genomic coverage derived from intrachromosomal reads only.

Our Hi-C experiment generated ~22% interchromosomal reads in CAST×J129, whereas ~55% of the reads in GM12878 were interchromosomal. In other words, 620 M paired-end reads out of 795 M were useful in CAST×J129, with a usable coverage of 30×. In humans, only 262 M paired-end reads out of 577 M were useful, resulting in a usable coverage of 17×. Thus, we have a lower usable coverage in humans. In our experience, the fraction of all reads that are intrachromosomal versus interchromosomal in a Hi-C experiment may vary between experiments and across cell types.

Analysis of HaploSeq data using HapCUT

We used the HapCUT²⁴ algorithm to perform the computational aspects of HaploSeq. This method was originally designed to work on conventional genome sequencing (WGS) or mate-pair sequencing data. HapCUT constructs a graph with the heterozygous variants as nodes and DNA fragment(s) connecting two nodes as edges. Therefore, only fragments with at least two heterozygous variants are useful for haplotype phasing. HapCUT extracts such ‘haplotype-informative’ fragments from a coordinate-sorted BAM file using a sorting method that stores each potential haplotype-informative read in a buffer until its mate is seen. We increased the buffer size to allow HapCUT handle large insert–sized proximity-ligation reads.

HapCUT uses a greedy max-cut heuristic to identify the haplotype solution for each connected component in the graph with the lowest score under the MEC scoring function. In particular, the original HapCUT algorithm used O(n) iterations to find the best cut for a component with n variants. Because Hi-C data resulted in chromosomal spanning haplotypes with a single large connected component, the default method took several days of computing time to phase the CAST×J129 genome. To reduce the computation time, we assessed the impact of reducing the number of max-cut iterations on the accuracy of phasing. For CAST×J129 system, we observed that increasing the number of max-cut iterations beyond 1,000 did not significantly improve the accuracy. For GM12878, we allowed up to 100,000 iterations. In summary, the higher the number of heterozygous variants in the largest connected component of the graph, the lower this parameter should be to make predictions within a reasonable time.

HapCUT uses the best cut for each component to update the current haplotype solution and the edge weights in the graph. This process is iterated multiple times until a combinatorial object function can no longer be improved. We used a maximum of 21 such iterations in CAST×J129 and 101 in GM12878 cells. Our parameters in GM12878 cells allowed HapCUT to obtain higher accuracy given the lower variant density and reduced sequence coverage compared to the mouse data. The modified version of HapCUT can be downloaded from https://sites.google.com/site/vibansal/software/hapcut. The source code is available as Supplementary Data.

Maximum insert size analysis

As previously mentioned the probability of a Hi-C read being in cis versus h-trans varies as a function of the distance between the two read pairs (Fig. 2c). At shorter genomic distances, the probability that an intrachromosomal read is in h-trans is very low. At large distances (>30 Mbp), this probability rises substantially and is in theory more likely to introduce erroneous connections for HapCUT to phase. To account for this, we used the Hi-C data for chromosomes 1, 5, 10, 15 and 19 in the CAST×J129 data and repeated haplotype reconstruction allowing variable maximum insert size values. We simply excluded any reads where the insert size between reads was greater than the allowable maximum insert size. We performed this analysis using the low variant density case, for this analysis because lower density was most amenable for applications in humans (Supplementary Fig. 4a,b). This step resulted in increase in accuracy of HaploSeq analysis with moderate reduction in resolution.

Insert size–dependent probability correction

A useful feature of the HapCUT algorithm is that it accounts for the base quality score at a variant location to calculate the score of a potential haplotype. In other words, if a sequencing read that links two variants and the base quality at one variant location is low, this read is given relatively lower weight by HapCUT in generating its final haplotype calls. Therefore, HapCUT can use this information to try to disregard potential sequencing errors from making erroneous haplotype connections. As we previously mentioned, in Hi-C data errors may also arise due to h-trans interactions, which are much more frequent than sequencing errors and show a distance-dependent behavior. Therefore, we attempted to account for the likelihood of an interaction being in cis versus h-trans based on the distance between the two reads. We used the CAST×J129 Hi-C data to identify reads that are in cis or h-trans. We binned the insert sizes into 50-kb bins and estimated the probability of a read being h-trans (#h-trans/ (#cis+#h-trans). We then used local regression (LOWESS) at 2% smoothing to predict h-trans probabilities at any given insert size. For every intrachromosomal read, we multiplied the cis probabilities (1 – h-trans) with the base qualities to account for the odds of this intrachromosomal read being a h-trans interaction. As a result, reads that are more likely to be h-trans are given lower weight by HapCUT in identifying the haplotype solution.

Adding h-trans interaction probabilities increases HaploSeq accuracy moderately, without having any affect on resolution. As a comparison, maximum insert size of 30 Mb had an error rate of 1.1% in chromosome 19 (Supplementary Fig. 4b). After adding h-trans probabilities, the error rate is 0.9% (Supplementary Fig. 6a), where error rate is defined as 1 – accuracy.

Local conditional phasing simulation

In order to study our ability to perform local phasing at different percentages of resolution, we performed a stepwise analysis. First, we generated seed haplotypes at different resolutions. Then, we used Beagle (v4.0)³⁷ to perform local phasing under the guidance of the seed haplotype. Finally, we checked accuracy of local phasing by comparing it to phasing information known a priori from 1000 Genomes Project.

To simulate seed haplotypes at different resolutions, we first simulated seed genotypes. We used different combinations of read length and coverage to obtain seed genotypes of various resolutions. In particular, we used Hi-C intra-chromosomal read starting positions from H1 and H1-derived cells (unpublished data) to generate pairs of reads of a given read length and coverage. This allowed us to maintain the Hi-C data structure and the observed distribution of insert sizes in the simulated data. To generate the seed genotype, we constructed a graph with nodes representing heterozygous variants in GM12878 (chromosome 1) and edges corresponding to reads that cover multiple variants. This graph is essentially a genotype graph because we don’t know the phasing yet. Hence, the whole point of this graph is to provide a two subset of variants: one that is a part of the seed genotype and other that is not (which are the gaps to be inferred by local phasing), based on the resolution and Hi-C data structure. We generated seed genotypes at our required parameters of read length and coverage to attain a specific resolution. We used these seed genotypes for both local phasing (Fig. 5a) and to study the minimal requirements for generating seed haplotypes of enough resolution (Fig. 5b,c). These two analyses were done independently and in both cases, we repeated generating seed genotypes and downstream analysis ten times to note the average results.

To perform local conditional phasing, we need an a priori haplotype system to check accuracy of our local conditional phasing. Because a priori haplotype information from the trio covers only a fraction of heterozygous variants, we decided to perform local phasing simulation only on the trio subset. Specifically, we required every variant that was part of either seed genotype or “gaps” to be part of the 1000 Genomes-phased trio. We converted seed genotypes to seed haplotypes using the trio information while keeping “gap” variants as unphased. We then used local phasing conditioned on the seed haplotype to infer phasing of the gap variants using Beagle. We allowed homozygous variants to assist Beagle in making better predictions from the Hidden Markov Model.

To perform neighborhood correction for a seed haplotype unphased variant, we collected three variants each from both upstream and downstream, that are phased in seed haplotype. Then we checked if there was 100% correlation between the phasing present in the seed haplotype to what is predicted by Beagle. This provides an estimate of how well Beagle could have performed in this “local” region. If there is a 100% match, we consider the variant as conditionally phased. If there is not a 100% match, we disregarded the unphased variant in the final haplotype. We tried other window sizes such as 5 and 10 (Supplementary Fig. 8) and found no improvement in accuracy.

Local conditional phasing in human GM12878 cells

We coupled HaploSeq analysis and local conditional phasing to increase resolution in GM12878 cells. Local conditional phasing was performed as described earlier on genotypes that are common between GM12878 (ref. 43) and population samples. In addition, as the seed haplotype is not 100% accurate, we marked the seed haplotype phased variants that did not agree with local phasing. These marked variants were made “unphased” as these could be potential errors from HaploSeq. Hence, apart from using neighborhood correction for deciding whether a gap variant needs to be locally phased (as in the simulation), we also used this information to mark variants in the seed haplotype that could be potentially erroneous. This allowed a minor increase in accuracy after local phasing (Supplementary Fig. 9).

Overall HaploSeq accuracy was estimated as the fraction of heterozygous variants correctly phased in the MVP block after local phasing (Table 2 and Supplementary Fig. 9). In particular, we used only the variants phased in trio to estimate accuracy. For local phasing in chrX, we made the male haploid genotypes homozygous.

GM12878 cells have a lower variant density than CAST×J129 and a lower coverage added more constraints on the prediction model resulting in a relatively higher HaploSeq error rate of 2%, when compared to 0.8% in low density CAST×J129 case. A usable coverage of 25–30× (as shown in Fig. 5b,c) could help gain accuracy and potentially cover more rare variants in the seed haplotype. Currently, about 16% of the variants are not locally phased due to their absence in population. In particular, we used rsIDs to match variants, as our predictions were in hg18 whereas the 1000 genomes data were in hg19 format. These could be phased either by additional Hi-C data or even conventional genome sequencing data, which can potentially link gap variants to variants in seed MVP block. The important aspect in HaploSeq analysis is the ability to form seed chromosome-spanning haplotype, which cannot be made from conventional genome sequencing or mate-pair or fosmids.

Fosmid simulations

To simulate fosmid-based sequencing (Fig. 3b,c and Supplementary Fig. 1b), we emulated fosmid clones as paired-end sequencing, with insert sizes close to 40 kb. We reasoned that this approach is easier to simulate and yet maintains the data structure that fosmids add to haplotype reconstruction. For example, our simulation produces haplotype blocks of up to 1 Mb in humans (Supplementary Fig. 1b), a size that has been reported by other groups^25–27.