article: Humans are a diploid species that inherit one set of chromosomes paternally and one homologous set of chromosomes maternally. Unfortunately, most human sequencing initiatives ignore this fact in that they do not directly delineate the nucleotide content of the maternal and paternal copies of the 23 chromosomes individuals possess (i. e. , they do not ‘phase’ the genome) often because of the costs and complexities of doing so. We compared 11 different widely-used approaches to phasing human genomes using the publicly available ‘Genome-In-A-Bottle’ (GIAB) phased version of the NA12878 genome as a gold standard. The phasing strategies we compared included laboratory-based assays that prepare DNA in unique ways to facilitate phasing as well as purely computational approaches that seek to reconstruct phase information from general sequencing reads and constructs or population-level haplotype frequency information obtained through a reference panel of haplotypes. To assess the performance of the 11 approaches, we used metrics that included, among others, switch error rates, haplotype block lengths, the proportion of fully phase-resolved genes, phasing accuracy and yield between pairs of SNVs. Our comparisons suggest that a hybrid or combined approach that leverages: 1. population-based phasing using the SHAPEIT software suite, 2. either genome-wide sequencing read data or parental genotypes, and 3. a large reference panel of variant and haplotype frequencies, provides a fast and efficient way to produce highly accurate phase-resolved individual human genomes. We found that for population-based approaches, phasing performance is enhanced with the addition of genome-wide read data; e. g. , whole genome shotgun and/or RNA sequencing reads. Further, we found that the inclusion of parental genotype data within a population-based phasing strategy can provide as much as a ten-fold reduction in phasing errors. We also considered a majority voting scheme for the construction of a consensus haplotype combining multiple predictions for enhanced performance and site coverage. Finally, we also identified DNA sequence signatures associated with the genomic regions harboring phasing switch errors, which included regions of low polymorphism or SNV density.
  Whole genome sequencing (WGS) of individual human genomes has been made possible by rapid advances in DNA sequencing technology (reviewed in [1]). In fact, an initial draft of an individual’s genome, which results in spelling out over 3. 2 billion base pairs, now costs approximately $1,000. The exploitation of WGS in association and clinical studies has led to some remarkable successes in the identification of disease-causing, overtly pathogenic variants [2] as well as the use of these variants in diagnostic, prognostic and pharmacogenetic assessment settings [3]. Unfortunately, due to inherent limitations in current sequencing technologies, such as short DNA sequence read lengths, a number of issues remain. For example, sequencing strategies are not error-free, making it possible that some nucleotides in an individual genome are erroneously assigned. This error varies greatly among sequencing technologies [1,4]. In addition, although virtually all sequencing technologies can capture simple single nucleotide variants (SNVs), the identification of some forms of structural variation (e. g. , moderately sized insertion/deletion variants (indels) and complex copy number variants (CNVs) ) have proven difficult to identify in a straightforward manner [5,6]. Steady improvements in sequencing laboratory protocols and computational workflows are beginning to overcome these issues [6], but they also contribute to another, perhaps more fundamental issue surrounding the fact that humans are diploid and possess two ‘genomes: ’ one inherited maternally and one paternally as a set of 23 chromosome pairs. Understanding how cis (i. e. , near each other on the same chromosome) and trans (far from each other or on opposite chromosomes) combinations of variants impact phenotypic expression is likely to be of crucial importance in human biology. This activity depends critically on having phase information, i. e. , knowing the unique nucleotide content of each of the two chromosomes making up the 23 pairs. Most genome sequencing strategies derive DNA sequence for the 23 chromosome pairs simultaneously from a pool of maternally and paternally-derived DNA (e. g. , from blood samples) and hence do not segregate the nucleotide content of each copy of a homologous chromosome pair. Further, most strategies seek to identify nucleotides at each position of a genome and, by comparing nucleotides of the sequenced genome to the sequence of nucleotides for a vetted reference genome, identify SNVs and possibly other forms of variation. In general, each SNV position is then recorded as having either one copy of a non-reference nucleotide and one copy of a reference nucleotide (i. e. , ‘heterozygous’), two copies of a non-reference nucleotide (i. e. , ‘homozygous non-reference’), two copies of the ‘reference’ nucleotide (i. e. , ‘homozygous reference’), or, rarely, two different non-reference nucleotides. Note that phase information is really only relevant for heterozygous loci (i. e. , where there is a difference in the nucleotide content of homologous chromosome pair). Obtaining phase (or haplotype) information is often non-trivial, as reviewed by us and other research groups [7–9]. As noted, there are a number of very important genetically-mediated phenomena for which characterization depends critically on phase information. For example, the phenomenon of compound heterozygosity, in which two different mutations, one impacting the function of paternally-inherited homologous copy of a gene and the other impacting the function of the maternally-inherited homologous copy of that gene, can only be distinguished from situations in which the two mutations reside on the same copy of the gene by having phase information [7,10]. This is particularly important for the highly polymorphic human leukocyte antigen (HLA) region of the genome. In addition, phase information is also required in dissecting many evolutionary and population genetics phenomena, such as haplotype frequency differences across populations and linkage disequilibrium (LD) patterns. For example, LD patterns around specific mutations are useful for determining the age of those mutations as well as imputing variants to individuals that have not been sequenced or genotyped at the relevant position in the sequence [11]. A number of strategies have been, or are being, developed to generate phase information [7,8]. Some of these strategies involve elaborate laboratory protocols, others rely on de novo assembly computational methods, and yet others make use of population haplotype frequency information. Sequencing read lengths and the depth (i. e. , coverage) of sequencing for a genome play large roles in the ability to phase a genome. By far the most straightforward, most often-used, and most-intuitive strategy for obtaining phase information from genome sequencing involves sequencing the parents (and/or other relatives) of an individual whose genome is to be phased. This way, Mendel’s laws can be used to essentially observe which variants were inherited from each of the parents. Unfortunately, even with parental sequence information there are problems that can plague obtaining phase information. For example, if the mother, father, and offspring are all heterozygous at a locus, the parental origin of the variants in the offspring cannot be determined, and the phase information for variants at that locus will remain ambiguous. We compared a variety of methods for phasing entire human genomes leveraging the publicly available genome of a European (i. e. , ‘CEU, ’ per community conventions) female individual NA12878. The NA12878 ‘Genome-In-A-Bottle (GIAB) ’ reference was constructed based on multiple-platform arbitrated genotype [12] and haplotype information by multiple research groups at Illumina, Inc. and Real Time Genomics, Inc. using a suite of techniques and sequence information available from the relatives of NA12878 [13,14]. To actually compare the performance of the different phasing strategies, we considered laboratory-based phasing methods by analyzing previously published datasets (i. e. , we used what was published, since we did not have access to the raw material and laboratory protocols to recreate them), as well as computational phasing strategies that we were able to put together based on available software. We recognize that the use of publicly available genomes could be biased since the research groups generating the data may have used many protocols and workflows that are not necessarily available to the public, but we argue that our analyses of those public domain phased genomes provides an ‘upper bound’ (i. e. , best case scenario) on their accuracy. We benchmarked the accuracy of the laboratory and computational phasing approaches using a variety of metrics including switch error rates, haplotype block lengths, and the proportion of fully phase-resolved genes. Our results suggest that hybrid computational approaches to phasing a single genome are accurate and cost-effective, and exhibit a performance that nearly matches or surpasses laboratory-based approaches.
  We compared 11 different approaches for phasing the NA12878 genome. The approaches considered are briefly described below with references provided where appropriate. We also refer the reader to websites associated with the commercial groups whose products enable some of these approaches. The descriptions below consider the laboratory-based methods first, followed by computational methods. Population-based phasing, such as those strategies implemented in Beagle, Eagle2, and SHAPEIT software packages, makes use of a haplotype reference panel for haplotype frequency information. For this study, the 1000GP phase 3 reference panel dataset was used and contains 2,503 individuals excluding individual NA12878. We also used the Haplotype Reference Consortium (HRC) panel dataset which contains 22,690 individuals excluding individual NA12878 [34]. The HRC reference panel included WGS datasets from over 20 studies to develop a large combined haplotype reference panel, currently with predominantly European ancestry. The HRC reference panel also included the 1000GP individuals. The 1000GP reference panel for Beagle and SHAPEIT runs were obtained from the corresponding project web sites of the phasing tools [35,36]. The HRC reference panel was obtained from the European Genome-Phenome Archive (EGA) (Dataset ID: EGAD00001002729) [37]. We considered five different metrics to assess phasing performance. All have been used in previous research to assess how well individual genomes can be phased. Each of these metrics is described briefly below. The datasets and scripts used for comparing the phasing methods is available at https: //github. com/ywchoi/phasing.
  A total of 11 phasing approaches were compared based on the NIST-GIAB reference genome NA12878. Laboratory-based experimental phasing previously reported from 10X Genomics, CPT-seq, fosmid sequencing, and Moleculo were collected from the literature and public datasets. Computational phasing results were generated in-house for the NA12878 genome using Beagle [32], Eagle2 [30], SHAPEIT [31], and HapCUT [33] with one or more combinations of the following as input data for the relevant computations: the haplotype reference panels from the 1000 Genomes Project (1000GP) with 2. 5k individuals or the haplotype references from the Haplotype Reference Consortium (HRC) with 23k individuals, conventional WGS reads generated from the Illumina or PacBio sequencing platforms, RNA-Seq reads, or parental genotype. To determine what, if any, impact the composition of a reference panel had on the strategies relying on population-based approaches, we examined the effects of the ancestral origins and size of the reference panel using SHAPEIT. The SERs were measured using the same number of individuals (n = 347) from each of well-known population ‘supergroups’ associated with the 1000GP data, namely AFR, AMR, EAS, EUR, and SAS, for the assessments of the SHAPEIT-based strategies. The EUR supergroup, when used as the reference panel alone, produced the lowest switch error rate (Fig 6A). We further investigated the effect of ancestral subpopulations on phasing accuracy, by dividing the EUR supergroup into well-known 1000GP subpopulations: CEU, FIN, GBR, IBS, and TSI, with each group including the same number of individuals (n = 85). We compared the phasing accuracy of the subpopulations against the supergroup EUR and the overall 1000GP mixed population at the same group size (Fig 6B). The SHAPEIT switch error rates obtained for CEU, EUR, and 1000GP were 1. 9,2. 0, and 2. 8%, respectively. Given that the reference genome NA12878 belongs to the EUR population and the CEU subgroup, our results suggest that a reference panel that contains matching ancestry with the query genome at the supergroup (i. e. EUR) level, and not necessarily at the subpopulation level (i. e. CEU), was sufficient to achieve excellent phasing performance. We also examined the effect of reference panel size by randomly subsampling from each of the 1000GP EUR (n = 502) and HRC (n = 22,690) panels. For both panels, the phasing accuracy was improved by increasing the reference panel size, and the improvement diminished as the panel size increases in all SHAPEIT approaches with and without WGS reads or parents’ genotype data. Specifically, when 50X PacBio reads were supplied, only ~500 HRC individuals were required to obtain the same accuracy as achievable by using the entire 23k HRC reference panel alone (i. e. ~0. 3% SER). Furthermore, when familial information from parental genotype data was used, the SER obtained was 0. 1% with 100 HRC individuals and as low as 0. 03% when using the entire panel of 23k HRC individuals (Fig 6C).
  Current human WGS strategies do not provide phase information by default. This creates issues for identifying certain clinically-meaningful phenomena, such as compound heterozygosity, as well as population-level phenomena such as haplotype diversity and LD patterns that could help resolve migratory patterns and mutation origins. We evaluated 11 different available strategies to resolve phase for human genomes. In general, these strategies fell into 3 different categories: those that rely purely on building haplotypic contigs from DNA sequencing reads or unique laboratory-based approaches through genome dilution, such as the 10X, CPT, fosmid pools, and Moleculo; those that rely on purely population-based strategies in which a reference panel of previously resolved haplotypes is leveraged to probabilistically assign the most likely haplotype pair to an individual based on their genotype data (i. e. Beagle, Eagle2, SHAPEIT); and those that rely on some combination of sequence read- and population-based constructs, and/or parent offspring trios genotype information. Our results suggest that: (1) the recently developed large comprehensive reference panel by the HRC greatly enhance phasing performance; and (2) population-based phasing performance can be greatly enhanced with family trio data that can match or surpass the precision of laboratory-based methods. In all, laboratory-based methods are more expensive (e. g. , $2,000 per genome using the 10X Genomics approach [40]), and often require protocols that add complexity to workflows needed to use them but are likely to come down in costs. For cost-effectiveness, a combination of computational methods and input data types as shown in this study likely worth considerations. In our view, by far the most efficient and least expensive current approaches to phasing are likely going to be based on computational methods leveraging large reference panels such as the HRC. Our analyses showed that population-based phasing performance can be improved if WGS/RNA read data are available. Long read-based strategies (e. g. PacBio) provide an additional advantage over other sequencing strategies for their ability to link kilobase-range long distance variants. However, the best improvement to population phasing performance, which involved a 10-fold reduction of errors and also a 10-fold increase in haplotype block lengths, was the addition of parental genotype data. We also proposed a majority voting approach to exploit the possibility of combining phasing strategies to improve the overall SNV site coverage and the number of phase-resolved genes, and also to build consensus from multiple predicted haplotypes to simplify downstream analysis. Traditional strategies to phasing genomes that involve sequencing the parents (or relatives) or a target individual are enjoying a resurgence because of that strategy’s ability to also help validate de novo or private mutations in clinical contexts. Although such a strategy appears costly in terms of the amount of sequencing that needs to be performed as it also depends on the availability of relatives, the continuous reduction in sequencing cost could favor the familial approach for improving the performance of computational phasing methods. With regard to the fosmid pool sequencing dataset, for comparison with the NA12878 GIAB gold standard and other phasing approaches we used the liftOver tool [26] to convert the genomic coordinates from the NA12878 fosmid study from NCBI36 (hg18) to GRCh37 (hg19). As the liftover process for NA12878 fosmid data may have led to loss of some haplotype blocks, the fosmid-based phasing result for the European individual ‘‘Max Planck One” (MP1) from a previous study is mentioned here for additional comparison. For the NA12878 individual, approximately 0. 48 million fosmid clones were sequenced (i. e. 32 pools with 15,000 clones per pool), which was estimated to provide 5X coverage of the genome with 77% of SNVs phased [10]. For the MP1 individual, deep sequencing of approximately 1. 44 million fosmid clones (96 pools) equivalent to 14X coverage of the genome had led to 99% SNVs and 81% genes phased [10]. We focused on phasing performance for the autosomal chromosomes. Because chromosome X has a number of specific features, such as variants within it exhibiting unique inheritance patterns, it does contain genomic features, including the pseudo-autosomal regions (PAR) and non-PAR segments, that would require special data handling and analysis modifications to be accommodated in our analyses. For NA12878, who is a female, we attempted a phasing comparison for the non-PAR region of chromosome X. The SER obtained with 10X and Moleculo was 0. 15 and 0. 60% respectively, which was about 2. 43 and 1. 19 times higher than the autosomes. For chromosome X, SHAPEIT phasing produced haplotype blocks with SER at 0. 3%. The addition of parental genotype reduced the SER to 0. 05% (S4 Table). Our analyses do raise some important questions. First, the accuracy of our results relies on the NIST-GIAB gold standard variants used in this study. Recently, a further improved high confidence NA12878 variant dataset (Platinum genotype calls) was developed and validated using extensive family information spanning 3 generations and 17 members, and incorporating 6 different variant calling pipelines [41]. For the NA12878 individual, a total of 3. 5 million SNVs (heterozygous or homozygous) were identified in the Platinum variant dataset containing 800,564 additional SNVs not included in the NIST-GIAB set, and conversely a set of 62,946 SNVs in the NIST-GIAB set are not included in the Platinum dataset [41]. Nevertheless, over 97. 7% SNVs of the NIST-GIAB variant set used in this study overlapped with this high confidence variant dataset and have high genotype concordance level (>99. 99%) [41]. Using the Platinum dataset [41] as a gold standard, we observed similar phasing performance for the different approaches tested (S3 Fig). Second, we obtained phased genomes deposited in the public domain by, e. g. 10X and Moleculo, and these were likely put together by the relevant research teams to showcase these technologies. Thus, the results of these strategies may have been pursued with some level of optimization that is not likely to be obtained with routine and large-scale use of these strategies. Finally, our analyses focused on SNVs and did not consider short insertion-deletion polymorphisms (indels) and structural variants such as moderate to large indels, copy number variants (CNVs), inversions, or other forms of structural variations (SVs). Such variations have a significant role to play in many diseases and are hence very important to consider. In addition, certain combinations of SVs and CNVs might be very important to resolve for population-based studies. For example, variants within different copies of a gene may help resolve the origins and recent evolution of those copies. Ultimately, since humans are indeed a diploid species, the ability to phase human genomes using fast and low-cost computational approaches will push diploid genome analysis to a new level.
summary: Humans are a diploid species that inherit one set of chromosomes paternally and one set of chromosomes maternally. Separating the nucleotide content of the maternally and paternally-derived chromosomes for an individual, i. e. , ‘phasing’ that individual’s genome, is not trivial with today’s sequencing technologies. This is in part due to the fact that most available sequencing technologies generate short sequencing reads that make it hard to assemble individual homologous chromosome pairs. Phase information can be crucial for putting into context the likely functional consequences of DNA sequence variants as well as certain evolutionary and population genetics phenomena. In order to assess the reliability of current sequencing-based phasing strategies, we compared 11 different approaches using a public domain reference genome as a test case. These phasing strategies included laboratory-based experimental techniques as well as purely computational approaches. Importantly, our comparisons show that a hybrid or combined approach that leverages population-based phasing via the SHAPEIT software suite works well and can be improved with the addition of genome-wide sequence read or parental genotype data.
keywords: sequencing techniques
  variant genotypes
  human genomics
  sequence assembly tools
  genetic mapping
  genome sequencing
  genome analysis
  molecular biology techniques
  dna
  research and analysis methods
  comparative genomics
  molecular biology
  nucleotide sequencing
  biochemistry
  haplotypes
  nucleic acids
  heredity
  genetics
  biology and life sciences
  genomics
  computational biology
title: Comparison of phasing strategies for whole human genomes