WIXLIB

Technical reports 2014 Nature America, Inc. All rights reserved.Figure 2 Overview of the CPT-seq workflow.There are three key steps: (I) indexedtransposition, (II) pooling, dilution andcompartmentalization, and (III) indexed PCR.A set of 96 different indexed transposomecomplexes are used to set up 96 independenttransposition reactions to create separatevirtual genomic partitions (step I). Transpositionreactions are pooled together, diluted tosubhaploid DNA content and split into 96compartments (step II). Upon removal of thetransposase with SDS, compartment-specificlibraries are generated using indexed PCR(step III). All samples are pooled togetherafter PCR and prepared for sequencing.(I) 35 min Indexedtransposome Indexedtransposome EDTA(II) 10 minPoolingphysical compartment, using a combinatorialtwo-step indexing scheme (Fig. 2). The firstindex is incorporated into the genomic DNAlibrary during transposition (defining the virDilution and compartmentalizationtual partitions), and the second index is incorporated during indexed PCR of each physicalcompartment. In this manner, a set of virtualpartitions, equal in number to the indexedtransposition reactions in the first step, can(III) 120 min SDSbe defined within each physical compartmentof the subsequent dilution step.Transposon Genomic TransposonAs a concrete example, a single genomicindex 1DNAindex 2sample is split into M 96 independentIndexed PCRPCRPCRtransposition reactions, each employingindex 1index 2uniquely indexed transposon adaptors.These 96 separately indexed and adaptorized genomic libraries (contiguitypreserved) are then pooled, diluted andredistributed into a set of N 96 discrete physical compartments method. After physical compartmentalization, the DNA-transposome(wells on a plate), with each physical compartment now having complexes are denatured by the addition of SDS detergent, and theM 96 virtual partitions. A second compartmental index (N 96) is content of each compartment is amplified with a pair of indexed PCRincorporated during PCR. In this manner, this 2-tier indexing process primers (96 different PCR indices).effectively creates a total of M N 96 96 9,216 virtual compartAfter PCR, combinatorially indexed libraries from all wells are comments whose library elements are demultiplexed after sequencing by bined and sequenced using a dual-indexing workflow to read all fourreading out the unique set of combinatorial indices. In addition to index sequences (two bipartite codes, each consisting of left and rightthe two tiers of indexing, index sequences on a single tier are created indices) and genomic DNA inserts (Supplementary Fig. 3). A standardfrom the combination of ‘left’ and ‘right’ index sequence combina- IVC plot (sequencing intensity versus cycle number; Supplementarytions (Supplementary Fig. 3)24. A set of only 8 12 oligonucleotides Fig. 4) confirmed the effectiveness of the four-primer sequencis used to create the 96 asymmetrically indexed transposomes, and ing strategy, with sequence reads through genomic DNA (read 1),an additional 8 12 oligonucleotide primers are used to create the 96 both bipartite indices (reads 2 and 3) and genomic DNA (read 4).asymmetrically compartmental PCR indices (Online Methods).In summary, genome-wide haplotype information can be capturedOur implementation of this concept consists of taking a genomic with a simple workflow consisting of three steps: (i) parallelized andDNA sample and aliquotting 1 ng into 96 different wells of a micro- indexed transposition of genomic DNA, (ii) pooling, dilution andtiter plate, with each well containing a Tn5 transposome mix with a physical compartmentalization of libraries after transposition andunique index (M 96; virtual partitions). After tagmentation, the 96 (iii) parallelized and indexed PCR.indexed genomic libraries (contiguity preserved) are pooled, dilutedand redistributed into 96 separate wells of a microtiter plate (N 96; CPT-seq haplotyping resultsphysical compartments). The original pool is diluted such that there As proof of concept, we applied this strategy to the haplotype-resolvedis approximately 3% haploid content per virtual partition or 3 copies whole-genome sequencing of a classical HapMap family trio18 (NA12878,of the genome per physical compartment (96 virtual partitions per NA12891 and NA12892). We prepared high-molecular-weight inputphysical compartment). Such low haploid content per virtual com- DNA for the trio samples using the Gentra PureGene kit (Qiagen) withpartment is required to avoid collisions between maternal and pater- an expected average size of 100–200 kb (Supplementary Fig. 5). For thenal copies from the same genomic region. Notably, given the virtual NA12878 sample, we also used a commercially available genomic DNAportioning, the amount of DNA per physical compartment is two source with low molecular weight (i.e., an average DNA size of 50 kb).orders of magnitude higher than in other dilution-based haplotyping All DNA samples were processed through the described workflow andapproaches, contributing to the robustness of amplification in our sequenced on four lanes of a HiSeq 2000, generating 80–130 Gb of dataNature GeneticsADVANCE ONLINE PUBLICATION

Technical reportsTable 1 Haplotyping results for the trio samplesNA12878NA12878DNA prepCoriellGentraRead length7651Number of reads (million)1,7491,622Number of HiSeq lanes44Mapped bases (Gb)13383DNA/partition (Mb)2158Informative island4576N50 (kb)SNPs covered (%)93.1598.46Phasing block N50 (kb)4901,485Long switches/Mb0.140.18Step I percent phased90.6496.81Step I accuracyLA: 99.97 LA: 99.96PA: 99.41 PA: 99.58Step II percent phased71.7192.84Step II accuracyLA: 99.96 LA: 99.96PA: 99.82 PA: 99.89Step III percent phased94.698.17Step III accuracyLA: 99.97 LA: 99.96PA: 99.83 PA: 99.88Intra-island :PA:LA:PA:LA:PA:Whole-genome phasing data are provided for the samples NA12891, NA12892 andNA12878 prepared using the Gentra PureGene Cell kit or acquired from Coriell(NA12878 only; Online Methods). The average DNA content of all partitions is shownin megabases and reflects the actual fraction of material added per virtual compartmentafter accounting for losses. In the ideal case, 3% haploid content ( 90 Mb) was addedper virtual partition. The analysis is carried out in three main steps: step I (ReFHapphasing), step II (removing single-linked and conflicting SNPs) and step III (addingmissing SNPs using 1000 Genomes Project data). The fraction of SNPs phased aswell as point-switch (PA) and long-switch (LA) accuracies are reported for each step ofthe analysis. N50 values are reported for both haplotyping islands and the haplotypingblocks obtained after phasing individual islands into longer contigs.mapping to the genome (Table 1). Sequence reads were demultiplexedinto 9,216 distinct partitions on the basis of their unique combinatorialindices. Data from all partitions were mapped to the reference humangenome (hg19). Schematic and representative genome coverage plotsfor the mapped reads from three representative indices are shown inFigure 3 and Supplementary Figure 6. Clusters or islands of reads wereobserved scattered across the genome, suggesting that these islands originated from a single molecule. Furthermore, a plot of nearest-neighbordistances between reads for a given index across a genomic regionexhibited a bimodal distribution (Fig. 3, bottom, and SupplementaryFig. 7). The reads from within an island (proximal reads) formed onepeak, and the reads between islands or sparse singleton reads (distalreads) formed the other peak. The coverage within each island (seeFig. 3, top, and Supplementary Fig. 8 for a schematic and representative plot, respectively) exhibited a ‘strobed’ pattern, with only 5–10%coverage. Despite this low intra-island coverage per indexed partition,the combined genome-wide coverage from all 9,216 partitions wasbetween 97 and 99% (Table 1). Using high-molecular-weight DNA, theN50 of informative islands was between 70 and 90 kb, and switching tolow-molecular-weight DNA decreased the island size to 45 kb. About50–60 Mb of genome, or 2% haploid equivalents, was observed pervirtual compartment (Table 1).Assignment of SNPs to their respective haplotypes was performed ina multi-step process (Online Methods and Supplementary Fig. 9)25.ReFHap26 phased 94–97% of the SNPs from the high-molecularweight DNA samples (Table 1). The N50 of the assembled haplotyping blocks for high-molecular-weight DNA samples ranged from 1.4to 2.3 Mb (Supplementary Fig. 5). We evaluated the yield and theaccuracy of the current phasing method by plotting the probabilitythat heterozygous SNP pairs were on the same phasing block as a Virtual compartment 1Virtual compartment XVirtual compartment 0Count 2014 Nature America, Inc. All rights reserved.DNA source600400Distal reads(inter-island)20000246Distance (bp) to next alignment (log10)Figure 3 Demonstration of haplotype read islands. Coverage plots areshown for three representative indices across part of chromosome 22.Reads from the same contiguous molecule display as islands of readclusters across the genome (middle) or as one mode of a bimodaldistribution from a nearest-neighbor plot of mapped reads (bottom;a representative distance plot from one index). Gray and black regionsin the middle panel represent regions of the chromosome that are,respectively, absent or present in a given PCR compartment. Only theblack regions (haplotyping islands) are covered by sequencing readsthat carry the index for that given physical compartment. Aligned readsare sorted on the basis of their genomic coordinates, and the distancebetween neighboring alignments from the same partition is recorded.A bimodal distribution is observed, with gray regions represented bythe distal (inter-island) subpopulation and the black regions or islandsrepresented by the proximal (intra-island) subpopulation. Breaks betweenthe islands imply that two neighboring islands do not necessarily belongto the same haplotype. A high ratio of the intra-island peak to the interisland peak indicates strong enrichment of the proximal regions of thegenome that are in the same haplotype phase. Representative intra-islandcoverage is shown in the top panel.function of the distance between them (Fig. 4a). For all pairs that wereon the same phasing block, the probability that a pair was phased correctly, again as a function of distance, was plotted (Fig. 4b). Phasingyield (the average genomic distance with 80% of the SNPs sharingthe phasing block) was in megabase-scale blocks, and a minimumaccuracy of 99.8% extended to 50-kb pairwise SNP distances (Fig. 4a,b).Analogous to other reports26, we separated switch errors intotwo categories—long switches and point switches—with the totalswitch error defined as the sum of short- and long-switch error. Longswitches were defined by a large-scale transition from one haplotypeto another (for example, 00000001111111, where each 0 or 1 represents one haplotype, either from the mother or father), whereas pointaDVANCE ONLINE PUBLICATIONNature Genetics

Technical reportsbPhasing 399.9299.9199.90605040302010605040302010Pairwise SNP distance (bp)5 106662 10105 15 1051052 104 5110 21044310 5110 11035106 105 21066105 15 1051052 104 5110 21044310 15 102 1033510 1300 2014 Nature America, Inc. All rights reserved.7005,0010 0,0015 0,0020 0,0025 0,00030,0035 0,0040 0,0045 0,0050 .8499.8299.8080SNPs in correct phase (%)805, 00010 0,0015 0,0020 0,00025,0030 0,0035 0,00040,0045 0,00050,000SNPs sharing phasing block (%)90Phasing accuracy1002aPairwise SNP distance (bp)Figure 4 Summary of phasing results. (a) Phasing yield. The probability that heterozygous SNP pairs are on the same phasing block as a function of thedistance between them. (b) Phasing accuracy. For all pairs that are on the same phasing block, the probability that a pair is phased correctly is plottedas a function of distance. Insets show the linear relationship of the percentage of SNPs sharing a phasing block (a) and the percentage of SNPs incorrect phase (b) versus pairwise SNP distances of 0–50 kb.switches were defined by local phasing errors that did not affect adjacent positions (for example, 000000100000). The ReFHap step had along-switch accuracy of 99.95–99.97% and a point-switch accuracyof 99.3–99.6% for the samples with four lanes of sequencing data.This level of coverage and accuracy could be achieved without anyimputation and population-level data.Phasing accuracy can be improved by removing SNPs that cannotbe unambiguously assigned to one haplotype and singletons, i.e., SNPsthat are covered by only one data point (Table 1 and SupplementaryFig. 9). This filter enables the creation of a high-quality haplotypingbackbone at the expected cost of lower coverage. Optionally, a subsetof missing SNPs can be imputed using data from the 1000 GenomesProject panel27 (Supplementary Fig. 9). The 1000 Genomes Projectdata were only used to perform imputation on SNPs that were notcovered (‘filling imputation’) but not to connect haplotyping blocks(‘stitching imputation’), as the latter results in a much higher longswitch error rate (Supplementary Fig. 10). Therefore, the N50 of theassembled haplotyping blocks does not change after the imputationstep, and there is minimal chance of introducing errors by creatingchimeric blocks from both parents at this step. The final switch accuracy was very high (with point-switch and long-switch accuraciesof 99.75% and 99.96%, respectively), and the major error mode wassingle point switches (Table 1). Depending on accuracy and coverage requirements, certain steps of the analysis pipeline can be eitherincluded or excluded.We also analyzed genome-wide phasing performance as a functionof sequencing depth. Sequencing data were down-sampled and phasedas previously described. Even with as little as 40–60 Gb of sequencingdata, an amount substantially less than that in previously reportedstudies15, about 95% of the SNPs were covered and phased with a longswitch accuracy of 99.96% or higher (Supplementary Table 1).One key motivation for haplotype-resolved genome sequencingrelates to the phasing of de novo mutations and compound heterozygous variations1. On the basis of our results (SupplementaryTable 2), 45 of the 48 de novo mutations previously described andvalidated for the NA12878 cell line28 were successfully phased here,which further demonstrates the advantage of haplotype-resolvedwhole-genome sequencing. Additionally, 33 de novo SNPs observedboth with CPT-seq and LFR (Long Fragment Read)13,15 wereNature GeneticsADVANCE ONLINE PUBLICATIONconcordant with LFR data for NA12878 and were verified with phasing data from the grandchildren NA12886 and NA12885. For eachde novo SNP, we located the phasing block of the SNP in CPT-seq andcompared the phase of the ten neighboring SNPs to that determinedwith the LFR method. We found all of the phases to be concordant.Furthermore, to demonstrate the phasing of compound heterozygousvariants, we examined the SNPs in genes described as putative compound heterozygotes in NA12878 by Kamphans et al.29. Of a total of10 SNP pairs and 1 trio (23 total SNPs), 19 were phased by CPT-seq,covering 6 pairs and the trio. All phased SNPs shared by both methodswere concordant.DISCUSSIONWe present a simple and robust workflow for capturing genome-widehaplotype information in three basic steps—contiguity-preservingtransposition, dilution and PCR (Fig. 2)—with an overall processingtime of less than 3 h. The strength of this platform relies on (i) the ability to transpose DNA with universal primers and index sequences whilemaintaining contiguity, taking full advantage of DNA quality withoutenforcing any size selection or physical constraint; (ii) universal primeramplification, simplifying downstream processes and contributing toa more uniform genome-wide representation as compared to randomprimer–based methods; and (iii) combinatorial indexing, enabling aparallelized, highly partitioned library construction process that createsthousands of indexed libraries useful in dilution haplotyping.Another critical aspect of this platform is the scalability of the combinatorial indexing scheme. As we demonstrated here, 96 indexedtransposition reactions and 96 indexed PCR reactions generated(using only 40 indexed oligonucleotides) 9,216 virtual compartments.The combinatorial indexing scheme minimized both the number ofphysical compartments required and the amounts of reagent used.The method can easily be adapted to a different number of virtualcompartments depending on the application and available sequencing throughput. For example, 48 96 and 384 384 versions ofour workflow are expected to generate approximately 4,600 and147,000 compartments, respectively. Although the method is basedon dilution, the concept of virtual compartments eliminates the needto dilute DNA to subhaploid content for each physical compartment,thereby minimizing the challenges of low-input amplification.

2014 Nature America, Inc. All rights reserved.Technical reportsThis unique feature allows multiple parental copies of the samegenomic region from a given sample to be present in each physicalcompartment as long as they are from different indexed transposition reactions.Unlike other dilution-based haplotyping methods that build upredundancy by whole-genome amplification, our approach transformsa set of subhaploid genomic molecules directly into a library, albeitwith low coverage (i.e., strobed reads; Fig. 3, top) of any moleculewithin each partition. This low coverage is a consequence of severallosses.

haplotype blocks (N50 .4-2.3 Mb). The rapid, scalable and cost-effective workflow could enable haplotype resolution to become routine in human genome sequencing. Most genomic studies thus far ignore the diploid nature of the human genome1. However, the context in which variation occurs on each indi-