HK1246901B - Method and computer system for analyzing genomic dna of organisms - Google Patents

Description

Method and computer system for analyzing genomic DNA of organisms

This application is a divisional application with an application date of April 13, 2012, application number “201280029331.7” (PCT application number PCT/US2012/033686), and invention name “Processing and Analysis of Complex Nucleic Acid Sequence Data”.

Cross-reference to related applications

This application claims the benefit of priority to U.S. Provisional Patent Application No. 61/517,196, filed April 14, 2011, which is hereby incorporated by reference in its entirety.

This application claims the benefit of priority to U.S. Provisional Patent Application No. 61/527,428, filed August 25, 2011, which is hereby incorporated by reference in its entirety.

This application claims the benefit of priority to U.S. Provisional Patent Application No. 61/546,516, filed October 12, 2011, which is hereby incorporated by reference in its entirety.

Background of the Invention

There is a need for improved techniques for analyzing complex nucleic acids, such as methods, in particular, for improving sequence accuracy and for analyzing sequences with substantial errors introduced through nucleic acid amplification.

Furthermore, there is a need for improved techniques for determining parental contributions to the genomes of higher organisms, ie, the human genome. Methods for haplotype phasing, including computational methods and experimental phasing, are reviewed in Browning and Browning, Nature Reviews Genetics 12:703-7014, 2011.

SUMMARY OF THE INVENTION

The present invention provides techniques for analyzing sequence information derived from complex nucleic acid sequencing (as defined herein) that result in unit type phasing, error reduction, and other features based on algorithms and analysis techniques developed in conjunction with long fragment read (LFR) technology.

According to one aspect of the present invention, a method for determining the sequence of a complex nucleic acid (e.g., a whole genome) of one or more organisms (i.e., an individual organism or a population of organisms) is provided. Such a method comprises: (a) receiving a plurality of reads of the complex nucleic acid at one or more computing devices; and (b) generating, using a computing device, an assembled sequence of the complex nucleic acid from the reads, the assembled sequence comprising less than 1.0, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.08, 0.07, 0.06, 0.05, or 0.04 false single nucleotide variants per megabase at a call rate of 70, 75, 80, 85, 90, or 95% or higher, wherein the method is performed by the one or more computing devices. In some aspects, a computer-readable non-transitory storage medium stores one or more sequential instructions comprising instructions that, when executed by one or more computing devices, cause the one or more computing devices to perform the steps of such a method.

According to one embodiment, wherein such a method involves haplotype phasing, the method further comprises identifying a plurality of sequence variants in the assembled sequence and phasing the sequence variants (e.g., 70, 75, 80, 85, 90, 95% or more of the sequence variants) to produce a phased sequence, i.e., a sequence in which the sequence variants are phased. Such phasing information can be used in the context of error correction. For example, according to one embodiment, such a method comprises identifying as an error a sequence variant that is inconsistent with the phasing of at least two (or three or more) phased sequence variants.

According to another such embodiment, in such methods, the step of receiving multiple read results of complex nucleic acids includes a computing device and/or its computer logic that receives multiple read results from each of a plurality of aliquots, each aliquot comprising one or more fragments of complex nucleic acids. Information about the aliquots providing such fragments can be used to correct errors or respond to bases that would otherwise be "no response." According to one such embodiment, such methods include a computing device and/or its computer logic that responds to bases at the positions of the assembly sequence based on preliminary base calls from the positions of two or more aliquots. For example, a method can include responding to bases at a certain position of the assembly sequence based on preliminary base calls from at least two, at least three, at least four, or more than four aliquots. In some embodiments, such methods can include identifying a base call as true if it exists in at least two, at least three, at least four, or more than four aliquots. In some embodiments, such methods can include identifying a base call as genuine if it is present in at least a majority (or at least 60%, at least 75%, or at least 80%) of the aliquots that make preliminary base calls to the position in the assembled sequence. According to another such embodiment, such methods include a computing device and/or computer logic thereof that identifies a base call as genuine if it is present three or more times in reads from two or more aliquots.

According to another such embodiment, the aliquot of the origin of the read result is determined by identifying the aliquot-specific tag (or aliquot-specific tag group) attached to each fragment. Optionally, such aliquot-specific tags include error correction or error detection codes (e.g., Reed-Solomon error correction codes). According to one embodiment of the present invention, after sequencing the fragment and the attached aliquot-specific tags, the resulting read result includes tag sequence data and fragment sequence data. If the tag sequence data is correct, that is, if the tag sequence matches the tag sequence for aliquot identification, or alternatively if the tag sequence data has one or more errors that can be corrected using an error correction code, the read result including such tag sequence data can be used for all purposes, particularly for a first computer method (e.g., performed on one or more computing devices), which requires tag sequence data and produces a first output, including but not limited to unit type phasing, sample multiplexing, library multiplexing, phasing, or any error correction method based on correct tag sequence data (e.g., error correction method based on identifying the origin aliquot of a specific read result). If the tag sequence is incorrect and cannot be corrected, the read results containing such incorrect tag sequence data are not discarded, but are used in a second computer method (e.g., performed by one or more computing devices) that does not require tag sequence data, includes but is not limited to positioning, assembly, and set-based statistics, and produces a second output.

According to another embodiment, such a method further comprises: a computing device and/or computer logic thereof that provides a first phased sequence of a region of a complex nucleic acid, wherein the region comprises short tandem repeats; a computing device and/or computer logic thereof that compares reads of the first phased sequence of the region (e.g., regular or mate-pair reads) with reads of a second phased sequence of the region (e.g., using sequence coverage); and a computing device and/or computer logic thereof that identifies a short tandem repeat expansion in one of the first phased sequence or the second phased sequence based on the comparison.

According to another embodiment, the method further comprises a computing device and/or computer logic thereof that obtains genotype data from at least one parent of the organism and generates an assembled sequence of the complex nucleic acid from the reads and the genotype data.

According to another embodiment, the method further comprises a computing device and/or computer logic thereof that performs steps comprising: aligning a plurality of the reads for a first region of the complex nucleic acid, thereby creating an overlap between the aligned reads; identifying N heterozygous candidates within the overlap; clustering the space of ^2N to ^4N possibilities or a selected subspace thereof, thereby creating a plurality of clusters; identifying two clusters with the highest density, each identified cluster comprising a substantially noise-free center; and repeating the aforementioned steps for one or more additional regions of the complex nucleic acid. The clusters identified for each region can define a contig, and these contigs can be matched to each other to form groups of contigs, one for each haplotype.

According to another embodiment, such methods further comprise providing an amount of a complex nucleic acid, and sequencing the complex nucleic acid to generate a read.

According to another embodiment, in such methods, the complex nucleic acid is selected from the group consisting of a genome, an exome, a transcriptome, a methylome, a mixture of genomes of different organisms, and a mixture of genomes of different cell types of an organism.

According to another aspect of the present invention, an assembled human genome sequence produced by any of the above methods is provided. For example, one or more computer-readable non-transitory storage media stores an assembled human genome sequence produced by any of the above methods. According to another aspect, a computer-readable non-transitory storage medium stores one or more sequences of instructions comprising instructions that, when executed by one or more computing devices, cause the one or more computing devices to perform any, some, or all of the above methods.

According to another aspect of the present invention, a method for determining the sequence of a whole human genome is provided, the method comprising: (a) receiving a plurality of reads of the genome at one or more computing devices; and (b) generating an assembled sequence of the genome from the reads using the one or more computing devices, the assembled sequence comprising less than 600 false heterozygous single nucleotide variants per gigabase at a genome call rate of 70% or greater; according to one embodiment, the assembled sequence of the genome has a genome call rate of 70% or greater and an exome call rate of 70% or greater. In some aspects, a computer-readable non-transitory storage medium stores one or more sequential instructions comprising instructions that, when executed by one or more computing devices, cause the one or more computing devices to perform any of the inventive methods described herein.

According to another aspect of the present invention, a method for determining the sequence of a whole human genome is provided, the method comprising: (a) receiving, at one or more computing devices, a plurality of reads from each of a plurality of aliquots, each aliquot comprising one or more fragments of a genome; and (b) generating, with the one or more computing devices, a phased assembly sequence of the genome from the reads, the assembly sequence comprising less than 1000 false single nucleotide variants per gigabase at a genome call rate of 70% or greater. In some aspects, a computer-readable non-transitory storage medium stores one or more sequential instructions comprising instructions that, when executed by one or more computing devices, cause the one or more computing devices to perform such a method.

BRIEF DESCRIPTION OF THE DRAWINGS

Figures 1A and 1B show examples of sequencing systems.

FIG. 2 shows an example of a computing device that can be used in or in conjunction with a sequencer and/or computer system.

Figure 3 shows the general architecture of the LFR algorithm.

Figure 4 shows pairwise analysis of adjacent heterozygous SNPs.

Figure 5 shows an example of selecting hypotheses and assigning scores to hypotheses.

Figure 6 shows the graph construction.

Figure 7 shows the graph optimization.

Figure 8 shows the contig alignment.

Figure 9 shows parent-assisted universal phasing.

Figure 10 shows the natural contig separation.

Figure 11 shows general phasing.

Figure 12 shows error detection using LFR.

FIG. 13 shows an example of a method for reducing the number of false negatives, wherein confident heterozygous SNP calls can be made regardless of how small the number of reads is.

Figure 14 shows the detection of CTG repeat expansion in human embryos using haplotype-resolved clonal coverage.

15 is a graph showing amplification of purified genomic DNA standards (1.031, 8.25, and 66 picograms [pg]) and 1 or 10 PVP40 cells using the multiple displacement amplification (MDA) protocol, as described in Example 1.

Figure 16 shows the data related to GC preference obtained by amplifying two MDA protocols. The average number of cycles across the entire plate was determined and deducted from each individual marker to calculate the "Δ cycle" number. The GC content of 1000 base pairs around each marker was plotted against the Δ cycle to indicate the relative GC preference of each sample (not shown). The absolute value of each Δ cycle was summed to create a "Δ sum" metric. Lower Δ sums and the relatively flat curve of the data relative to GC content produced a better presented full-genome sequence. The Δ sums were 61 (for our MDA method) and 287 (for SurePlex amplified DNA), indicating that our protocol produces a much smaller GC preference than the SurePlex protocol.

Figure 17 shows the genome coverage of sample 7C and 10C.Use the 10 megabase moving average of the 100 kilobase coverage window of standardization to cover drawing with respect to the haploid genome.The dotted lines at copy number 1 and 3 places represent haploid and triploid copy number respectively.These two embryos are males, and have haploid copy number for X and Y chromosome.In these samples, do not find other loss or acquisition of full chromosome or chromosome large segment.

Figure 18 is a schematic diagram of an embodiment of the barcode adapter design for the method of the present invention. The LFR adapter is composed of a unique 5' barcode adapter, a common 5' adapter and a common 3' adapter. The common adapters are all designed to have 3' dideoxynucleotides that cannot be connected to the 3' fragments, which eliminates the formation of adapter dimers. After connection, the closed portion of the adapter is removed and replaced with an unblocked oligonucleotide. The remaining nick is resolved by subsequently nicking with a Taq polymer and connecting with a T4 ligase.

Figure 19 shows a cumulative GC coverage plot. Cumulative GC coverage was plotted for LFR and standard libraries to compare GC preference differences. For sample NA19240 (a and b), three LFR libraries (replicate 1, replicate 2, and 10 cells) and one standard library were plotted for both the entire genome (c) and coding portion only (d). In all LFR libraries, coverage loss in high GC regions was evident, which was more pronounced in coding regions (b and d) that contained a higher proportion of GC-rich regions.

Figure 20 shows a comparison of haplotype performance between genome assemblies. The variant responses of the standard assembly library and the LFR assembly library were combined and used as loci for phasing, except for the specified cases. The LFR phasing rate is based on the calculation of parental phased heterozygous SNPs. *For those individuals without parental genomic data (NA12891, NA12892 and NA20431), the phasing rate was calculated by dividing the number of phased heterozygous SNPs by the number of heterozygous SNPs expected to be true (the number of SNPs attempted to be phased - 50,000 expected errors). The N50 calculation is based on the total assembly length of all contigs relative to the NCBI component 36 (component 37 in the case of NA1924010 cells and high coverage and NA20431 high coverage) human reference genome. Since all DNA is denatured into single strands and dispersed on a 384-well plate, the haploid fragment coverage is 4 times larger than the number of cells. The insufficient amount of starting DNA explains the lower phasing efficiency in the NA20431 genome. Samples with #10 cells were measured using coverage of wells containing more than 10 cells (which may be a result of these cells being in various stages of the cell cycle at the time of collection). Phasing rates ranged from 84% to 97%.

Figure 21 shows the LFR unit typing algorithm. (a) Variant extraction: Variant is extracted from the read results of the tagged aliquot. The 10-base Reed-Solomon code ensures that the label recovery can be achieved through error correction. (b) Connectivity evaluation of heterozygous SNP pairs: For each heterozygous SNP pair within a certain neighborhood, the matrix of shared aliquots is calculated. Ring 1 is the total heterozygous SNPs on a chromosome. Ring 2 is the total heterozygous SNPs in the neighborhood of the ring 1 heterozygous SNP on the chromosome. This neighborhood is restricted by the expected number of heterozygous SNPs and the expected fragment length. (c) Graph generation: An undirected graph is generated in which nodes correspond to heterozygous SNPs and connections correspond to the orientation and strength of the best hypothesis of the relationship between those SNPs. (As used in this article, a "node" is a data [data item or data object] that can have one or more numerical values, which represent base calls or other sequence variations (such as heterozygosity or indels) in a polynucleotide sequence.) The direction is binary. Figure 21 depicts the flipped and unflipped relationships between pairs of heterozygous SNPs, respectively. The strength is defined by using fuzzy logic operations on the elements of the shared aliquot matrix. (d) Graph optimization: The graph is optimized via minimum spanning tree operations. (e) Overlap generation: Each subtree is simplified into overlapping groups, which is performed by keeping the first heterozygous SNP unchanged and flipping or not flipping the other heterozygous SNPs on the subtree based on their path to the first heterozygous SNP. It is arbitrary to assign parent 1 (P1) and parent 2 (P2) to each overlapping group. Gaps in the full chromosome tree define the boundaries of different subtrees/overlap groups on the chromosome. (f) Mapping LFR overlapping groups to parental chromosomes: Using parental information, mother or father labels are placed on the P1 and P2 unit types of each overlapping group.

Figure 22 shows the haplotype discordance between replicate LFR libraries. The two replicate libraries from samples NA12877 and NA19240 were compared at all shared phased heterozygous SNP loci. This is a comprehensive comparison because most phased loci are shared between the two libraries.

Figure 23 shows the error reduction achieved by LFR. The standard library heterozygous SNP response alone, and the combination with the LFR response, are all independently phased by repeating the LFR library. Generally, LFR introduces about 10 times more false positive variant responses. This is most likely due to the random incorporation of incorrect bases during the multiple displacement amplification based on phi29. Importantly, if the heterozygous SNP response is required to be phased and to be visible in three or more independent wells, the reduction in error is significant, and the result is also better than the standard library without error correction. LFR can also remove errors from the standard library, which improves the response accuracy by about 10 times.

Figure 24 shows the LFR response of the position of no response. In order to prove the potential of LFR to rescue the position of no response, three example positions that were not responded (no response) by standard software were selected on chromosome 18. By phasing them with C/T heterozygous SNPs as part of the LFR overlapping group, these positions can be partially or completely responded. The distribution of shared holes (those holes with at least one reading result for each base in the paired bases; there are 16 pairs of bases for a pair of assessed loci) allows the three N/N positions to be responded to A/N, C/C and T/C responses again, and C-A-C-T and T-N-C-C are defined as unit types. Using the information of the hole allows LFR to accurately respond to those alleles that have only as few as 2-3 reading results (about 3 times less than the case of no hole information) in 2-3 expected holes.

Figure 25 shows the number of genes with multiple adverse variants in each analyzed sample.

Figure 26 shows genes with allelic differential expression and SNPs that alter TFBSs in NA20431. In a non-exhaustive list of genes demonstrated to have significant allelic differential expression, six genes were found to have SNPs that alter TFBSs that correlated with the observed expression differences between alleles. All positions are given relative to NCBI build 37. "CDS" stands for coding sequence, and "UTR3" stands for 3' untranslated region.

Detailed Description of the Invention

As used herein and in the appended claims, the singular forms "a," "an," and "the" include plural referents unless the context clearly dictates otherwise. Thus, for example, reference to "a polymerase" refers to one reagent or a mixture of such reagents, and reference to "the method" includes reference to equivalent steps and/or methods known to those skilled in the art, and so forth.

Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. All publications mentioned herein are incorporated herein by reference for the purpose of describing and disclosing the devices, compositions, formulations, and methods that are described in the publications and that can be used in connection with the presently described methods.

Where a range of values is provided, it is understood that each intervening value, to the tenth of the unit of the lower limit, unless the context clearly dictates otherwise, between the upper and lower limits of the stated range, and any other stated or intervening value in the stated range is encompassed within the invention. The upper and lower limits of these smaller ranges may independently be included, and the smaller ranges are also encompassed within the invention, subject to any express exclusion in the stated range. Where a stated range includes one or both of the limits, ranges excluding either of those included limits are also encompassed within the invention.

In the following description, a number of specific details are listed to provide a more thorough understanding of the present invention. However, it will be apparent to those skilled in the art that the present invention can be practiced without one or more of these specific details. In other cases, well-known features and procedures known to those skilled in the art have not been described to avoid obscuring the present invention.

While the invention has been described primarily with reference to specific embodiments, it is contemplated that other embodiments will become apparent to those skilled in the art upon reading this disclosure, and it is intended that such embodiments be encompassed within the methods of the invention.

Sequencing systems and data analysis

In some embodiments, sequencing of a DNA sample (e.g., such as a sample representing a whole human genome) can be performed by a sequencing system. Two examples of sequencing systems are shown in FIG1 .

Figures 1A and 1B are block diagrams of an example sequencing system 190, which is configured to implement the techniques and/or methods for nucleic acid sequence analysis according to the embodiments described herein. The sequencing system 190 may include or be associated with multiple subsystems, such as, for example, one or more sequencers such as a sequencer 191, one or more computer systems such as a computer system 197, and one or more data repositories such as a data repository 195. In the embodiment shown in Figure 1A, the multiple subsystems of the system 190 may be communicatively connected via one or more networks 193, which may include packet switching or other types of network infrastructure devices (e.g., routers, switches, etc.) configured to facilitate information exchange between remote systems. In the embodiment shown in Figure 1B, the sequencing system 190 is a sequencing device, in which multiple subsystems (e.g., such as a sequencer 191, a computer system 197, and possibly a data repository 195) are components that are coupled for communication and/or operation and integrated within the sequencing device.

In some operational contexts, the data repository 195 and/or computer system 197 of the embodiments shown in Figures 1A and 1B can be configured within a cloud computing environment 196. In a cloud computing environment, storage devices including data repositories and/or computing devices including computer systems can be allocated and instantiated as utilities and used as needed; thus, the cloud computing environment provides as a service infrastructure (e.g., physical and virtual machines, raw/block storage, firewalls, load balancers, aggregators, networks, storage clusters, etc.), a platform (e.g., computing devices and/or solution stacks that can include operating systems, programming language execution environments, database servers, web servers, application servers, etc.), and software necessary to implement any storage-related and/or computing tasks (e.g., applications, application programming interfaces or APIs, etc.).

Note that in various embodiments, the techniques described herein may be implemented by a variety of systems and devices including some or all of the above-described subsystems and components (e.g., such as sequencers, computer systems, and data repositories) in various configurations and form factors; as such, the example embodiments and configurations shown in Figures 1A and 1B should be viewed in an illustrative and non-limiting sense.

Sequencer 191 configurations and can be operated as to accept the target nucleic acid 192 that is derived from biological sample fragment, and target nucleic acid is implemented to order-checking.Any suitable machine that can implement order-checking can be used, wherein such machine can use various sequencing technologies, it includes but not limited to by hybridization order-checking, by connecting order-checking, by synthetic order-checking, single molecule order-checking, optical sequence detection, electromagnetic sequence detection, voltage variation sequence detection and be suitable for producing any other now known or later technology of reading result order-checking result from DNA.In a plurality of embodiments, sequencer can order-check target nucleic acid, and can produce and read result order-checking result, it can comprise or do not comprise breach, and can be or not pairing-to (or paired end) reading result.As shown in Figure 1A and 1B, sequencer 191 order-check target nucleic acid 192, and obtain reading result order-checking result 194, it is transmitted to be stored in one or more data storage repositories 195 and/or processed by one or more computer systems 197 with (temporarily and/or lasting).

The data repository 195 can be executed on one or more storage devices (e.g., hard drives, optical disks, solid-state drives, etc.), which can be configured as disk arrays (e.g., such as SCSI arrays), storage clusters, or any other suitable storage device configuration. The storage devices of the data repository can be configured as internal/integrated components of the system 190 or as external components attachable to the system 190 (e.g., such as external hard drives or disk arrays) (e.g., as shown in FIG. 1B ), and/or can be communicatively interconnected in a suitable manner, such as, for example, a grid, a storage cluster, a storage area network (SAN), and/or a network attached storage (NAS) (e.g., as shown in FIG. 1A ). In various embodiments and implementations, the data repository can be executed on the storage devices as one or more file systems that store information in files, as one or more databases that store information in data records, and/or as any other suitable data storage configuration.

The computer system 197 may include one or more computing devices comprising a general purpose processor (e.g., a central processing unit or CPU), memory, and computer logic 199, which, together with configuration data and/or operating system (OS) software, may implement some or all of the techniques and methods described herein, and/or may control the operation of the sequencer 191. For example, any of the methods described herein (e.g., for error correction, cell type phasing, etc.) may be implemented in whole or in part by a computing device comprising a processor that may be configured to execute logic 199 for implementing the various methods of the method. In addition, although the method steps may be presented as numbered steps, it should be understood that the steps of the methods described herein may be implemented simultaneously (e.g., performed in parallel by a cluster of computing devices) or in a different order. The functionality of the computer logic 199 may be implemented in a single integrated module (e.g., in an integrated logic) or may be combined in two or more software modules that may provide some additional functionality.

In some embodiments, computer system 197 can be a single computing device. In other embodiments, computer system 197 can include multiple computing devices that can communicate and/or be operably interconnected in a grid, cluster, or cloud computing environment. Such multiple computing devices can be configured in different form factors such as computing nodes, blades, or any other suitable hardware configurations. For these reasons, the computer system 197 in Figures 1A and 1B should be viewed in an illustrative and non-restrictive sense.

2 is a block diagram of an example computing device 200 as part of a sequencer and/or computer system that may be configured to execute instructions for implementing various data processing and/or control functionality.

In FIG2 , computing device 200 includes several components interconnected directly or indirectly via one or more system buses, such as bus 275. Such components may include, but are not limited to, keyboard 278, persistent storage device 279 (e.g., such as a fixed disk, solid-state disk, optical disk, etc.), and display adapter 282, to which one or more display devices (e.g., such as an LCD monitor, a flat-panel monitor, a plasma screen, etc.) may be coupled. Peripheral devices and input/output (I/O) devices (which are coupled to I/O controller 271) may be connected to computing device 200 by a variety of means known in the art, including, but not limited to, one or more serial ports, one or more parallel ports, and one or more universal serial buses (USB). External interface 281 (which may include a network interface card and/or a serial port) may be used to connect computing device 200 to a network (e.g., such as the Internet or a local area network (LAN)). External interface 281 may also include a number of input interfaces that may receive information from various external devices, such as, for example, a sequencer or any component thereof. Interconnection via system bus 275 allows one or more processors to communicate with each other. Processor (e.g., CPU) 273 communicates with each connected component and executes instructions (and/or controls their execution) from system memory 272 and/or storage device 279, as well as information exchange between the various components. System memory 272 and/or storage device 279 may be embodied as one or more computer-readable, non-transitory storage media that store sequences of instructions executed by processor 273 and other data. Such computer-readable, non-transitory storage media include, but are not limited to, random access memory (RAM), read-only memory (ROM), electromagnetic media (e.g., such as hard drives, solid-state drives, thumb drives, floppy disks, etc.), optical media such as compact disks (CDs) or digital versatile disks (DVDs), flash memory, etc. Various data values and other structured or unstructured information may be output from one component or subsystem to another, presented to a user via display adapter 282 and a suitable display device, transmitted over a network via external interface 281 to a remote device or remote data repository, or stored (temporarily and/or permanently) on storage device 279.

Any method and functionality implemented by computing device 200 can be performed in a module or integrated manner using hardware and/or computer software in a logical form. As used herein, "logic" refers to a group of instructions that can be operated as the data of implementing one or more functionalities and/or returning one or more result forms or the data used by other logical elements when executed by one or more processors (e.g., CPU) of one or more computing devices. In multiple embodiments and implementations, any given logic can be performed as one or more software components executable by one or more processors (e.g., CPU), as one or more hardware components such as application specific integrated circuits (ASIC) and/or field programmable gate arrays (FPGA), or as any combination of one or more software components and one or more hardware components. The software component of any particular logic can be performed without limitation as a standalone software application, as a client in a client-server system, as a server in a client-server system, as one or more software modules, as one or more function libraries, and as one or more static and/or dynamic link libraries. During execution, the instructions of any particular logic may be embodied as one or more computer processes, threads, fibers, and any other suitable runtime entities, which may be instantiated on the hardware of one or more computing devices and may be allocated computing resources, which may include, but are not limited to, memory, CPU time, storage space, and network bandwidth.

Techniques and algorithms used for the LFR process

Base calling

General methods for sequencing a target nucleic acid using the compositions and methods of the invention are described herein and, for example, in U.S. Patent Application Publication No. 2010/0105052-A1; Published Patent Application Nos. WO2007120208, WO2006073504, WO2007133831, and US2007099208, and U.S. Patent Application Nos. 11/679,124; 11/981,761; 11/981,661; 11/981,605; 11/981,793; 11/981,804; 11/451,691; 11/981,607; 11/981,767; 11/982,467; 11/451,692; 11/541 ,225;11/927,356;11/927,388;11/938,096;11/938,106;10/547,214;11/981,730;11/981,685;11/981,797;11/934,695;11/934,697;11/934,703;1 2/265,593; 11/938,213; 11/938,221; 12/325,922; 12/252,280; 12/266,385; 12/329,365; 12/335,168; 12/335,188; and 12/361,507, which are incorporated herein by reference in their entirety for all purposes. See also Dr. Mana et al., Science 327, 78-81, 2010. Long fragment read (LFR) methods have been disclosed in U.S. Patent Application Nos. 12/816,365, 12/329,365, 12/266,385, and 12/265,593 and U.S. Patent Nos. 7,906,285, 7,901,891, and 7,709,197, which are hereby incorporated by reference in their entireties. Further details and improvements are provided herein.

In some embodiments, data extraction can rely on two types of image data: a bright field image of all DNB positions on the partitioned surface, and a fluorescent image set obtained during each sequencing cycle. Data extraction software can be used to identify all objects with bright field images, and then for each such object, the software can be used to calculate the average fluorescence value of each sequencing cycle. For any given cycle, there are four data points, which correspond to four images taken at different wavelengths to inquire whether the base is A, G, C or T. These raw data points (also referred to as "base calls" in this article) are merged to produce a discontinuous read result sequencing result for each DNB.

A computing device can assemble an identification base group to provide sequence information about a target nucleic acid and/or identify the presence of a specific sequence in the target nucleic acid. For example, a computing device can assemble an identification base group according to the techniques and algorithms described herein by executing various logics; examples of such logic are software code written in any suitable programming language such as Java, C++, Perl, Python, and any other suitable conventional and/or object-oriented programming language. When executed in the form of one or more computer processes, such logic can read results, write, and/or otherwise process structured and unstructured data, which can be stored in various structures on persistent storage and/or in volatile memory; examples of such storage structures include, but are not limited to, files, tables, database records, arrays, lists, vectors, variables, memory and/or processor registers, persistent and/or memory data objects exemplified from object-oriented classes, and any other suitable data structures. In some embodiments, the identified bases are assembled into a complete sequence by comparing overlapping sequences obtained from multiple sequencing cycles implemented on multiple DNBs. As used herein, the term "complete sequence" refers to the sequence of a portion or the entire genome and a portion or the entire target nucleic acid. In other embodiments, the assembly method implemented by one or more computing devices or its computer logic utilizes and can be used for " piecing together " overlapping sequences to provide the algorithm of complete sequence.In other embodiments, reference table is used for assisting the sequence of identifying to be assembled into complete sequence.Can use the reference table of existing sequencing data compilation about selecting organism.For example human genome data can be via National Center for Biotechnology Information at ftp.ncbi.nih.gov/refseq/release, or visit via J.Craig Venter Institute at www.jcvi.org/researchhuref/.The subset of whole human genome information or human genome information can be used for creating the reference table for specific sequencing inquiry.In addition, specific reference table can comprise from the empirical data of specific group, the genetic sequence construction of the mankind with specific ethnicity, geographical tradition, religion or culture limit group, because the variation in human genome can make reference data tilt with the information origin wherein contained. Exemplary methods for responding to variations in a polynucleotide sequence compared to a reference polynucleotide sequence and for polynucleotide sequence assembly (or reassembly) are provided, for example, in U.S. Patent Publication No. 2011-0004413, entitled “Method and System for Calling Variations in a Sample Polynucleotide Sequence with Respect to a Reference Polynucleotide Sequence” (which is incorporated herein by reference for all purposes).

In any embodiment of the invention discussed herein, nucleic acid templates and/or DNB populations can include many target nucleic acids to substantially cover the entire genome or the entire target polynucleotide. As used herein, "substantially covering" means that the amount of nucleotides (i.e., target sequence) analyzed contains an equivalent of at least two copies of the target polynucleotide, or in another aspect, at least 10 copies, or in another aspect, at least 20 copies, or in another aspect, at least 100 copies. The target polynucleotide can include DNA fragments comprising genomic DNA fragments and cDNA fragments and RNA fragments. Guidance for reconstructing the steps of the target polynucleotide sequence can be found in the following references, which are incorporated by reference: Lander et al, Genomics, 2: 231-239 (1988); Vingron et al, J. Mol. Biol., 235: 1-12 (1994); and similar references.

In some embodiments, four images are generated for each interrogation position of the complex nucleotide sequence being sequenced, one for each color dye. The resulting intensity of each of the positions and four colors of each of the positions of each point in the image is determined by adjusting the crosstalk between the dye and the background intensity. The quantitative model can be fitted to the resulting four-dimensional data set. A quality score is used to respond to a given point base, and the quality score reflects how well the four intensities fit the model.

Base calling of four images per field of view can be performed in several steps by one or more computing devices or their computer logic. First, the image intensities are corrected for background using a modified morphological "image open" operation. Since the positions of the DNBs line up with the camera pixel positions, intensity extraction is accomplished as a simple readout of the pixel intensities from the background-corrected images. These intensities are then corrected for several sources of both optical and biological signal crosstalk, as described below. The corrected intensities are then passed to a probabilistic model, which ultimately produces a set of four possible base call outcomes for each DNB. Several metrics are then combined using a pre-fitted logistic regression to calculate a base call score.

Intensity Correction: Several sources of biological and optical crosstalk are corrected using a linear regression model implemented as computer logic executed by one or more computing devices. Linear regression is superior to deconvolution methods, which are computationally more expensive and produce results of similar quality. Sources of optical crosstalk include filter band overlap between the four fluorescent dye spectra and lateral crosstalk between adjacent DNBs due to light diffraction at their close proximity. Biological sources of crosstalk include incomplete washing from previous cycles, probe synthesis errors and probe "slippage" that contaminates signals at adjacent positions, and incomplete anchor extension when interrogating bases "outside" (further away from the anchor) of the anchor. Linear regression is used to determine the portion of the DNB intensity that can be predicted using the intensity of any adjacent DNB or the intensity from previous cycles or other DNB positions. The portion of the intensity that can be explained by these sources of crosstalk is then subtracted from the initial extracted intensity. To determine the regression coefficient, the intensity on the left side of the linear regression model needs to consist primarily of "background" intensity, that is, the intensity of DNBs that will not respond to the given base being regressed. This requires a pre-calling step using the initial intensities. Once a DNB is selected that does not have a specific base call (with reasonable confidence), the computing device or its computer logic performs a simultaneous regression of the sources of crosstalk:

Neighbor DNB crosstalk was corrected using the above regression. Additionally, each DNB was corrected for its specific neighborhood using a linear model involving all neighbors in all available DNB positions.

Base Call Probability: Calling a base with maximum intensity does not result in different shapes for the background intensity distributions of the four bases. To account for these possible differences, a probability model is developed based on an empirical probability distribution of background intensities. Once the intensities are corrected, the computing device or its computer logic pre-calls some DNBs (DNBs that pass a certain confidence threshold) with maximum intensity and uses these pre-call DNBs to drive the background intensity distribution (the intensity distribution of DNA that does not call a given base). After obtaining this distribution, the computing device can calculate the tail probability under the distribution for each DNB, which describes the empirical probability that the intensity is the background intensity. Therefore, for each DNB and each of the four intensities, the computing device or its logic can obtain and store its probability of being background. The computing device can then calculate the probabilities of all possible base calls using these probabilities. Possible base call outcomes need to also describe points that can be doubly or generally multiply occupied by DNBs, or not occupied by DNBs. Combining the calculated probabilities with their prior probabilities (lower priors for multiply occupied or empty points) yields the probabilities of 16 possible outcomes:

These 16 probabilities can then be combined to obtain a reduced set of four probabilities for the four possible base calls. That is:

Score calculation: Logistic regression is used to obtain the score calculation formula. A computing device or its computer logic fits a logistic regression to the positioning results of the base calls using several metrics as input. The metrics include the probability ratio between the call base and the next highest base, the intensity of the call base, an indicator variable for the identity of the call base, and a metric describing the overall clustering quality of the field. All metrics are converted to the log-odds-ratio between the coordinated and uncoordinated calls as collinear. Cross-validation is used to improve the model. A logit function with the final logistic regression coefficient is used to calculate the resulting score.

Positioning and assembly

In another embodiment, the read result data is encoded in a compressed binary form and includes both the base of the call and a quality score. The quality score is associated with the base accuracy. Analysis software logic, including sequence assembly software, can use the score to determine the contribution of evidence from each base with a read result.

Reads can be "gapped" due to DNB structure. Gap size varies (typically +/- 1 base) due to the inherent variability of enzyme digestion. Due to the random access nature of cPAL, reads can occasionally have unread bases ("no calls") in otherwise high-quality DNBs. Read pairs are paired.

The positioning software logic that can compare read result data and reference sequence can be used for the data location produced by sequencing method described herein.When being executed by one or more computing devices, this type of positioning logic generally can allow small changes relative to reference sequence, such as those changes caused by each genome variation, read result error or not reading result base.This characteristic often allows directly rebuilding SNP.In order to support the larger variation of transfer, including large-scale structural variation or dense variant zone, each arm of DNB can be separately positioned, and spouse (mate) pairing constraint is applied after comparison.

As used herein, the term "sequence variant" or simply "variant" includes any variant, including but not limited to substitution or replacement of one or more bases; insertion or deletion of one or more bases (also known as "indels"); inversions; transitions; duplications or copy number variations (CNVs); trinucleotide repeat expansions; structural variations (SVs; e.g., intrachromosomal or interchromosomal rearrangements, e.g., translocations); and the like. In a diploid genome, "heterozygosity" or "het" is two different alleles of a particular gene in a gene pair. The two alleles can be different mutants or wild-type alleles paired with a mutant. The present method can also be used in analyzing non-diploid organisms, whether such organisms are haploid/haploid (N=1, where N=the haploid number of chromosomes) or polyploid or aneuploid.

In some embodiments, assembly of sequence read results can utilize software logic that supports DNB read result structures (paired, gapped read results with non-responsive bases) to produce diploid genome assemblies, which in some embodiments can be used to generate sequence information for the LFR method of the present invention for phasing heterozygous sites.

The method of the present invention can be used to reconstruct new segments that are not present in the reference sequence. In some embodiments, the following algorithm can be used, which utilizes a combination of evidence (Bayesian) reasoning and an algorithm based on a de Bruijin graph. In some embodiments, a statistical model that is empirically calibrated for each data set can be used, allowing all read result data to be used without pre-filtering or data trimming. Large-scale structural changes (including but not limited to deletions, translocations, etc.) and copy number changes can also be detected by adjusting the paired read results.

Phasing LFR data

Figure 3 depicts the main steps in LFR data phasing. These steps are as follows:

(1) Graph Construction Using LFR Data : One or more computing devices or their computer logic generate an undirected graph in which vertices represent heterozygous SNPs and edges represent connections between those heterozygous SNPs. Edges are composed of directions and connection strengths. One or more computing devices can store such graphs in storage structures including, but not limited to, files, tables, database records, arrays, lists, vectors, variables, memory and/or processor registers, persistent and/or memory data objects instantiated from object-oriented classes, and any other suitable transient and/or persistent data structures.

(2) Graph construction using mate pair data : Step 2 is similar to step 1, where the connection is based on the mate pair data instead of the LFR data. To perform the connection, a DNB must be found with two heterozygous SNPs of interest in the same read (same arm or mate arm).

(3) Graph combination: The computing device or its computer logic representation of each of the above graphs is performed via an NxN sparse matrix, where N is the number of candidate heterozygous SNPs on the chromosome. Two nodes can have only one connection in each of the above methods. In the case of combining the two methods, two nodes can have up to two connections. Therefore, the computing device or its computer logic can use a selection algorithm to select a connection as the selected connection. For these studies, it was found that the quality of the mate pair data was significantly inferior to the quality of the LFR data. Therefore, only the connections derived from LFR were used.

(4) Graph pruning: A series of heuristics are designed and applied by a computing device to the stored graph data to remove some erroneous connections. More precisely, a node must satisfy the condition of at least two connections in one direction and one connection in the other direction; otherwise, it is eliminated.

(5) Graph Optimization: The computing device or its computer logic optimizes the graph by generating a minimum spanning tree (MST). The power function is set to -|strength|. During this process, lower strength edges are eliminated when possible due to competition with stronger paths. Thus, the MST provides a natural selection of the strongest and most reliable connections.

(6) Contig Creation: Once the minimum spanning tree is generated and/or stored in a computer-readable medium, the computing device or its logic may reorient all nodes, obtaining a constant number of nodes (here, the first node). This first node is the anchor node. For each node, the computing device then searches for a path to the anchor node. The orientation of the test node is the aggregate of the orientations of the edges on the path.

(7) Universal phasing: After the above steps, the computing device or its logic phases each of the contigs created in the previous steps. Here, in contrast to phasing, this part of the results is called pre-phased, indicating that this is not the final phasing. Since the first node is arbitrarily chosen as the anchor node, the phasing of the entire contig does not have to be consistent with the parent chromosome. For universal phasing, several heterozygous SNPs for which trio information is available on the contig are used. These trios of heterozygous SNPs are then used to identify the alignment of the contigs. At the end of the universal phasing step, all contigs have been appropriately labeled and can therefore be considered to be whole chromosome contigs.

Contig generation

In order to produce overlapping groups, for each heterozygous SNP pair, computing device or its computer logic tests two hypotheses: forward direction and reverse direction.Forward direction means that two heterozygous SNPs are connected with the same direction that they were originally listed (initially in alphabetical order).Reverse direction means that two heterozygous SNPs are connected with the reverse order of their initial listing.Fig. 4 depicts the paired analysis to adjacent heterozygous SNPs, which involves classifying forward and reverse directions into heterozygous SNP pairs.

Each direction can have digital support, has shown the validity of corresponding hypothesis.This support is the function of 16 units of the connection matrix shown in Figure 5, and this Figure 5 shows the example of hypothesis selection, and assigns score to it.In order to simplify the function, 16 variables are simplified into 3: power (power) 1, power 2 and impurity (impurity). Power 1 and power 2 are the two highest value units corresponding to each hypothesis. Impurity is the ratio of the sum of all other units (rather than 2 corresponding to the hypothesis) and the sum of the units in the matrix. The selection between two hypotheses is carried out based on the sum of the corresponding units. The hypothesis with higher sum is the winning hypothesis. The following calculation is only used to distribute the intensity of the hypothesis. Strong hypothesis is the hypothesis with high numerical value for power 1 and power 2 and low numerical value for impurity.

The three metrics Power 1, Power 2, and Impurity are fed into a fuzzy inference system ( FIG. 6 ) to reduce their effects to a single numerical value—a score—between 0 and 1 (inclusive). A fuzzy inference system (FIS) is implemented as computer logic that can be executed by one or more computing devices.

The linking operation was performed for each heterozygous SNP pair within a reasonable distance up to the expected contig length (eg, 20-50 Kb).Figure 6 shows the graph construction, depicting some exemplary links and intensities for three neighboring heterozygous SNPs.

The rules of the fuzzy inference engine are defined as follows:

(1) If power 1 is small and power 2 is small, the score is very small.

(2) If power 1 is medium and power 2 is small, the score is small.

(3) If power 1 is medium and power 2 is medium, then the score is medium.

(4) If power 1 is large and power 2 is small, the score is medium.

(5) If power 1 is large and power 2 is medium, the score is large.

(6) If power 1 is large and power 2 is large, the score is very large.

(7) If the impurity is smaller, the score is larger.

(8) If the impurity is medium, the score is small.

(9) If the impurities are large, the score will be very small.

For each variable, the definitions of small, medium, and large are different and determined by its specific membership function. After exposing the fuzzy inference system (FIS) to each set of variables, the input set's contribution to the rules is propagated to the fuzzy logic system, and a single (defuzzified) number is generated as output: the score. This score is limited to between 0 and 1, with 1 indicating the highest quality.

After applying FIS to each node, computing device or its computer logic builds whole picture. Fig. 7 shows the example of this figure. Node is colored according to the direction of winning hypothesis. The intensity of each connection is derived by applying FIS to the heterozygous SNP of interest. Once constructing preliminary figure (top figure of Fig. 7), computing device or its computer logic optimizes this figure (bottom figure of Fig. 7) and simplifies it into tree. This optimization process is completed by generating minimum span tree (MST) from initial figure. MST guarantees the unique path from each node to any other node.

FIG7 illustrates graph optimization. In this application, the first node on each contig is used as an anchor node, and all other nodes are oriented relative to the node. Depending on the orientation, each hit will have to flip or flip to match the orientation of the anchor node. FIG8 illustrates the contig alignment method for a given example. At the end of this method, phased contigs are obtained.

At this point in the quantitative method, the two haplotypes are separated. Although it is known that one of these haplotypes is from template, and one is from father, it is not known which one is from which parent. In the next step of phasing, a computing device or its computer logic attempts to classify the correct parental tag (maternal/father) into each haplotype. This process is called universal phasing. In order to do this, it is necessary to know the connection between at least several heterozygous SNPs (on overlapping groups) and parents. This information can be obtained by performing a three-person group (maternal-father-offspring) phasing. Using a triple sequencing genome, some loci with known parental connections are identified, more specifically when at least one parent is homozygous. Then, a computing device or its computer logic uses these connections to classify the correct parental tag (maternal/father) into the entire overlapping group, that is, to implement parent-assisted universal phasing (Fig. 9).

To ensure high accuracy, the following were implemented: (1) when possible (e.g., in the case of NA19240), triplet information was obtained from multiple sources (e.g., in-house and 1000 Genomes) and combinations of such resources were used; (2) contigs were required to contain at least two known triply phased loci; (3) contigs with a series of triplet mismatches in a row (indicating a segment error) were eliminated; and (4) contigs with a single triplet mismatch at the end of a triplet locus (indicating a potential segment error) were eliminated.

Figure 10 shows natural contig separation. Regardless of whether parental data is used, contigs often discontinue beyond a certain point in nature. Reasons for contig separation are: (1) greater than normal DNA fragmentation or lack of amplification in certain regions, (2) low heterozygous SNP density, (3) poly-N sequences on the reference genome, and (4) DNA repetitive regions (which are prone to mislocalization).

Figure 11 shows universal phasing. One of the main advantages of universal phasing is the ability to obtain complete chromosome "contigs." This is possible because each contig (after universal phasing) carries a haplotype with the correct parental label. Therefore, all contigs carrying the label maternal can be placed on the same haplotype; and similar operations can be performed for paternal contigs.

Another major advantage of the LFR method is the ability to significantly improve the accuracy of heterozygous SNP responses. Figure 12 shows two examples of error detection derived from the LFR method. Figure 12 (left side) shows the first example, in which the connection matrix does not support the hypothesis of any expectation. This indicates that one of heterozygous SNPs is not actually a heterozygous SNP. In this example, the A/C heterozygous SNP is actually a homozygous locus (A/A), which is marked as a heterozygous locus by the assembler error. This error can be identified and eliminated or (in this case) corrected. Figure 13 (right side) shows the second example, in which the connection matrix of this situation supports both hypotheses simultaneously. This is the unreal sign of heterozygous SNPerozygous response.

A "healthy" heterozygous SNP linkage matrix is one with only two high cells (at the expected heterozygous SNP positions, i.e., not on a straight line). All other possibilities point to potential problems and can be eliminated or used to generate alternate base calls for the locus of interest.

Another advantage of the LFR method is the ability to respond to heterozygous SNPs with weaker support (e.g., where it is difficult to locate DNBs due to bias or mismatch rates). Since the LFR method requires additional constraints on heterozygous SNPs, the threshold that heterozygous SNP responses require in non-LFR assemblers can be lowered. Figure 13 shows an example of this situation, in which confident heterozygous SNP responses can be made despite a small amount of reads. In Figure 13 (right side), under normal circumstances, a low number of supportive reads would prevent any assembler from confidently responding to the corresponding heterozygous SNPs. However, since the connection matrix is "clean," heterozygous SNP responses can be more confidently attributed to these loci.

Annotation of SNPs in splice sites

Introns in transcribed RNA need to be spliced out before they become mRNA. Information about splicing is embodied in the sequences of these RNAs and is based on consistency. Mutations in the splice site consensus sequence are the cause of many human diseases (Faustino and Cooper, Genes Dev. 17: 419-437, 2011). Most splice sites conform to simple consensus sequences at fixed positions around the exons. At this point, a program for annotating splice site mutations has been developed. In this program, a consensus splice position model (www.life.umd.edu/labs/mount/RNAinfo) is used. A search is performed on the pattern: CAG|G ("|" indicates the start of the exon) in the 5' end region of the exon and MAG|GTRAG ("|" indicates the end of the exon) in the 3' end region of the same exon. Here, M = {A, C}, R = {A, G}. In addition, the splicing consensus positions were classified into two categories: Type I, where consistency with the model is 100% required; and Type II, where consistency with the model is maintained in greater than 50% of cases. It is speculated that SNP mutations in Type I positions will cause missed splicing, while SNPs in Type II positions will only reduce the efficiency of the splicing event.

The program logic for annotating splice site mutations includes two parts. In part I, a file containing a model position sequence from the input reference genome is generated. In part 2, the SNPs from the sequencing project are compared with these model position sequences and any type I and type II mutations are reported. The program logic is exon-centric, replacing the intron-centric (for ease of analyzing the genome). For a given exon, in its 5' end, we look for a consensus "cAGg" (for positions -3, -2, -1, 0. 0 means the beginning of the exon). Uppercase letters mean type I positions, while lowercase letters mean type II positions). In the 3' end of the exon, a search is performed for a consensus "magGTrag" (for position sequences -3, -2, -1, 0, 1, 2, 3, 4). Exons released from genomes that do not meet these requirements are only ignored (approximately 5% of all cases). These exons fall into other minor categories of consensus splice sites and are not investigated by the program logic. Any SNP from the sequenced genome is compared with the model sequence at these genomic positions. Any mismatch in type I is reported. Mismatches in type II positions are reported if the mutation deviates from consensus.

The program logic described above detects most bad splice site mutations. Reported bad SNPs are undoubtedly problematic. However, there are many other bad SNPs that cause splicing problems that this program doesn't detect. For example, the human genome contains many introns that don't conform to the consistency mentioned above. Furthermore, branch point mutations within introns can also cause splicing problems. These splice site mutations are not reported.

Annotate SNPs that affect transcription factor binding sites (TFBSs). The JASPAR model was used to find TFBSs from the released human genome sequence (build 36 or build 37). JASPAR Core is a collection of matrix-modeled TFBS position frequency data for 130 vertebrates (Bryne et al., Nucl. Acids Res. 36: D102-D106, 2008; Sandelin et al., Nucl. Acids Res. 23: D91-D94, 2004). These models were downloaded from the JASPAR website ( http://jaspar.genereg.net/cgi-bin/jaspar_db.pl?rm=browse&db=core&tax_group=vertebrates ). These models were converted into position weight matrices (PWMs) using the following formula: wi = log2[(fi + pNi1/2)/(Ni + Ni1/2)/p], where: fi is the frequency of observation for a particular base at position 1; Ni is the total number of observations at that position; and p is the background frequency of the current nucleotide, which defaults to 0.25 (bogdan.org.ua/2006/09/11/position-frequency-matrix-to-position-weight-matrix-pfm2pwm.html; Wasserman and Sandelin, Nature Reviews, Genetics 5: P276-287, 2004). A specialized program, Mast (meme.sdsc.edu/meme/mast-intro.html), was used to search for TFBS sites within a sequence segment within a genome. The program was run to extract TFBS sites in the reference genome. The steps are summarized as follows: (i) For each gene with an mRNA, extract the [-5000, 1000] putative TFBS-containing regions from the genome, where 0 is the mRNA start position. (ii) Run a Mast search of the putative TFBS-containing sequences for all PWM models. (iii) Select those hits above a given threshold. (iv) For regions with multiple or overlapping hits, select only the single hit, i.e., the hit with the highest Mast search score.

By means of the TFBS model hit with reference to the genome generated and/or stored in a suitable computer readable medium, a computing device or its computer logic can identify the SNPs located in the hit region. These SNPs can affect the model and the hit score changes. A second program is written to calculate this type of change in hit score because the segment containing the SNP is run into the PWM model twice, once for reference, and a second time for the segment with the SNP replacement. SNPs that cause the segment hit score to decline by more than 3 are accredited as bad SNPs.

Selection of genes with two bad SNPs. Genes with bad SNPs were classified into two categories: (1) those affecting the AA sequence of transcription; and (2) those affecting the transcription binding site. For AA sequence effects, the following SNP subcategories were included:

(1) Nonsense or non-stop mutations . These mutations result in either a truncated protein or an extended protein. In either case, the function of the protein product is completely lost or less efficient.

(2) Splice site variants . These mutations cause the splice sites of introns to be destroyed (for those positions that require a certain nucleotide to be 100% according to the model) or severely reduced (for those positions that require a certain nucleotide to be greater than 50% according to the model. SNPs cause the splice site nucleotide to mutate to another nucleotide that is less than 50% identical, as predicted by the splice site consensus sequence model). These mutations have the potential to produce proteins that are truncated, lack exons, or have a severely reduced amount of protein product.

(3) Polyphen2 annotation of AA variation . For SNPs that cause changes in the protein amino acid sequence rather than its length, Polyphen2 (Adzhubei et al., Nat. Methods 7:248-249, 2010) was used as the main annotation tool. Polyphen2 annotated SNPs as "benign", "unknown", "possibly damaging" and "probably damaging". Both "possibly damaging" and "probably damaging" were identified as bad SNPs. These classification assignments by Polyphen2 were based on the structural predictions of the Polyphen2 software.

For transcription binding site mutations, a model maximum score (maxScore) of 75% was used based on the reference genome as a screen for TFBS binding sites. Any model hits with a maximum score of <= 75% in the region were removed. For those remaining hits, if the SNP caused the hit score to drop by more than 3, it was considered a deleterious SNP.

Two categories of genes are reported. Category 1 genes are those with at least 2 bad AA-affecting mutations. These mutations can be all on a single allele (category 1.1) or spread across 2 unique alleles (category 1.2). Category 2 genes are a superset of the category 1 set. Category 2 genes are genes that contain at least 2 bad SNPs, regardless of whether it is AA-affecting or TFBS site-affecting. However, the requirement is that at least 1 SNP is AA-affecting. Category 2 genes are those in category 1, or those with 1 deleterious AA mutation and 1 or more deleterious TFBS-affecting variants. Category 2.1 means that all of these deleterious mutations come from a single allele, while category 2.2 means that the deleterious SNPs come from two unique alleles.

The foregoing techniques and algorithms are applicable to methods for sequencing complex nucleic acids, optionally in combination with LFR processing prior to sequencing (LFR combined with sequencing may be referred to as "LFR sequencing"), which are described in detail below. Such methods for sequencing complex nucleic acids may be implemented by one or more computing devices that execute computer logic. An example of such logic is software code written in any suitable programming language such as Java, C++, Perl, Python, and any other suitable conventional and/or object-oriented programming language. When executed in the form of one or more computer processes, such logic may read results, write, and/or otherwise process structured and unstructured data, which may be stored in a plurality of structures on persistent storage and/or in volatile memory; examples of such storage structures include, but are not limited to, files, tables, database records, arrays, lists, vectors, variables, memory and/or processor registers, persistent and/or memory data objects instantiated from object-oriented classes, and any other suitable data structures.

Improving accuracy in long-read sequencing

In DNA sequencing using certain long-read technologies (e.g., nanopore sequencing), long (e.g., 10-100 kb) read lengths are available, but generally have higher false negative and false positive rates. The ultimate accuracy of sequences from such long-read technologies can be significantly enhanced using cell type information (complete or partial phasing) according to the following general approach.

First, a computing device or its computer logic compares the reads to each other. A large number of heterozygous calls are expected to be present in the overlap. For example, if two to five 100 kb fragments overlap by at least 10%, this results in a >10 kb overlap, which can roughly translate into 10 heterozygous loci. Alternatively, each long read is compared to a reference genome, which implicitly provides a multiple alignment of the reads.

Once multiple read results are aligned, overlaps can be considered. The fact that the overlap includes a large number of heterozygous loci (e.g., N=10) can be adjusted to account for heterozygous combinations. This combination results in a larger space of haplotype probabilities (4N or 4^ ^N ; if N=10, then ^4N =about 1 million). Of all these ^4N points in N-dimensional space, only two points are expected to contain biologically feasible information, i.e., those corresponding to two haplotypes. In other words, there is a noise suppression rate of ^4N /2 (here 1e6/2 or about 500,000). In fact, most of this ^4N space is degenerate, especially because the sequences are already aligned (and therefore similar), and also because each locus generally does not carry more than 2 possible bases (if it is truly heterozygous). Therefore, the lower bound of this space is actually ^2N (if N=10, then ^2N =about 1000). Therefore, the noise suppression rate can be only ^2N /2 (here 1000/2=500), which is still quite impressive. As the number of false positives and false negatives increases, the size of the space expands from ^2N to ^4N , which in turn leads to higher noise suppression rates. In other words, as noise increases, it is automatically suppressed more. Therefore, the expected output product retains only a very small (and fairly constant) amount of noise, almost independent of the input noise. (The tradeoff is the yield loss in noisier conditions). Of course, these suppression rates are changed if: (1) the errors are systematic (or other data idiosyncrasies), (2) the algorithm is not optimal, (3) the overlap is short, or (4) the coverage redundancy is small. N can be any integer greater than 1, such as 2, 3, 5, 10, or more.

The following methods can be used to improve the accuracy of long-read sequencing methods that may have a large initial error rate.

First, a computing device or its computer logic compares several reads, for example, 5 reads or more, such as 10-20 reads. Assume the reads are approximately 100 kb and the shared overlap is 10%, resulting in a 10 kb overlap among the 5 reads. Also assume there is heterozygosity every 1 kb. Therefore, there will be a total of 10 heterozygosities in this common region.

Then, computing device or its computer logic fills in part (for example, only non-zero elements) or the entire matrix of the alpha10 possibilities (wherein alpha is between 2 and 4) of the above-mentioned 10 candidate heterozygosities. In one implementation, only 2 of the alpha10 units of this matrix should be high density (for example, as measured by a threshold value, which can be predetermined or dynamic). These are the units corresponding to true heterozygosity. These two units can be considered to be essentially noiseless centers. The remainder should contain almost 0 and occasionally 1 affiliation, especially when the error is not systematic. If the error is systematic, there can be a clustering event (for example, having more than only a 3rd unit of 0 or 1), which makes the task more difficult. However, even in this case, the cluster affiliation of a false cluster should be significantly weaker than the cluster affiliation of two expected clusters (for example, as measured by absolute or relative amounts). The trade-off in this case is that the starting point should include more multiple sequences of comparisons, which is directly related to having a longer reading result or larger coverage redundancy.

The above steps assume that two feasible clusters are observed between overlapping read results. For a large number of false positives, the situation will not be like this. If this is the case, in alpha dimensional space, the two expected clusters will be blurred, that is, instead of being a single point with high density, they will be the blurred clusters of M points around the unit of interest, where these units of interest are the noiseless centers at the center of the cluster. This enables the clustering method to capture the position of the expected point, despite the fact that the accurate sequence is not presented in each read result. Cluster events can also occur when the cluster is fuzzy (i.e., more than two centers can be arranged), but in a similar manner to the description above, for diploid organisms, scores (such as the total count of cluster units) can be used to distinguish between weaker clusters and two real clusters. Two real clusters can be used to create overlapping groups to multiple regions, as described herein, and overlapping groups can be matched in two groups to form unit types to the larger region of complex nucleic acids.

Finally, the computing device or its computer logic can use population-based (known) haplotypes to increase confidence and/or provide additional guidance in finding true clusters. One way to implement this approach is to assign a weight to each observed haplotype and assign a smaller, but non-zero, value to unobserved haplotypes. By doing so, a preference is achieved for natural haplotypes that have been observed in the population of interest.

Read results using tag sequence data with uncorrected errors

As discussed herein, according to one embodiment of the present invention, a sample of complex nucleic acid is divided into multiple aliquots (e.g., wells in a multi-well plate), amplified, and fragmented. Aliquot-specific tags are then connected to the fragments to identify the aliquots from which specific fragments of the complex nucleic acid originate. Optionally, the tags contain error correction codes, such as Reed-Solomon error correction (or error detection) codes. When sequencing the fragments, both the tags and the fragments of the complex nucleic acid sequence are sequenced. If there is an error in the tag sequence and it is impossible to identify the aliquots from which the fragments originate, or if the sequence is corrected using an error correction code, the entire read result can be abandoned, resulting in a large amount of sequence data loss. It should be noted that the read result containing correct and corrected tag sequence data is highly accurate, but has a low yield, while the read result containing tag sequence data that cannot be corrected is low accurate, but has a high yield. Instead, such sequence data is used for methods different from those that require such data to identify the aliquots of origin by relying on the combination of specific tags with specific aliquots. Examples of methods requiring reads with correct (or corrected) tag sequence data include, but are not limited to, sample or library multiplexing, phasing or error correction, or any other method requiring correct (or corrected) tag sequences. Examples of methods that can employ reads with uncorrected tag sequence data include any other method, including, but not limited to, mapping, reference-based and local reassembly, and set-based statistics (e.g., allele frequencies, positions of de novo mutations, etc.).

Convert long reads into virtual LFRs

Algorithms designed for LFR (including phasing algorithms) can be used for long read results by assigning random virtual tags (with consistent distribution) to each (10-100kb) long fragment. Virtual tags have the benefit of enabling a truly consistent distribution for each code. Due to the differences in the merged codes and the differences in the decoding efficiency of the codes, LFR cannot achieve this consistency level. A ratio of 3:1 (and up to 10:1) can be easily observed in the representation of any two codes in LFR. However, the virtual LFR method results in a true 1:1 ratio between any two codes.

Methods for sequencing complex nucleic acids

Overview

According to one aspect of the present invention, methods for sequencing complex nucleic acids are provided. According to certain embodiments of the present invention, methods for sequencing very small amounts of such complex nucleic acids (e.g., 1 pg to 10 ng) are provided. Even after amplification, such methods produce assembled sequences characterized by high response rates and accuracy. According to other embodiments, aliquot sampling is used to identify and eliminate errors in the sequencing of complex nucleic acids. According to another embodiment, LFR is used in combination with complex nucleic acid sequencing.

Unless otherwise indicated, practice of the present invention can adopt the conventional techniques and descriptions of organic chemistry, polymer technology, molecular biology (including recombinant technology), cell biology, biochemistry and immunology, which are within the skill of the art. Such conventional techniques include polymer array synthesis, hybridization, connection and use of markers to detect hybridization. Specific illustrations of suitable techniques can be provided by reference to the examples below. However, other equivalent conventional methods can certainly also be used. Such routine techniques and descriptions can be found in standard laboratory manuals, such as Genome Analysis: A Laboratory Manual Series (Vols. I-IV), Using Antibodies: A Laboratory Manual, Cells: A Laboratory Manual, PCR Primer: A Laboratory Manual, and Molecular Cloning: A Laboratory Manual (all from Cold Spring Harbor Laboratory Press), Stryer, L. (1995) Biochemistry (4th Ed.) Freeman, New York, Gait, “Oligonucleotide Synthesis: A Practical Approach” 1984, IRL Press, London, Nelson and Cox (2000), Lehninger, Principles of Biochemistry 3rd Ed., W.H. Freeman Pub., New York, N.Y. and Berg et al. (2002) Biochemistry, 5th Ed., W.H. Freeman Pub., New York, N.Y., all of which are incorporated herein by reference in their entirety for all purposes.

General methods for sequencing a target nucleic acid using the compositions and methods of the invention are described herein and in, for example, U.S. Patent Publication Nos. 2010/0105052 and 2007099208 and U.S. Patent Application Nos. 11/679,124 (published as US 2009/0264299); 11/981,761 (US 2009/0155781); 11/981,661 (US 2009/0005252); 11/981,605 (US 2009/0011943); 11/981,793 (US 2009-0118488); 11/451,691 (US 2007/0099208); 11/981,607 (US 2008/0234136); 11/981,767 (US 2009/0137404); 11/982,467 (US 2009/0137414); 11/451,692 (US 2007/0072208); 11/541,225 (US 2010/0081128; 11/927,356(US2008/0318796); 11/927,388(US 2009/0143235); 11/938,096(US 2008/0213771); 11/938,106(US 2008/0171331);10/547,214(US 2007/0037152); 11/981,730 (US 2009/0005259); 11/981,685 (US 2009/0036316); 11/981,797 (US 2009/0011416); 11/934,695 (US 2009/0075343); 11/934,697 (US 2009/0111705); 11/934,703 (US 2009/0111706); 12/265,593 (US 2009/0203551); 11/938,213 (US 2009/0105961);11/938,221(US 2008/0221832); 12/325,922 (US 2009/0318304); 12/252,280 (US 2009/0111115); 12/266,385 (US 2009/0176652); 12/335,168 (US 2009/0311691); 12/335,188 (US 2009/0176234); 12/361,507 (US 2009/0263802), 11/981,804 (US 2011/0004413); and 12/329,365; published international patent application numbers WO2007120208, WO2006073504, and WO2007133831, all of which are incorporated herein by reference in their entirety for all purposes. Exemplary methods for responding to variations in a polynucleotide sequence compared to a reference polynucleotide sequence and for polynucleotide sequence assembly (or reassembly) are provided, for example, in U.S. Patent Publication No. 2011-0004413, (App. No. 12/770,089), which is incorporated herein by reference in its entirety for all purposes. See also Drmanac et al., Science 327, 78-81, 2010. Co-pending related application Nos. 61/623,876, entitled "Identification Of DNA Fragments And Structural Variations", is also incorporated by reference in its entirety and for all purposes.

This method includes extracting and fragmenting the target nucleic acid from the sample. The fragmented nucleic acid is used to generate a target nucleic acid template, which generally includes one or more adapters. The target nucleic acid template is subjected to an amplification method to form a nucleic acid nanoball, which is generally arranged on a surface. The nucleic acid nanoballs of the present invention are subjected to sequencing applications, generally via sequencing by ligation technology, including a combined probe anchor ligation ("cPAL") method, which is described in more detail below. cPAL and other sequencing methods can also be used to detect specific sequences, such as single nucleotide polymorphisms ("SNPs") in the nucleic acid constructs of the present invention (which include nucleic acid nanoballs and linear and circular nucleic acid templates). The above-mentioned patent application and the cited article by Drmanac et al. provide additional detailed information on the following: for example, preparing nucleic acid templates, including adapter design, inserting adapters into genomic DNA fragments to generate circular library constructs; amplifying such library constructs to generate DNA nanoballs (DNBs); generating arrays of DNBs on a solid support; cPAL sequencing; etc., which are used in combination with the methods disclosed herein.

As used herein, the term "complex nucleic acid" refers to a large population of different nucleic acids or polynucleotides. In certain embodiments, the target nucleic acid is genomic DNA; exome DNA (a subset of whole-genome DNA enriched for transcribed sequences, containing a set of exons in a genome); transcriptome (i.e., a set of all mRNA transcripts generated in a cell or cell population, or cDNA generated by such mRNA), methylome (i.e., a population of methylation sites in a genome and a methylation pattern); microbiome (microbiome); a mixture of genomes of different organisms, a mixture of genomes of different cell types of an organism; and other complex nucleic acid mixtures (examples include, but are not limited to, microbiome, xenografts, solid tumor biopsies including both normal cells and tumor cells, etc.) comprising a large number of different nucleic acid molecules, including subsets of complex nucleic acids of the aforementioned types. In one embodiment, such complex nucleic acids have an entire sequence comprising at least one gigabase (Gb) (a diploid human genome comprises approximately 6Gb sequences).

Non-limiting examples of complex nucleic acids include "circulating nucleic acids" (CNA), which are nucleic acids that circulate in human blood or other body fluids (e.g., including but not limited to lymph, fluid, ascites, milk, urine, feces, and bronchial lavage) and can be distinguished as cell-free (CF) or cell-associated nucleic acids (reviewed in Pinzani et al., Methods 50:302-307, 2010) (e.g., circulating fetal cells in the bloodstream of an expectant mother (see, e.g., Kavanagh et al., J. Chromatol. B 878:1905-1911, 2010) or circulating tumor cells (CTCs) from the bloodstream of a cancer patient (see, e.g., Allard et al., Clin Cancer Res. 10:6897-6904, 2004)). Another example is genomic DNA from a single cell or a small number of cells, such as, for example, a small number of cells from a biopsy (e.g., fetal cells obtained from a blastocyst trophectoderm biopsy; cancer cells from a needle aspirate of a solid tumor; etc.). Another example is pathogens in tissue, blood or other body fluids, such as bacterial cells, viruses or other pathogens, etc.

As used herein, the term "target nucleic acid" (or polynucleotide) or "nucleic acid of interest" refers to any nucleic acid (or polynucleotide) suitable for processing and sequencing by the methods described herein. Nucleic acids can be single-stranded or double-stranded and can include DNA, RNA or other known nucleic acids. Target nucleic acids can be those of any organism, including but not limited to viruses, bacteria, yeast, plants, fish, reptiles, amphibians, birds and mammals (including but not limited to mice, rats, dogs, cats, goats, sheep, cattle, horses, pigs, rabbits, monkeys and other non-human primates and humans). Target nucleic acids can be obtained from an individual or multiple individuals (i.e., a population). The sample from which nucleic acid is obtained can contain nucleic acids from a mixture of cells or even organisms, such as: a human saliva sample comprising human cells and bacterial cells; a mouse xenograft comprising mouse cells and cells from a transplanted human tumor; and the like.

The target nucleic acid can be unamplified or can be amplified by any suitable nucleic acid amplification method known in the art. The target nucleic acid can be purified according to methods known in the art to remove cellular and subcellular impurities (lipids, proteins, carbohydrates, nucleic acids different from those to be sequenced, etc.), or they can be unpurified, i.e., include at least some cellular and subcellular impurities, including but not limited to intact cells that are damaged to release their nucleic acids for processing and sequencing. Target nucleic acids can be obtained from any suitable sample using methods known in the art. Such samples include but are not limited to: tissue, isolated cells or cell cultures, body fluids (including but not limited to blood, urine, serum, lymph, saliva, anal and vaginal secretions, sweat and semen); air, agriculture, water and soil samples, etc. In one aspect, nucleic acid constructs of the present invention are formed from genomic DNA.

The high coverage of shotgun sequencing is desirable because it can overcome the errors in base calls and assembly. As used in this article, for any given position in assembly sequence (assembled sequence), the term "sequence coverage surplus", "sequence coverage" or only "covering" means the quantity of the readings representing the position. It can be calculated with N x L/G from the length (G) of the initial genome, the number of readings (N) and the average reading length (L). Coverage can also be directly calculated by carrying out the counting of bases to each reference position. For full genome sequence, coverage is represented by the mean value of all bases in the assembly sequence. Sequence coverage is the average number of times (as described above) that base is read. It is often represented by "multiple coverage", for example "40 times of coverage", which means that each base represents with an average of 40 readings in the final assembly sequence.

As used herein, the term "call rate" means the percentage comparison of bases that are fully responded in a complex nucleic acid, typically with reference to a suitable reference sequence, such as, for example, a reference genome. Therefore, for the whole human genome, a "genome call rate" (or simply "call rate") is the percentage of bases that are fully responded in the human genome relative to the whole human genome reference. An "exome call rate" is the percentage of bases that are fully responded in the exon group relative to the exon group reference. The exome sequence can be obtained by sequencing a portion of the genome enriched by a plurality of known methods for selectively capturing a target genomic region from a DNA sample. Alternatively, the exome sequence can be obtained by sequencing a whole human genome comprising an exome sequence. In this way, the whole human genome sequence can have both a "genome call rate" and an "exome call rate". There is also a "raw read result call rate" that reflects the number of bases assigned A/C/G/T, rather than the total number of bases attempted. (Occasionally, the term "coverage" is used instead of "call rate," but the meaning will be apparent from the context).

Preparation of fragments of complex nucleic acids

Nucleic acid isolation. Target genomic DNA is isolated using conventional techniques, such as those disclosed in Sambrook and Russell, Molecular Cloning: A Laboratory Manual, cited above. In some cases, particularly if small amounts of DNA are employed in a particular step, it may be advantageous to provide a carrier DNA, such as an unrelated circular synthetic double-stranded DNA, to be mixed and used with the sample DNA whenever only small amounts of sample DNA are available and there is a risk of loss, e.g., through nonspecific binding to container walls.

According to some embodiments of the invention, genomic DNA or other complex nucleic acids are obtained from a single cell or a small number of cells, with or without purification.

Long fragments are desirable for LFR. Long fragments of genomic nucleic acid can be isolated from cells by many different methods. In one embodiment, the cells are lysed and the complete nuclear precipitate is obtained by a gentle centrifugation step. Then, the genomic DNA is released via proteinase K and RNase digestion for several hours. The material can be processed to reduce the concentration of residual cell waste, for example, by dialysis for a period of time (i.e., 2-16 hours) and/or dilution. Since this type of method does not require the use of many destructive methods (such as ethanol precipitation, centrifugation, and vortex oscillation), genomic nucleic acid remains intact to a large extent, producing most fragments with a length exceeding 150 kilobases. In some embodiments, the length of the fragment is about 5 to about 750 kilobases. In other embodiments, the length of the fragment is about 150 to about 600, about 200 to about 500, about 250 to about 400, and about 300 to about 350 kilobases. The minimum fragment that can be used for LFR is a fragment (about 2-5kb) containing at least two heterozygosities, and there is no maximum theoretical size, although the fragment length can be limited due to the shearing derived from the initial nucleic acid preparation operation. Techniques that produce larger fragments result in fewer aliquots being required, and those that produce shorter fragments may require more aliquots.

In case the DNA is separated and before being sampled in single wells in equal parts, it is carefully fragmented to avoid the loss of material, particularly from the sequence at each fragment end, because the loss of such material can result in the breach in the final genome assembly. In one embodiment, sequence loss is avoided by using a rare nicking enzyme that creates a polymerase, such as the start site of the phi29 polymerase, at a distance of approximately 100kb from each other. Because the polymerase creates a new DNA chain, it displaces the old chain, which creates overlapping sequences near the polymerase start site. Therefore, very few sequence deletions are arranged.

Controlled use of 5' exonucleases (before or during amplification, for example, by MDA) can promote multiple replications of initial DNA from a single cell, thus minimizing the growth of early errors through copy replication.

In other embodiments, long DNA fragments are isolated and manipulated in a manner that minimizes shearing or adsorption of DNA to the container, including, for example, separation of cells in agarose or oil in agarose gel plugs, or use of specially coated tubes and plates.

In some embodiments, further replication of fragmented DNA from a single cell before aliquoting can be achieved by ligating adapters to single-stranded priming overhangs and using adapter-specific primers and phi29 polymerase to generate two copies of each long fragment. This can generate the equivalent of four cells' worth of DNA from a single cell.

Fragmentation. The target genomic DNA is then fractionated or fragmented to the desired size by conventional techniques including enzymatic digestion, shearing, or sonication, the latter two of which are particularly useful in the present invention.

The fragment size of the target nucleic acid can vary depending on the source target nucleic acid and the library construction method used, but for standard whole genome sequencing, the length of such fragments is generally in the range of 50 to 600 nucleotides. In another embodiment, the length of the fragment is 300 to 600 or 200 to 2000 nucleotides. In yet another embodiment, the length of the fragment is 10-100, 50-100, 50-300, 100-200, 200-300, 50-400, 100-400, 200-400, 300-400, 400-500, 400-600, 500-600, 50-1000, 100-1000, 200-1000, 300-1000, 400-1000, 500-1000, 600-1000, 700-1000, 700-900, 700-800, 800-1000, 900-1000, 1500-2000, 1750-2000 and 50-2000 nucleotides. Longer fragments can be used for LFR.

In other embodiments, fragments of a specific size or within a specific size range are isolated. Such methods are well known in the art. For example, gel fractionation can be used to generate a population of fragments of a specific size within a certain base pair range, for example, 500 base pairs + 50 base pairs.

In many cases, enzymatic digestion of the extracted DNA is not required because the shearing forces generated during cracking and extraction can generate fragments in the desired range. In other embodiments, restriction endonucleases can be used to generate shorter fragments (1-5 kb) by enzymatic fragmentation. In another embodiment, about 10 to about 1,000,000 genome equivalents of DNA ensure that the fragment colony covers the entire genome. A library containing nucleic acid templates generated from such colonies of overlapping fragments can thus include a target nucleic acid, the sequence of which, once identified and assembled, can provide a major part or the entire sequence of the entire genome.

In some embodiments of the present invention, fragments are prepared using a controlled random enzymatic ("CoRE") fragmentation method. CoRE fragmentation is an enzymatic endpoint assay and has the advantages of enzymatic fragmentation (such as the ability to use it on lower amounts and/or volumes of DNA) without many of its disadvantages (including sensitivity to variations in substrate or enzyme concentrations and sensitivity to digestion time).

In one aspect, the present invention provides a fragmentation method referred to herein as controlled random enzymatic (CoRE) fragmentation, which can be used alone or in combination with other mechanical and enzymatic fragmentation methods known in the art. CoRE fragmentation involves a series of three enzymatic steps. First, the nucleic acid is subjected to an amplification method that is carried out in the presence of dNTPs doped with a certain ratio of deoxyuracil ("dU") or uracil ("U") to cause the dUTP or UTP substitution at a defined and controllable ratio of T positions in the two chains of the amplification product. Any suitable amplification method can be used in this step of the present invention. In certain embodiments, multiple displacement amplification (MDA) in the presence of dNTPs doped with dUTP or UTP at a defined ratio with dTTP is used to produce an amplification product with dUTP or UTP substituted into certain points on the two chains.

After amplification and insertion of the uracil module, the uracil is then typically removed via a combination of UDG, EndoVIII, and T4PNK to create single base gaps with functional 5' phosphate and 3' hydroxyl termini. Single base gaps are created at an average interval defined by the frequency of U in the MDA product. That is, the higher the amount of dUTP, the shorter the resulting fragments. As will be appreciated by those skilled in the art, other techniques that result in the selective replacement of nucleotides with modified nucleotides that can similarly produce cleavage, such as chemical or other enzymatically susceptible nucleotides, may also be used.

Treatment of a nicked nucleic acid with a polymerase having exonuclease activity causes the nick to "translate" or "shift" along the length of the nucleic acid until the nicks on opposite strands converge, thereby creating a double-strand break, which produces a relative population of double-stranded fragments of relatively homogeneous size. The exonuclease activity of a polymerase such as Taq polymerase cuts a short DNA strand close to the nick, while the polymerase activity "fills in" the nick and subsequent nucleotides in the strand (actually, Taq moves along the strand, using its exonuclease activity to excise bases and add the same base, with the result that the nick shifts along the strand until the enzyme reaches the end).

Because the size distribution of double-stranded fragments is a result of the ratio of dTTP to dUTP or UTP used in the MDA reaction, rather than the duration or extent of the enzymatic treatment, this CoRE fragmentation method results in a high degree of fragmentation reproducibility, generating a population of double-stranded nucleic acid fragments that are all similar sizes.

Fragment End Repair and Modification. In certain embodiments, after fragmentation, the target nucleic acids are further modified to prepare them for insertion into multiple adaptors according to the methods of the invention.

After physical fragmentation, the target nucleic acid generally has a combination of blunt ends and overhangs and a combination of phosphate radical and hydroxyl chemistry at the end. In this embodiment, the target nucleic acid is treated with several enzymes to create a blunt end with specific chemistry. In one embodiment, a polymerase and dNTPs are used to fill any 5' single strands of the overhang to create a blunt end. A polymerase with 3' exonuclease activity (generally but not always the same enzyme as the 5' active enzyme, such as T4 polymerase) is used to remove the 3' overhang. Suitable polymerases include but are not limited to T4 polymerase, Taq polymerase, Escherichia coli DNA polymerase 1, Klenow fragment, reverse transcriptase, phi29 related polymerases, including derivatives of wild-type phi29 polymerase and such polymerases, T7 DNA polymerase, T5 DNA polymerase, RNA polymerase. These techniques can be used to generate blunt ends, which can be used for a variety of applications.

In another optional embodiment, the end chemistry is altered to prevent the target nucleic acids from ligating to each other. For example, in addition to polymerases, protein kinases can also be used to create blunt ends by utilizing their 3' phosphatase activity to convert 3' phosphate groups into hydroxyl groups. Such kinases can include, but are not limited to, commercially available kinases such as T4 kinase, as well as non-commercial kinases that have the desired activity.

Similarly, phosphatases can be used to convert terminal phosphate groups into hydroxyl groups. Suitable phosphatases include, but are not limited to, alkaline phosphatase (including calf intestinal phosphatase), Antarctic phosphatase, apyrase, pyrophosphatase, inorganic (yeast) thermostable inorganic pyrophosphatase, and the like, which are known in the art.

These modifications prevent the target nucleic acids from being connected to each other in the subsequent steps of the method of the present invention, thus ensuring that during the step in which the adapter (and/or adapter arm) is connected to the end of the target nucleic acid, the target nucleic acid will be connected to the adapter and not to other target nucleic acids. The target nucleic acid can be connected to the adapter in a desired direction. The modified ends avoid unwanted structures in which the target nucleic acids are connected to each other and/or the adapters are connected to each other. The direction in which each adapter-target nucleic acid is connected can also be controlled via the terminal chemistry of the control adapter and the target nucleic acid. Such modifications can prevent the creation of nucleic acid templates containing different fragments connected with unknown structures, thus reducing and/or eliminating errors in sequence identification and assembly that can be derived from such unwanted templates.

The DNA can be denatured following fragmentation to generate single-stranded fragments.

Amplification. In one embodiment, after fragmentation (and indeed before or after any step outlined herein), an amplification step can be applied to the fragmented nucleic acid population to ensure that all fragments of sufficiently large concentrations are available for subsequent steps. According to one embodiment of the present invention, a method for sequencing a small amount of complex nucleic acids, including those of higher organisms, is provided, wherein such complex nucleic acids are amplified to generate enough nucleic acids for sequencing by the methods described herein. The sequencing methods described herein provide highly accurate sequences with a high response rate using a gene equivalent as starting material in the case of sufficient amplification. Note that cells contain approximately 6.6 picograms (pg) of genomic DNA. The full genome or other complex nucleic acids of a small amount of cells from a single cell or organism (including higher organisms such as humans) can be implemented by the methods of the present invention. Sequencing of complex nucleic acids of higher organisms can be achieved using 1 pg, 5 pg, 10 pg, 30 pg, 50 pg, 100 pg or 1 ng of complex nucleic acid as starting material, which is amplified by any nucleic acid amplification method known in the art to generate, for example, 200 ng, 400 ng, 600 ng, 800 ng, 1 μg, 2 μg, 3 μg, 4 μg, 5 μg, 10 μg or more of complex nucleic acid. We also disclose nucleic acid amplification protocols that minimize GC bias. However, the need for amplification and subsequent GC bias can be further reduced by isolating only one cell or a small number of cells, culturing them for a sufficient time under appropriate culture conditions known in the art, and using the progeny of one or more starting cells for sequencing.

Such amplification methods include, but are not limited to: multiple displacement amplification (MDA), polymerase chain reaction (PCR), ligation chain reaction (sometimes called oligonucleotide ligase amplification, OLA), cycling probe technology (CPT), strand displacement assay (SDA), transcription-mediated amplification (TMA), nucleic acid sequence-based amplification (NASBA), rolling circle amplification (RCA) (for circularized fragments), and invasive cleavage technology.

Amplification can be performed after fragmentation or before or after any of the steps outlined herein.

MDA Amplification Protocol with Reduced GC Bias In one aspect, the present invention provides a method for preparing samples wherein about 10 Mb of DNA per aliquot is faithfully amplified, for example, about 30,000-fold based on the amount of starting DNA, prior to library construction and sequencing.

According to one embodiment of the LFR method of the present invention, LFR begins by treating genomic nucleic acid, typically genomic DNA, with a 5' exonuclease to create 3' single-stranded protrusions. Such single-stranded protrusions serve as MDA start sites. The use of exonucleases also eliminates the need for a pre-amplification heat or alkaline denaturation step and does not introduce bias into the fragment population. In another embodiment, alkaline denaturation is combined with a 5' exonuclease treatment, which results in a greater reduction in bias than seen with either treatment alone. The DNA treated with the 5' exonuclease and optionally alkaline denaturation is then diluted to a subgenomic concentration and dispersed among a plurality of aliquots, as discussed above. After being divided into aliquots, for example, between a plurality of wells, the fragments in each aliquot are amplified.

In one embodiment, multiple displacement amplification (MDA) based on phi29 is used. Many studies have examined the scope of unwanted amplification bias, background product formation, and chimeric artifacts introduced via phi29-based MDA, but many of these shortcomings have occurred under extreme amplification conditions (greater than 1 million times). Typically, LFR uses substantially lower amplification levels and starts with long DNA fragments (e.g., about 100 kb), which produces effective MDA and more acceptable amplification bias levels and other amplification-related issues.

We have developed improved MDA scheme to overcome the problems relevant with the MDA using various additives (such as DNA modifying enzymes, sugars and/or chemicals, such as DMSO), and/or reduce, improve or replace the different components of MDA reaction conditions to further improve the scheme. In order to minimize the chimera, following reagent may also be included, which is used to reduce the availability of the displacement single-stranded DNA that plays an incorrect template role for the DNA chain (which is the common mechanism of chimera formation) that is extended. The main source of the coverage preference introduced by MDA is caused by the amplification difference between the GC-rich region and the AT-rich region. This can be corrected by using the different reagents in the MDA reaction and/or by regulating primer concentration to produce an environment uniformly triggered in all %GC intervals of genome. In some embodiments, random hexamer is used in triggering MDA. In other embodiments, other primer designs are utilized to reduce the preference. In other embodiments, the use of a 5' exonuclease before or during MDA can help initiate successful priming with low bias, particularly with longer (i.e., 200 kb to 1 Mb) fragments that can be used to sequence regions characterized by long segment duplications (i.e., in some cancer cells) and complex repeats.

In some embodiments, improved, more efficient fragmentation and ligation steps are used that reduce the number of MDA amplification rounds required to prepare a sample by up to 10,000-fold, which further reduces MDA-derived bias and chimera formation.

In some embodiments, the MDA reaction is designed to introduce uracil into the amplified product to prepare for CoRE fragmentation. In some embodiments, a standard MDA reaction using random hexamers is used to amplify the fragments in each well; alternatively, random 8-mer primers can be used to reduce amplification bias (e.g., GC bias) in the fragment population. In other embodiments, several different enzymes can also be added to the MDA reaction to reduce amplification bias. For example, low concentrations of non-processive 5' exonucleases and/or single-stranded binding proteins can be used to create binding sites for 8-mers. Chemical agents such as betaine, DMSO, and trehalose can also be used to reduce bias.

After amplifying the fragments in each aliquot, optionally, the amplified product can be subjected to another round of fragmentation. In some embodiments, the CoRE method is used to further fragment the fragments in each aliquot after amplification. In such embodiments, the MDA amplification of the fragments in each aliquot is designed to incorporate uracil into the MDA product. Each aliquot containing the MDA product is treated with a mixture of uracil DNA glycosylase (UDG), DNA glycosylase-lyase endonuclease VIII and T4 polynucleotide kinase to excise the uracil base and create a single base gap with a functional 5' phosphate root and 3' hydroxyl group. Double-stranded blunt-end breaks are caused by nick translation using a polymerase such as Taq polymerase, which generates ligatable fragments of a size range that depends on the dUTP concentration added in the MDA reaction. In some embodiments, the CoRE method used involves removing uracil by phi29 polymerization and chain displacement. Fragmentation of the MDA product can also be achieved via ultrasonic treatment or enzymatic treatment. Enzymatic treatments that can be used in this embodiment include, but are not limited to, DNase I, T7 endonuclease I, micrococcal nuclease, and the like.

After the MDA product is fragmented, the ends of the resulting fragments can be repaired. Many fragmentation techniques can generate ends with overhangs and ends with functional groups that cannot be used for subsequent ligation reactions, such as 3' and 5' hydroxyl groups and/or 3' and 5' phosphate groups. Having fragments that are repaired to have blunt ends can be useful. It is also desirable to modify the ends to add or remove phosphate and hydroxyl groups, thereby preventing "polymerization" of the target sequence. For example, phosphatase can be used to eliminate phosphate groups so that all ends contain hydroxyl groups. Each end can then be selectively changed to allow for connection between the desired components. One end of the "activated" fragment can then be treated with alkaline phosphatase. The fragments can then be tagged with adapters to identify fragments from the same aliquot in the LFR method.

The fragments in each aliquot are tagged. After amplification, the DNA in each aliquot is tagged, thereby identifying the aliquot from which each fragment originated. In other embodiments, the amplified DNA in each aliquot can be further fragmented before tagging with adapters so that fragments from the same aliquot all contain the same tag; see, for example, US 2007/0072208, which is incorporated herein by reference.

According to one embodiment, the adapters are designed in two sections: one section is common to all wells and the blunt ends directly connect the fragments using the methods further described herein. The "common" adapter is added as two adapter arms: one arm is a blunt end that connects to the 5' end of the fragment, and the other arm is a blunt end that connects to the 3' end of the fragment. The second section of the tagging adapter is a "barcode" section that is unique to each well. This barcode is generally a unique nucleotide sequence, and the same barcode is given to each fragment in a particular well. In this way, when the tagged fragments from all wells are recombined for sequencing applications, fragments from the same well can be identified via the identification barcode adapter. The barcode is connected to the 5' end of the common adapter arm. The common adapter and the barcode adapter can be connected to the fragments sequentially or simultaneously. As will be described in more detail herein, the ends of the common adapter and the barcode adapter can be modified so that each adapter segment will connect in the correct orientation and to the correct molecule. Such modifications prevent "polymerization" of adapter segments or fragments by ensuring that the fragments cannot ligate to each other and that adapter segments can only ligate in the exemplified orientation.

In another embodiment, a three-segment design is utilized for the adapters used to tag the fragments in each well. This embodiment is similar to the barcode adapter design described above, except that the barcode adapter segment is divided into two segments. This design allows for a large number of possible barcodes, which is achieved by allowing the combination of barcode adapter segments by connecting different barcode segments together to form a complete barcode segment generation. This combined design provides a larger complete set of possible barcode adapters while reducing the number of full-size barcode adapters that need to be generated. In another embodiment, 8-12 base pair error correction barcodes are used to achieve unique identification of each aliquot. In some embodiments, the same number of adapters as the wells are used (384 and 1536 in the above non-limiting example). In another embodiment, the cost associated with generating adapters is reduced by a new combined labeling method based on two groups of 40 half-barcode adapters.

In one embodiment, library construction involves the use of two different adapters. The A and B adapters are easily modified to contain different half-barcode sequences to produce thousands of combinations. In another embodiment, a barcode sequence is incorporated into the same adapter. This can be achieved by dividing the B adapter into two parts, each of which has a half-barcode sequence separated by a common protruding sequence for connection. The two tag components each have 4-6 bases. An 8-base (2x 4 base) tag set can uniquely label 65,000 aliquots. An additional base (2x 5 bases) allows error detection, and 12 base tags (2x 6 bases, 12 million unique barcode sequences) can be designed to allow substantial error detection and correction using the Reed-Solomon design in 10,000 or more aliquots (U.S. patent application 12/697,995, published as US 2010/0199155, which is incorporated herein by reference). Both 2 x 5 base and 2 x 6 base tags can include the use of degenerate bases (ie, "wild cards") to achieve optimal decoding efficiency.

After the fragments in each well are tagged, all the fragments are combined or pooled to form a single population. These fragments can then be used to generate nucleic acid templates or library constructs for sequencing. The nucleic acid templates generated from these tagged fragments can be identified as belonging to a specific well based on the barcode tag adapter attached to each fragment.

Long fragment read (LFR) technology

Overview

Individual human genomes are diploid in nature, with half of the homologous chromosomes originating from each parent. The context in which mutations occur on each individual chromosome can have a profound impact on the expression and regulation of genes and other transcriptional regions of the genome. In addition, determining whether two potentially harmful mutations occur within one or both alleles of a gene has extremely important clinical implications.

The present method for whole genome sequencing lacks the ability of assembling parent chromosomes separately and describing the background (haplotype) of variation co-occurrence in a cost-effective manner. Simulation experiments have shown that chromosome level haplotype determination requires allele linkage information between at least 70-100kb range. This can not be achieved by using the prior art of amplified DNA, which is limited to readings less than 1000 bases due to the loss of linkage information in the long DNA molecules and order-checking that are difficult to consistently amplify. Pairing technology can provide the equivalent of extended reading length, but is limited to less than 10kb due to the low efficiency (due to the difficulty of the circular DNA longer than several kb) of generating this type of DNA library. This method also requires extreme reading coverage to contact all heterozygotes.

If single-molecule sequencing of DNA fragments larger than 100 kb is feasible, it could be used for haplotype determination when the accuracy of single-molecule sequencing is high and the detection/instrumentation costs are low. This is very difficult to achieve with high yield for short molecules, let alone 100 kb fragments.

The recent human genome sequencing has been implemented on short read result length (<200bp), highly parallelized systems, starting with hundreds of nanograms of DNA. These technologies are remarkable in generating large amounts of data quickly and economically. Unfortunately, short read results often paired with small pairing gap sizes (500bp-10kb) eliminate most of the SNP phase information beyond several thousand bases (McKernan et al., Genome Res.19:1527,2009). In addition, it is very difficult to maintain longer DNA fragments in the absence of multiple processing steps of fragmentation due to shearing.

Currently, three human genomes, namely those of J. Craig Venter (Levy et al., PLoS Biol. 5:e254, 2007) (one of Gujarati descent in India (HapMap sample NA20847; Kitzman et al., Nat. Biotechnol. 29:59, 2011) and two of European descent (Max Planck One [MP1]; Suk et al., Genome Res., 2011; genome.cshlp.org/content/early/2011/09/02/gr.125047.111.full.pdf; and HapMap Sample NA 12878; Duitama et al., Nucl. Acids Res. 40:2041-2053, 2012)), have been sequenced and assembled as diploids. All involve cloning long DNA fragments into constructs using a method similar to bacterial artificial chromosome (BAC) sequencing used during the construction of the human reference genome (Venter et al., Science 291: 1304, 2001; Lander et al., Nature 409: 860, 2001). Although these methods generate long phased contigs (N50s of 350 kb [Levy et al., PLoS Biol. 5: e254, 2007], 386 kb [Kitzman et al., Nat. Biotechnol. 29: 59-63, 2011], and 1 Mb [Suk et al., Genome Res. 21: 1672-1685, 2011]), they require large amounts of initial DNA, extensive library processing, and are too expensive to be used in routine clinical settings.

In addition, whole chromosome haplotype determination has been demonstrated via direct isolation of metaphase chromosomes (Zhang et al., Nat. Genet. 38:382-387, 2006; Ma et al., Nat. Methods 7:299-301, 2010; Fan et al., Nat. Biotechnol. 29:51-57, 2011; Yang et al., Proc. Natl. Acad. Sci. USA 108:12-17, 2011). These methods are excellent for long-range haplotype determination, but have not yet been used for whole genome sequencing and require preparation and isolation of whole metaphase chromosomes, which can be challenging for some clinical samples.

The LFR method overcomes these limitations. LFR includes DNA preparation and tagging together with associated algorithms and software to enable accurate assembly of separate sequences of parental chromosomes (i.e., complete haplotype determination) in diploid genomes at significantly reduced experimental and computational costs.

LFR is based on the physical separation of long fragments of genomic DNA (or other nucleic acids) between multiple different aliquots, so that there is a low probability of any given region of the genome of both the maternal and paternal components present in the same aliquots. By placing a unique identifier in each aliquot and summing up multiple aliquots for analysis, DNA sequence data can be assembled into a diploid genome, for example, the sequence of each parental chromosome can be measured. LFR does not require the fragment of complex nucleic acid to be cloned into a vector, as in the unit type determination method using large fragments (such as BAC) libraries. LFR does not need to directly separate the individual chromosomes of an organism. Finally, LFR can be implemented to individual organisms and does not require an organism colony to realize unit type phasing.

As used herein, the term "vector" refers to a plasmid or viral vector into which a foreign DNA segment is inserted. A vector is used to introduce foreign DNA into a suitable host cell, where the vector and the inserted foreign DNA replicate due to the presence of, for example, a functional origin of replication or an autonomously replicating sequence in the vector. As used herein, the term "clone" refers to the insertion of a DNA segment into a vector and the replication of the vector with the inserted foreign DNA in a suitable host cell.

LFR can be used with the sequencing methods discussed in detail herein, and more generally as a pre-processing method with any sequencing technology known in the art, including both short read and longer read methods. LFR can also be used in conjunction with various types of analysis, including, for example, analysis of transcriptomes, methylomes, and the like. Because it requires very little input DNA, LFR can be used to sequence one or a small number of cells and determine unit types, which can be particularly important for cancer, prenatal diagnostics, and personalized medicine. This can facilitate the identification of familial genetic diseases, among others. By making it possible to distinguish responses from two sets of chromosomes in a diploid sample, LFR also allows for higher confidence responses of low-coverage variants and non-variant positions. Other applications of LFR include parsing full-length sequencing of extensive rearrangements and alternatively spliced transcripts in cancer genomes.

LFR can be used to process and analyze complex nucleic acids, including but not limited to genomic DNA, purified or unpurified, including cells and tissues that are gently disrupted to release such complex nucleic acids without shearing and excessive fragmentation of such complex nucleic acids.

In one aspect, LFR generates a virtual read length of about 100-1000 kb in length.

In addition, LFR can also significantly reduce the computational requirements and associated costs of any short read result technology. Importantly, LFR eliminates the need to extend the length of the read result sequencing result (if it reduces the overall yield). Another benefit of LFR is that the errors or questionable base calls that can be derived from current sequencing technology are substantially (10 to 1000 times) reduced, typically 1 per 100kb, or 30,000 false positive calls per human chromosome genome, and a similar number of undetected variants per human genome. This significant reduction in error minimizes the need for tracking the construction of detection variants and facilitates the use of human genome sequencing for diagnostic applications.

In addition to being applicable to all sequencing platforms, LFR-based sequencing can be applied to any application, including but not limited to the study of structural rearrangements in cancer genomes, global methylome analysis, including haplotypes of methylation sites, and even de novo assembly applications for metagenomics or new genome sequencing of complex polyploid genomes, such as those found in plants.

In contrast to the consensus sequence of only parental or related chromosomes, LFR provides the ability to obtain the true sequence of each chromosome (despite its high similarity and the presence of long repeats and segmental duplications). To generate such data, sequence continuity is generally established over long DNA ranges, such as 100 kb to 1 Mb.

Yet another aspect of the invention includes software and algorithms for efficiently utilizing LFR data for whole chromosome haplotype and structural variation mapping and false positive/negative error correction to fewer than 300 errors per human chromosome.

In yet another aspect, the LFR technology of the present invention reduces the complexity of the DNA in each aliquot by 100-1000 fold, depending on the aliquot and the number of cells used. Complexity reduction and unit type separation in long DNAs greater than 100 kb can help to more efficiently and cost-effectively (up to 100-fold reduction in cost) assemble and detect all variations in human and other diploid genomes.

The LFR method described herein can be used as a pre-treatment step for sequencing diploid genomes using any sequencing method known in the art. In other embodiments, the LFR method described herein can be used on many sequencing platforms, including, for example, but not limited to, polymerase-based synthetic sequencing (e.g., HiSeq 2500 system, Illumina, San Diego, CA), ligation-based sequencing (e.g., SOLiD 5500, Life Technologies Corporation, Carlsbad, CA), ion semiconductor sequencing (e.g., ion PGM or ion proton sequencer, Life Technologies Corporation, Carlsbad, CA), zero-mode waveguide (e.g., PacBio RS sequencer, Pacific Biosciences, Menlo Park, CA), nanopore sequencing (e.g., Oxford Nanopore Technologies Ltd., Oxford, United Kingdom), pyrophosphate sequencing (e.g., 454 Life Sciences, Branford, CT), or other sequencing technologies. Some of these sequencing technologies are short read technologies, but other technologies produce longer reads, such as GS FLX+ (454 Life Sciences; up to 1000 bp), PacBio RS (Pacific Biosciences; approximately 1000 bp), and nanopore sequencing (Oxford Nanopore Technologies Ltd.; 100 kb). For unit type phasing, longer reads are advantageous, requiring much less computation, although they tend to have higher error rates, and errors in such long reads may need to be identified and corrected according to the methods outlined herein before unit type phasing.

According to one embodiment of the present invention, the basic steps of LFR include: (1) dividing long fragments of complex nucleic acid (e.g., genomic DNA) into aliquots, each aliquot containing one genome equivalent of DNA; (2) amplifying the genomic fragments in each aliquot; (3) fragmenting the amplified genomic fragments to create short fragments of a size suitable for library construction (e.g., about 500 bases in length in one embodiment); (4) tagging the short fragments to allow identification of the aliquot from which the short fragment originated; (5) pooling the tagged fragments; (6) sequencing the pooled, tagged fragments; and (7) analyzing the resulting sequence data to locate and assemble the data and obtain haplotype information. According to one embodiment, LFR uses 384-well plates with 10-20% haploid genomes in each well, resulting in a theoretical 19-38x physical coverage of both the maternal and paternal alleles for each fragment. The 19-38x redundancy of the initial DNA ensures complete genome coverage and high variant calling and phasing accuracy. LFR avoids the need for subcloning of complex nucleic acid fragments into vectors or for isolation of individual chromosomes (eg, metaphase chromosomes), and it can be fully automated, making it suitable for high-throughput, cost-effective applications.

We have also developed techniques using LFR for error reduction and other purposes detailed herein. LFR methods have been disclosed in U.S. Patent Application Nos. 12/816,365, 12/329,365, 12/266,385, and 12/265,593, and U.S. Patent Nos. 7,906,285, 7,901,891, and 7,709,197, all of which are hereby incorporated by reference in their entirety.

As used in this article, the term "haplotype" means the allele combination transmitted together at the adjacent position (locus) on the chromosome, or alternatively, a group of sequence variants statistically associated on the single chromosome of the chromosome pair. Each individual has two sets of chromosomes, i.e., a father and another mother. Usually, DNA sequencing only produces genotype information, i.e., the sequence of the disordered alleles along the DNA segment. Haplotype is inferred to be genotype and the alleles in each disordered pair are divided into two different sequences each called haplotype. Haplotype information is necessary for many different types of genetic analysis (including disease association studies and inferring colony ancestors).

As used herein, the term "phasing" (or resolution) means that sequence data are classified into two groups of parental chromosomes or haplotypes. Haplotype phasing refers to the problem of accepting a group of genotypes of an individual or a population (i.e., more than one individual) as input and outputting a pair of haplotypes (one is paternal and the other is maternal) for each individual. Phasing can involve the sequence data of a region of the genome, or as little as two sequence variants in a read result or overlapping group, which can be referred to as local phasing or microphasing. It can also involve the phasing of larger overlapping groups (generally including about 10 or more sequence variants) or even whole genome sequences, which can be referred to as "universal phasing." Optionally, sequence variants are phased during genome assembly.

Aliquot and sample multiple genome equivalents of complex nucleic acids

The LFR method is based on randomly physically dividing the genome in long fragments into multiple aliquots, so that each aliquot contains one haploid genome. As the fraction of the genome in each pool decreases, the statistical probability of having corresponding fragments from two parental chromosomes in the same pool decreases significantly.

In some embodiments, 10% of the genome equivalents are sampled in equal portions into each well of a multiwell plate. In other embodiments, 1% to 50% of the genome equivalents of the complex nucleic acid are sampled in equal portions into each well. As noted above, the number of aliquots and genome equivalents can depend on the number of aliquots, initial fragment size, or other factors. Optionally, double-stranded nucleic acids (e.g., human genomes) are denatured before aliquoting; in this way, single-stranded complements can be distributed into different aliquots. According to one embodiment, each aliquot contains 2, 4, 6, or more copies (or complements) of the majority of the chains of the complex nucleic acid (or 2, 4, 6, or more complements if double-stranded nucleic acids are denatured before aliquoting).

For example, at 0.1 genome equivalents per aliquot (approximately 0.66 picograms or pg of DNA at approximately 6.6 pg per human genome), there is a 10% probability that two fragments will overlap, and a 50% probability that those fragments will originate from different parental chromosomes; this yields a 95% overall probability that the base pairs in the aliquots are non-overlapping, i.e., a particular aliquot will be uninformative for a given fragment because the aliquot contains fragments originating from both maternal and paternal chromosomes. Uninformative aliquots can be identified because the sequence data derived from such aliquots contain an increased amount of "noise," that is, impurities in the matrix of heterozygous pair connections. The fuzzy interference system (FIS) allows for robustness against some degree of impurities, i.e., it can perform correct ligations despite the presence of impurities (up to a certain degree). Even smaller amounts of genomic DNA can be used, particularly in the context of microdroplets or nanodroplets or emulsions, where each droplet can contain one DNA fragment (e.g., a single 50 kb fragment of genomic DNA or approximately 1.5 x 10-5 genome equivalents). Even at 50% genome equivalents, most aliquots will be informative. At higher levels, such as 70% genome equivalents, informative wells can be identified and used. According to one aspect of the invention, 0.000015, 0.0001, 0.001, 0.01, 0.1, 1, 5, 10, 15, 20, 25, 40, 50, 60, or 70% genome equivalents of complex nucleic acid are present in each aliquot.

It will be appreciated that the dilution factor may depend on the initial size of the fragment. That is, using a gentle technique to separate genomic DNA, fragments of approximately 100 kb may be obtained, which are then sampled in equal portions. Techniques that tolerate larger fragments result in fewer aliquots, while techniques that generate shorter fragments may require more dilution.

We have successfully performed all six enzymatic steps in the same reaction without DNA purification, which facilitates miniaturization and automation and makes it feasible to adapt LFR to a wide variety of platforms and sample preparation methods.

According to one embodiment, each aliquot is contained in a separate well of a multi-well plate (e.g., a 384-well plate). However, any suitable type of container or system known in the art can be used to hold the aliquots, or the LFR method can be implemented using droplets or emulsions, as described herein. According to one embodiment of the present invention, the volume is reduced to sub-microliter levels. In one embodiment, automated pipetting methods can be used in a 1536-well format.

In general, as the number of aliquots increases, for example, to 1536, and the percentage of the genome decreases to about 1% haploid genomes, the statistical support for the haplotype increases significantly because the occasional presence of both maternal and paternal haplotypes in the same well decreases. Thus, a large number of small aliquots with negligible mixed haplotype frequencies per aliquot allows for the use of fewer cells. Similarly, longer fragments (e.g., 300 kb or longer) help bridge segments lacking heterozygous loci.

Nanoliter (nl) dispensing tools (e.g., Hamilton Robotics Nano pipetting heads, TTP LabTech Mosquito, etc.) that provide 50-100nl contactless pipetting can be used for rapid and low-cost pipetting to generate dozens of genomic libraries in parallel. The increased number of aliquots (compared to 384-well plates) results in a significant reduction in genomic complexity per well, which reduces overall computational costs by more than 10-fold and improves data quality. In addition, the automation of this method increases throughput and reduces the hands-on cost of generating libraries.

LFR using smaller aliquot volumes (including microdroplets and emulsions)

Even further cost reductions and other advantages can be achieved using microdroplets. In some embodiments, LFR is performed with combinatorial tagging in emulsions or microfluidic devices. Volumes down to picoliter levels in 10,000 aliquots can achieve even greater cost reductions due to lower reagent and computational costs.

In one embodiment, LFR uses 10 microliters (μl) volumes of reagent for each well in a 384-well format. For example, such volume can be reduced to by using a commercial automated pipetting method in a 1536-well format. Further volume reduction can be achieved using a nanoliter (nl) dispensing tool (such as Hamilton Robotics Nano pipetting head, TTP LabTech Mosquito, etc.) providing 50-100nl contactless pipetting, which can be used for fast and low-cost pipetting to produce tens of genomic libraries in parallel. Increasing the number of aliquots results in a larger reduction in genomic complexity in each well, which reduces overall computational cost and improves data quality. In addition, the automation of this method improves throughput and reduces the cost of producing a library.

In other embodiments, an 8-12 base pair error correction barcode is used to achieve unique identification of each aliquot.In some embodiments, the same number of adapters as wells are used.

In another embodiment, a novel combinatorial tagging method is used, which is based on two groups of 40 half-barcode adapters. In one embodiment, library construction involves the use of two different adapters. The A and B adapters are easily modified to each contain a different half-barcode sequence to produce thousands of combinations. In another embodiment, a barcode sequence is incorporated into the same adapter. This can be achieved by dividing the B adapter into two parts, each of which has a half-barcode sequence separated by a common protruding sequence for connection. The two tag components each have 4-6 bases. An 8-base (2x 4 base) tag group can uniquely tag 65,000 aliquots. An additional base (2x 5 bases) allows error detection, and 12 base tags (2x 6 bases, 12 million unique barcode sequences) can be designed to allow substantial error detection and correction using the Reed-Solomon design in 10,000 or more aliquots. In exemplary embodiments, both 2 x 5-base and 2 x 6-base tags are employed, including the use of degenerate bases (ie, "wild-cards") to achieve optimal decoding efficiency.

Volumes down to the picoliter level (e.g., in 10,000 aliquots) can achieve even greater reductions in reagent and computational costs. In some embodiments, this level of cost reduction and extensive aliquot sampling is achieved by combining the LFR method with combinatorial labeling into an emulsion or microfluidic device. The ability to perform all enzymatic steps in the same reaction without DNA purification facilitates the ability to miniaturize and automate this method and results in adaptability to a wide variety of platforms and sample preparation methods.

In one embodiment, the LFR method is used in conjunction with an emulsion-based device. The first step in adapting LFR to an emulsion-based device is to prepare an emulsion reagent with barcode-tagged combinatorial adapters that have a single unique barcode per droplet. Two sets of 100 half-barcodes are sufficient to uniquely identify 10,000 aliquots. However, increasing the number of half-barcode adapters to over 300 allows for random addition of barcoded droplets to be combined with sample DNA with a low probability that any two aliquots contain the same barcode combination. Combinatorial barcode adapter droplets can be generated and stored as reagents in a single tube for thousands of LFR libraries.

In one embodiment, the present invention is scaled up from 10,000 to 100,000 or more aliquot libraries. In other embodiments, the LFR method is adapted for such scale-up by increasing the number of initial semi-barcoded adapters. These combined adapter droplets are then fused one-to-one with droplets containing ligation-ready DNA representing less than 1% of a haploid genome. Using conservative estimates of 1 nl per droplet and 10,000 droplets, this represents a total volume of 10 μl for the entire LFR library.

Recent studies have also suggested improved GC bias and reduced background amplification after amplification (eg, by MDA) obtained by reducing reaction volumes to nanoliter sizes.

There are currently several types of microfluidic devices (such as those sold by Advanced Liquid Logic, Morrisville, NC) or pico/nanodroplets (such as RainDance Technologies, Lexington, MA) with pico/nanodroplet generation, fusion (3000/ seconds) and collection capabilities, and can be used in such embodiments of LFR. In other embodiments, using improved nanopipette or acoustic droplet ejection technology (such as LabCyte Inc., Sunnyvale, CA) or using a microfluidic device capable of handling up to 9216 single reaction wells (such as those produced by Fluidigm, South San Francisco, CA), about 10-20 nanoliters are dropped in a 3072-6144 or more form (still a cost-effective total MDA volume of 60 μl, without losing the ability to calculate cost savings or to sequence the genomic DNA from a small amount of cells) in a plate or on a glass slide. Increasing the number of aliquots results in a greater reduction in genomic complexity within each well, which reduces overall computational cost and improves data quality. Additionally, automation of this method increases throughput and reduces the cost of generating libraries.

Amplification

According to one embodiment, the LFR method begins with a short treatment of genomic DNA with a 5' exonuclease to create a 3' single-stranded overhang that serves as an MDA start site. The use of an exonuclease eliminates the need for a pre-amplification heat or alkaline denaturation step and does not introduce bias into the fragment population. Alkaline denaturation can be combined with a 5' exonuclease treatment, which results in a further reduction in bias. The DNA is then diluted to a subgenomic concentration and aliquoted. After aliquoting, the fragments in each well are amplified, for example, using the MDA method. In certain embodiments, the MDA reaction is a modified phi29 polymerase-based amplification reaction, although another known amplification method can be used.

In some embodiments, the MDA reaction is designed to introduce uracil into the amplified product. In some embodiments, the standard MDA reaction utilizing random hexamer is used to amplify the fragment in every hole. In many embodiments, different from random hexamer, random 8 aggressiveness primers are used to reduce the amplification preference in the fragment colony. In other embodiments, several different enzymes can also be added to the MDA reaction to reduce the amplification preference. For example, the non-carrying 5 ' exonuclease and/or single-stranded binding protein of low concentration can be used to create a binding site for 8 aggressiveness. Chemical agents such as betaine, DMSO and trehalose can also be used to reduce the preference via similar mechanisms.

Fragmentation

According to one embodiment, after the DNA in each well is amplified, the amplified product is subjected to a round of fragmentation. In some embodiments, the fragments in each well are further fragmented using the CoRE method described above after amplification. To use the CoRE method, the MDA reaction used to amplify the fragments in each well is designed to incorporate uracil into the MDA product. Fragmentation of the MDA product can also be achieved via ultrasonication or enzymatic treatment.

If the CoRE method is used to fragment the MDA product, each well containing amplified DNA is treated with a mixture of uracil DNA glycosylase (UDG), DNA glycosylase-lyase endonuclease VIII, and T4 polynucleotide kinase to excise the uracil base and create a single base gap with a functional 5' phosphate and 3' hydroxyl group. Nick translation using a polymerase such as Taq polymerase results in double-stranded blunt-end breaks, which generate ligatable fragments of a size range that depends on the concentration of dUTP added in the MDA reaction. In some embodiments, the CoRE method used involves the removal of uracil by phi29 polymerization and strand displacement.

After fragmentation of the MDA product, the ends of the resulting fragments can be repaired. Such repair may be necessary because many fragmentation techniques can generate ends with overhangs and ends with functional groups that cannot be used for subsequent ligation reactions, such as 3' and 5' hydroxyl groups and/or 3' and 5' phosphate groups. In many aspects of the present invention, it may be useful to have fragments that are repaired to have blunt ends, and in some cases, it may be desirable to change the end photochemistry so that there is no correct phosphate and hydroxyl group orientation, thereby preventing "polymerization" of the target sequence. Control of the end chemistry can be provided using methods known in the art. For example, in some cases, the use of phosphatase eliminates all phosphate groups so that all ends contain hydroxyl groups. Then, one end of the fragment can be "activated" by alkaline phosphatase treatment. Then, each end can be selectively changed to allow connection between the desired components. Then, one end of the fragment can be "activated", in some embodiments by treatment with alkaline phosphatase.

Following fragmentation and optional end repair, the fragments are tagged with adapters.

Tagging

Typically, tag adapter arms are designed in two segments: one segment is common to all wells and blunt-ended to directly ligate fragments using methods described further herein. The second segment is unique to each well and contains a "barcode" sequence so that when the contents of each well are combined, the fragments from each well can be identified.

According to one embodiment, the "common" adapter is added as two adapter arms: one arm is a blunt end that connects to the 5' end of the fragment, and the other arm is a blunt end that connects to the 3' end of the fragment. The second section of the tagging adapter is a "barcode" segment that is unique to each well. This barcode is generally a unique nucleotide sequence, and the same barcode is given to each fragment in a particular well. In this way, when the tagged fragments from all wells are recombined for sequencing applications, fragments from the same well can be identified by identifying the barcode adapter. The barcode is connected to the 5' end of the common adapter arm. The common adapter and the barcode adapter can be connected to the fragments sequentially or simultaneously. The ends of the common adapter and the barcode adapter can be modified so that each adapter segment will connect in the correct orientation and to the correct molecule. Such modifications prevent "polymerization" of adapter segments or fragments by ensuring that fragments cannot connect to each other and that adapter segments can only connect in the illustrated orientation.

In another embodiment, a three-segment design is utilized for the adapters used to tag the fragments in each well. This embodiment is similar to the barcode adapter design described above, except that the barcode adapter segment is divided into two segments. This design allows for a large number of possible barcodes by allowing combinatorial barcode adapter segments to be generated by connecting different barcode segments together to form a complete barcode segment. This combinatorial design provides a larger repertoire of possible barcode adapters while reducing the number of full-size barcode adapters that need to be generated.

According to one embodiment, after the fragments in each well are tagged, all fragments are combined to form a single population. These fragments can then be used to generate nucleic acid templates of the present invention for sequencing. The nucleic acid templates generated from these tagged fragments can be identified as originating from a specific well based on the barcode tag adapter attached to each fragment. Similarly, after sequencing the tag, the genomic sequence attached to it can also be identified as originating from that well.

In some embodiments, the LFR methods described herein do not include multiple levels or layers of fragmentation/aliquoting, as described in U.S. Patent Application No. 11/451,692, filed June 13, 2006, which is incorporated herein by reference in its entirety for all purposes. That is, some embodiments utilize only one round of aliquoting and also allow for the re-combining of aliquots for a single array, rather than using a different array for each aliquot.

LFR using one or a small number of cells as a source of complex nucleic acids

According to one embodiment, the LFR method is used to analyze the genome of a single cell or a small number of cells. The method used to isolate DNA in this case is similar to the method described above, but can take place in a smaller volume.

As discussed above, it is possible to achieve long fragments from cell separation of genomic nucleic acid by a variety of different methods. In one embodiment, by cell lysis, and with gentle centrifugation steps by complete nuclear precipitation. Then, via proteinase K and RNAse digestion for several hours, genomic DNA is released. In some embodiments, material can be processed to reduce the concentration of residual cell waste, and this type of processing is well known in the art and can include but is not limited to dialysis for a period of time (i.e. 2-16 hours) and/or dilution. Because this type of method for separating nucleic acid does not relate to many destructive methods (such as ethanol precipitation, centrifugation and vortex oscillation), genomic nucleic acid remains intact to a great extent, produces most fragments with a length exceeding 150 kilobases. In some embodiments, the length of the fragment is about 100 to about 750 kilobases. In other embodiments, the length of the fragment is about 150 to about 600, about 200 to about 500, about 250 to about 400 and about 300 to about 350 kilobases.

Once the DNA is isolated and before it is aliquoted into individual wells, the genomic DNA must be carefully fragmented to avoid loss of material, particularly loss of sequence from the ends of each fragment, as loss of such material can result in gaps in the final genome assembly. In one instance, sequence loss was avoided by using a rare nicking enzyme that creates start sites for polymerases, such as phi29 polymerase, at a distance of approximately 100 kb from each other. As the polymerase creates a new DNA chain, it displaces the old chain, with the end result being overlapping sequences near the polymerase start site, resulting in very few sequence deletions.

In some embodiments, the controlled use of a 5' exonuclease (before or during the MDA reaction) can promote multiple replications of the initial DNA from a single cell, thereby minimizing the growth of early errors through copy replication.

In one aspect, method of the present invention produces the quality genomic data from unicellular.Assuming that there is no DNA loss, there is the benefit of using the equivalent amount DNA starting from a large amount of preparation instead of a small amount of cells (10 or less).Start with less than 10 cells and ensure the consistent coverage in the long fragment of any given region of genome to basically all DNA accurate equal-division sampling.Start with 5 or less cells to allow 4 times or larger coverage of every 100kb DNA fragment in each aliquot and do not make the total reading result number increase to be higher than 120Gb (20 times of coverage of 6Gb diploid genome). However, a large amount of aliquots (10,000 or more) and longer DNA fragments (>200kb) are even more important for order-checking from a few cells, because for any given sequence, only have the same number of starting cells as overlapping fragments, and the appearance of overlapping fragments from two parent chromosomes in an aliquot can be a devastating information loss.

LFR is perfectly suited to this problem because it produces excellent results starting with only about 10 cells, equivalent to the starting input genomic DNA, and even a single cell will provide enough DNA to implement LFR. Generally, the first step in LFR is low-bias whole genome amplification, which can be particularly useful for single-cell genomic analysis. Due to DNA chain breaks and DNA loss during processing, even single-molecule sequencing methods may require a certain level of DNA amplification from a single cell. The difficulty in sequencing single cells comes from trying to amplify the entire genome. Studies using MDA on bacteria have suffered from the loss of roughly half of the genome in the final assembled sequence and a considerable amount of variation in the coverage of those sequencing intervals. This can be explained in part by the initial genomic DNA having nicks and chain breaks, which cannot be replicated at the ends and are thus lost during the MDA method. LFR provides a solution to this problem by creating long overlapping fragments of the genome before MDA. According to one embodiment of the present invention, in order to achieve this, a gentle method is used to separate genomic DNA from cells. The largely intact genomic DNA is then gently treated with a common nicking enzyme to generate a genome with semi-random nicks. The strand displacement ability of phi29 is then used to polymerize from the nicks, creating very long (>200 kb) overlapping fragments. These fragments are then used as starting templates for LFR.

Methylation analysis using LFR

In yet another aspect, the method and composition of the present invention are used for genome methylation analysis. Several methods are currently available for global genome methylation analysis. One method involves bisulfate treatment of genomic DNA and sequencing of the genome portion obtained by fragmentation of repetitive elements or by methylation-specific restriction enzymes. This technology produces information about overall methylation, but does not provide locus-specific data. The next higher resolution level uses DNA arrays and is limited to the number of features on the chip. Finally, the highest resolution and most expensive method requires bisulfate treatment, followed by sequencing of the entire genome. Using LFR, it is possible to sequence all bases of the genome, and assemble a complete diploid genome with digital information about the methylation level (i.e., 5-base sequencing) of each cytosine position in the human genome. In addition, LFR allows 100kb or larger methylation sequence blocks to be connected to cell type sequencing, providing methylation cell type determination, which is information that cannot be achieved with any currently available method.

In one non-limiting exemplary embodiment, methylation status is obtained in a method in which genomic DNA is first aliquoted and denatured for MDA. Next, the DNA is treated with bisulfite (i.e., the step requiring denatured DNA). The remaining preparation follows the methods described in, for example, U.S. application serial numbers 11/451,692 filed June 13, 2006 and 12/335,168 filed December 15, 2008, each of which is hereby incorporated by reference in its entirety for all purposes and in particular for all teachings relating to nucleic acid analysis of fragment mixtures according to long fragment read technology.

In one aspect, MDA amplifies each strand of a specific fragment, which independently generates 50% reads for any given cytosine position that are unaffected by bisulfite (i.e., the base opposite cytosine, guanine, is unaffected by bisulfite) and 50% that provide methylation status. The reduced DNA complexity of each aliquot facilitates precise positioning and assembly of less informative, typically 3-base (A, T, G) reads.

Bisulfite treatment has been reported to fragment DNA. However, denaturation and careful titration of the bisulfite buffer can avoid extensive fragmentation of genomic DNA. A 50% conversion of cytosine to uracil can be tolerated in LFR, which reduces exposure of DNA to bisulfite and minimizes fragmentation. In some embodiments, some degree of fragmentation after aliquoting is acceptable as it does not affect haplotype determination.

Using LFR to analyze cancer genomes

It has been proposed that more than 90% of cancers contain significant losses or gains in human genomic regions, known as aneuploidy, and some individual cancers have been observed to contain more than 4 copies of some chromosomes. This increased complexity of the copy number of chromosomes and intrachromosomal regions makes sequencing cancer genomes substantially more difficult. The ability of LFR technology to sequence and assemble very long (>100 kb) genomic fragments makes it perfectly suited for sequencing complete cancer genomes.

Error reduction by sequencing target nucleic acid in multiple aliquots

According to one embodiment, even if LFR-based phasing is not implemented and standard sequencing methods are used, the target nucleic acid is divided into multiple aliquots, each containing a certain amount of target nucleic acid. In each aliquot, the target nucleic acid is fragmented (if fragmentation is required) and the fragments are tagged with aliquot-specific tags (or aliquot-specific tag sets) before amplification. Alternatively, when processing tissue samples, one or more cells can be distributed to each of multiple aliquots, followed by cell disruption, fragmentation, tagged with aliquot-specific tags, and amplification. In either case, the DNA amplified from each aliquot can be sequenced separately or combined and sequenced after combining. One advantage of this method is that errors introduced due to amplification (or other steps occurring in each aliquot) can be identified and corrected. For example, a base call (e.g., identifying a specific base, such as A, C, G, or T) at a specific position in the sequence data (e.g., relative to a reference) can be accepted as true if the base call is present in sequence data from two or more aliquots (or other threshold number of aliquots), or in a substantial majority of expected aliquots (e.g., in at least 51, 70, or 80%), where the denominator can be limited to aliquots having a base call at the specific position. A base call can include an allele that changes heterozygosity or potential heterozygosity. A base call at a specific position can be accepted as false if it is present in only one aliquot (or other threshold number of aliquots), or in a substantial minority of aliquots (e.g., less than 10, 5, or 3 aliquots, or as measured by relative numbers, such as 20 or 10%). The threshold value can be predetermined or dynamically determined based on the sequencing data. If a base call at a particular position is not present in a substantial minority and in a substantial majority of expected aliquots (e.g., in 40-60%), it can be converted/accepted as a "no call." In some embodiments and implementations, multiple parameters (e.g., in distributions, probabilities, and/or other functions or statistics) can be used to characterize what can be considered a substantial minority or a substantial majority of aliquots. Examples of such parameters include, but are not limited to, one or more of the following: the number of base calls that identify a particular base; the coverage or total number of calling bases at a particular position; the number and/or identity of unique aliquots that generate sequence data including a particular base call; the total number of unique aliquots that generate sequence data containing at least one base call at a particular position; the reference base at a particular position; and the like. In one embodiment, a combination of the above parameters for a particular base call can be input into a function to determine a score (e.g., probability) for the particular base call. The score can then be compared to one or more threshold values as part of determining whether the base call is accepted (e.g., above a threshold), erroneous (e.g., below a threshold), or a no call (e.g., if all scores for the base calls are below a threshold). The determination of a base call can depend on the scores of other base calls.

As a basic example, if base call A is present in more than 35% (example of a score) of the aliquots containing reads at a position of interest, and base call C is present in more than 35% of these aliquots, and the other base calls each have a score of less than 20%, then the position can be considered heterozygous for A and C, possibly subject to other criteria (such as a minimum number of aliquots containing reads at the position of interest). Thus, each score can be input into another function (such as a heuristic, which may use comparisons or fuzzy logic) to provide a final determination of the base call for that position.

As another example, a specific number of aliquots containing base calls can be used as a threshold. For example, when analyzing cancer samples, there may be a low prevalence of somatic mutations. In such cases, a base call may appear in less than 10% of the aliquots covering the position, but the base call may still be considered correct, possibly subject to other criteria. Thus, various embodiments may use absolute or relative numbers, or both (e.g., as input to a comparison or fuzzy logic). Furthermore, such numbers of aliquots can be input to a function (as mentioned above), along with a threshold corresponding to each number, and the function can provide a score that can also be compared to one or more thresholds to make a final determination regarding the base call at a particular position.

Other examples of error correction functions relate to sequence errors in the original reads, which result in putative variant responses that are inconsistent with other variant responses and their haplotypes. If 20 reads of variant A are present in 9 and 8 aliquots belonging to the corresponding haplotypes, and 7 reads of variant G are present in 6 wells (5 or 6 of which are shared with aliquots with A reads), then logic can reject variant G as a sequencing error because for a diploid genome, only one variant can reside in one position in each haplotype. Variant A is substantially more supported by readings, while G reads substantially follow the aliquots of A reads, indicating that they are most likely to be generated due to erroneous reads G rather than A. If G reads are almost exclusively in aliquots separated from A, then this can indicate that G reads are mislocalized or that they are from contaminating DNA.

Identification of expansions in regions with short tandem repeats

Short tandem repeats (STRs) in DNA are segments of DNA with a strongly periodic pattern. STRs occur when a pattern of two or more nucleotides repeats, and the repeated sequences are directly adjacent to each other; the repetition can be complete or incomplete, meaning there can be several base pairs that do not match the periodic motif. Generally, the pattern ranges from 2 to 5 base pairs (bp) in length. STRs are often located in noncoding regions, such as introns. Short tandem repeat polymorphisms (STRPs) occur when homologous STR loci differ in the number of repeats between individuals. STR analysis is often used to determine genetic profiles for forensic purposes. STRs present in gene exons can represent hypermutation regions associated with human disease (Madsen et al, BMC Genomics 9:410, 2008).

In the human genome (and the genomes of other organisms), STRs include trinucleotide repeats, such as CTG or CAG repeats. Trinucleotide repeat expansions, also known as triplet repeat expansions, are caused by slippage during DNA replication and are associated with certain diseases classified as trinucleotide repeat disorders, such as Huntington's disease. Generally, the larger the expansion, the more likely it is to cause disease or increase the severity of the disease. This characteristic leads to the "early onset" feature seen in trinucleotide repeat disorders, that is, a trend of decreasing age of disease onset and increasing symptom severity over successive generations of affected families due to the expansion of these repeats. Identifying trinucleotide repeat expansions can be used to accurately predict the age of onset and disease progression for trinucleotide repeat disorders.

Using next-generation sequencing methods, expansions of STRs, such as trinucleotide repeats, can be difficult to identify. Such expansions cannot be located and may be absent or underrepresented in the library. Using LFR, it is possible to see a significant decrease in sequence coverage in STR regions. For example, regions with STRs will characteristically have lower coverage levels than regions without such repeats, and if expansions are present in such regions, there will be a substantial decrease in coverage in such regions, which can be observed in a plot of coverage versus position in the genome.

Figure 14 shows the example of the detection of CTG repeat expansion in the affected embryo.LFR is used to measure the parental haplotype of the embryo.In the figure of the clone coverage to position of the mean standardization, the haplotype with expanded CTG repeat does not have or has a very small amount of DNB that passes through the expansion zone, resulting in the reduction of coverage in the region.Reduction can also be detected in the combined sequence coverage of the two haplotypes; However, the decline of a haplotype may be more difficult to identify.For example, if the sequence coverage is an average of about 20, the region with the expansion zone will have a significant decline, for example, if the affected haplotype has 0 coverage in the expansion zone, then it will decline to 10.Like this, a 50% decline will occur.However, if the sequence coverage of the two haplotypes is compared, then coverage is 10 in the normal haplotype, and 0 in the affected haplotype, which is a decline of 10, but the overall percentage decline is 100%. Alternatively, one can analyze the relative amounts, which for the combined sequence coverage is 2:1 (normal versus coverage in the expanded region), but is 10:0 (haplotype 1 versus haplotype 2), which is either infinity or 0 (depending on how the ratios are formed), so a larger difference.

Diagnostic uses of sequence data

The sequence data generated using the methods of the present invention can be used for a wide variety of purposes. According to one embodiment, the sequencing methods of the present invention are used to identify sequence variations in complex nucleic acid sequences (e.g., whole genome sequences), which provide, for example, information about a characteristic or medical condition of a patient or embryo or fetus, such as the sex of the embryo or fetus, or the presence or prognosis of a disease with a genetic component (including, for example, cystic fibrosis, sickle cell anemia, Marfan syndrome, Huntington's disease, and hemochromatosis, or various cancers, such as breast cancer). According to another embodiment, the sequencing methods of the present invention are used to provide sequence information starting with 1-20 cells from a patient (including, but not limited to, a fetus or embryo) and assessing the patient's characteristics based on the sequence.

Cancer Diagnostics

Whole genome sequencing is a valuable tool in assessing the genetic basis of disease. Many diseases (such as cystic fibrosis) are known to have a genetic basis.

One application of whole genome sequencing is to understand cancer. The most important impact of next-generation sequencing on cancer genomics is the ability to resequence, analyze, and compare matched tumor and normal genomes of a single patient and multiple patient samples of a given cancer type. Using whole genome sequencing, a full range of sequence variations can be considered, including germline susceptibility loci, somatic single nucleotide polymorphisms (SNPs), small insertion and deletion (indel) mutations, copy number variations (CNVs), and structural variants (SVs).

Typically, a cancer genome is composed of the patient's germline DNA onto which somatic genomic changes have been superimposed. Somatic mutations identified by sequencing can be classified as "driver" or "passenger" mutations. So-called driver mutations are those that directly contribute to tumor progression by conferring a growth or survival advantage on the cell. Passenger mutations encompass neutral somatic mutations that have been acquired during errors in cell division, DNA replication, and repair; these mutations can be acquired when the cell is phenotypically normal or after neoplastic changes are apparent.

Historically, attempts have been made to elucidate the molecular mechanisms of cancer, and several "driver" mutations or biomarkers have been identified, such as HER2/neu2. Based on this type of gene, therapeutic regimens have been developed to specifically target tumors with known genetic changes. The best limiting example of this approach is the targeting of HER2/neu by trastuzumab (Herceptin) to breast cancer cells. However, cancer is not a simple single-cause disease, but rather a combination of genetic changes that can be different between individuals. Therefore, these other interferences to the genome can make some drug regimens become invalid for some individuals.

Cancer cells for whole genome sequencing can be obtained from a whole tumor biopsy (including microbiopsies of a small number of cells), cancer cells isolated from a patient's bloodstream or other bodily fluid, or any other source known in the art.

Preimplantation genetic diagnosis

One application of the method of the present invention is for preimplantation genetic diagnostics. About 2 to 3% of babies are born with some type of major birth defect. The risk of some problems due to genetic segregation of genetic material (chromosomes) increases with maternal age. About 50% of the chances of these types of problems are due to Down syndrome, which is a third copy of chromosome 21 (trisomy 21). The other half are due to other types of chromosomal abnormalities, including trisomy, point mutations, structural variations, copy number variations, and the like. Many of these chromosomal problems result in severely affected babies or even those that do not survive to delivery.

In medicine and (clinical) genetics, preimplantation genetic diagnosis (PGD or PIGD) (also known as embryo screening) refers to a procedure performed on embryos before implantation, and sometimes even on oocytes before fertilization. PGD can allow parents to avoid selective pregnancy termination. The term preimplantation genetic screening (PGS) is used to refer to a procedure that does not look for a specific disease, but uses PGD technology to identify embryos that are at risk due to, for example, a genetic condition that can cause a disease. The procedure performed on sex cells before fertilization may instead be referred to as a method of oocyte selection or sperm selection, although the method and purpose partially overlap with PGD.

Preimplantation genetic profiling (PGP) is a kind of assisted reproductive technology to implement the method for embryo selection, and described embryo seems to have the maximum chance of successful pregnancy.When being used for the women of late maternal age and for the patient of repeated in vitro fertilization (IVF) failure, mainly implement PGP as for detecting chromosomal abnormality such as aneuploidy, reciprocal translocation and Robertsonian translocation and other abnormal screening such as chromosomal inversion or deletion.In addition, PGP can be to genetic marker inspection feature, including multiple disease states.PGP uses the principle of back, because known many chromosomal genetic explanations most of pregnancy loss cases, and the human embryo of larger proportion is aneuploidy, and the selective replacement of euploid embryo should improve the chance of successful IVF treatment.Whole genome sequencing provides comprehensive chromosome analysis method, such as array comprehensive genomic hybridization (aCGH), quantitative PCR and SNP microarray methods alternative.For example, whole whole genome sequencing can provide information about single base variation, insertion, deletion, structural variation and copy number variation.

Since PGD can be performed on cells from different developmental stages, the biopsy protocol varies accordingly. Biopsy can be performed at all preimplantation stages, including but not limited to unfertilized and fertilized oocytes (for polar bodies, PBs), three-day cleavage stage embryos (for blastomeres), and blastocysts (for trophectoderm cells).

In view of the above detailed description of the present invention, according to one aspect of the present invention, a method for sequencing complex nucleic acids of an organism (e.g., a mammal such as a human, whether a single organism or a population comprising more than one individual) is provided, such method comprising: (a) aliquoting a sample of the complex nucleic acid to produce a plurality of aliquots, each aliquot containing a certain amount of the complex nucleic acid; (b) sequencing the amount of the complex nucleic acid from each aliquot to produce one or more read results from each aliquot; and (c) assembling the read results from each aliquot to produce an assembled sequence of the complex nucleic acid, which contains no more than 1, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.08, 0.06, 0.04 or less pseudo single nucleotide variants per megabase at a response rate of 70, 75, 80, 85, 90 or 95% or greater. If the complex nucleic acid is a mammalian (eg, human) genome, optionally, the assembled sequence has a genome call rate of 70% or greater and an exome call rate of 70, 75, 80, 85, 90, or 95% or greater. According to one embodiment, the complex nucleic acid comprises at least 1 gigabase.

According to one embodiment of such methods, the complex nucleic acid is double-stranded, and the method comprises separating the single strands of the double-stranded complex nucleic acid prior to aliquoting.

According to another embodiment, such methods include fragmenting the amount of complex nucleic acid in each aliquot to generate fragments of the complex nucleic acid. According to one embodiment, such methods further include tagging the fragments of the complex nucleic acid in each aliquot with an aliquot-specific tag (or aliquot-specific set of tags), wherein the aliquot-specific tag (or aliquot-specific set of tags) allows identification of the aliquot from which the tagged fragments originated. In one embodiment, such tags are polynucleotides, including, for example, tags comprising an error correction code or error correction codes, including but not limited to Reed-Solomon error correction codes.

According to another embodiment, such methods comprise pooling the aliquots prior to sequencing.

According to another embodiment of such methods, the sequence comprises base calls at sequence positions, and such methods comprise identifying a base call as genuine if it originates from two or more aliquots, or from three or more reads originating from two or more aliquots.

According to another embodiment, such methods include identifying a plurality of sequence variants in the assembled sequence and phasing the sequence variants.

According to another embodiment of such methods, the sample of complex nucleic acid comprises 1 to 20 cells of an organism or genomic DNA isolated from cells, which may be purified or unpurified. According to another embodiment, the sample comprises 1 pg to 100 ng, such as 1 pg, 6 pg, 10 pg, 100 pg, 1 ng, 10 ng, or 100 ng of genomic DNA, or 1 pg to 1 ng, or 1 pg to 100 pg, or 6 pg to 100 pg. For reference purposes, a single human cell contains approximately 6.6 pg of genomic DNA.

According to another embodiment, such methods comprise amplifying said amount of complex nucleic acid in each aliquot.

According to another embodiment of such methods, the complex nucleic acid is selected from the group consisting of a genome, an exome, a transcriptome, a methylome, a mixture of genomes of different organisms, a mixture of genomes of different cell types of an organism, and subsets thereof.

According to another embodiment of such a method, the assembled sequence has a coverage of 80x, 70x, 60x, 50x, 40x, 30x, 20x, 10x, or 5x. Lower coverage can be used with longer reads.

According to another aspect of the present invention, an assembled sequence of a mammalian complex nucleic acid is provided that contains less than 1 spurious single nucleotide variant per megabase at a call rate of 70% or greater.

According to another aspect of the present invention, a method for sequencing complex nucleic acids from an organism is provided, the method comprising: (a) providing a sample comprising 1 pg to 10 ng of complex nucleic acids; (b) amplifying the complex nucleic acids to generate amplified nucleic acids; and (c) sequencing the amplified nucleic acids to generate a sequence having a complex nucleic acid call rate of at least 70%. According to one such method, the complex nucleic acids are unpurified. According to another embodiment, the method comprises amplifying the complex nucleic acids by multiple displacement amplification. According to another embodiment, the method comprises amplifying the complex nucleic acids by at least 10, 100, 1000, 10,000, or 100,000-fold or more. According to another embodiment of the method, the sample comprises 1 to 20 cells (or cell nuclei) comprising complex nucleic acids. According to another embodiment, the method comprises lysing cells (or nuclei) comprising complex nucleic acids and cellular impurities, and amplifying the complex nucleic acids in the presence of the cellular impurities. According to another embodiment of the method, the cells are circulating non-blood cells from the blood of a higher organism. According to another embodiment of such a method, the assembled sequence has a call rate of 70, 75, 80, 85, 90 or 95% or more. According to another embodiment of such a method, the sequence comprises 2, 1, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.08, 0.06, 0.04 or less pseudo single nucleotide variants per megabase. According to another embodiment, such a method further comprises: aliquoting the sample to generate a plurality of aliquots, each aliquot comprising an amount of complex nucleic acid; amplifying the amount of complex nucleic acid in each aliquot to generate amplified nucleic acid in each aliquot; sequencing the amplified nucleic acid from each aliquot to generate one or more reads from each aliquot; and assembling the reads to generate a sequence. According to another embodiment, such methods further include: fragmenting the amplified nucleic acid in each aliquot to generate fragments of amplified nucleic acid in each aliquot; and tagging the fragments of amplified nucleic acid in each aliquot with aliquot-specific tags to generate tagged fragments in each aliquot. According to another embodiment of such methods, a base call at a sequence position is accepted as true if it is present in reads from two or more aliquots, or more strictly, occurs three or more times in reads from two or more aliquots. According to another embodiment, such methods further include identifying sequence variations in the sequence that provide information about a characteristic of the organism (e.g., a medical condition). According to another embodiment, the cells are circulating non-blood cells from the blood (or other sample) of a higher organism, including but not limited to fetal cells from maternal blood and cancer cells from the blood of a patient with cancer. According to another embodiment of the present invention, the complex nucleic acid is circulating nucleic acid (CNA). Thus, characteristics of the organism to be assessed can include but are not limited to the presence of and information about cancer (regardless of whether the organism is pregnant), and information about the sex or genetics of a fetus carried by the pregnant individual. For example, such methods can be used to identify single base variations, insertions, deletions, copy number variations, structural changes, or rearrangements associated with disease likelihood, medical diagnosis, or prognosis, etc. According to another embodiment of the present invention, a method for assessing the genetic status of an embryo (e.g., sex, parentage, the presence or absence of a genetic abnormality, or a genotype associated with a disease predisposition, etc.) is provided, comprising: (a) providing approximately 1-20 embryonic cells; (b) obtaining an assembly sequence generated by sequencing the genomic DNA of the cells, wherein the assembly sequence has a response rate of at least 80%; and (c) comparing the assembly sequence with a reference sequence to assess the genetic status of the embryo. For example, such methods can be used to identify single base variations, insertions, deletions, copy number variations, structural changes, or rearrangements associated with disease likelihood, medical diagnosis, or prognosis, etc. According to another embodiment, a method for assessing the genetic status of an embryo (e.g., sex, parentage, presence or absence of a genetic abnormality or a genotype associated with a disease predisposition, etc.) is provided, comprising: (a) providing approximately 1-20 embryonic cells; (b) obtaining an assembled sequence generated by sequencing the genomic DNA of the cells, wherein the assembled sequence has a response rate of at least 80% of the embryonic genome; and (c) comparing the assembled sequence with a reference sequence to assess the genetic status of the embryo.

According to another aspect of the present invention, an assembled whole human genome sequence is provided, comprising no more than 1 false single nucleotide variant per megabase and a call rate of at least 70%, wherein the sequence is generated by sequencing 1 pg-10 ng of human genomic DNA.

According to another aspect of the present invention, a method for phasing genomic sequence variants of an individual organism comprising a plurality of chromosomes is provided, the method comprising: (a) providing a sample comprising a mixture of vector-free fragments of each of the plurality of chromosomes; (b) sequencing the vector-free fragments to generate a genomic sequence comprising a plurality of sequence variants; and (c) phasing the sequence variants. According to one embodiment, such a method comprises phasing at least 70, 75, 80, 85, 90, or 95% or more of the sequence variants. According to another embodiment of such a method, the genomic sequence has a call rate of at least 70% of the genome. According to another embodiment of such a method, the sample comprises 1 pg to 10 ng of genome, or 1 to 20 cells of the individual organism. According to another embodiment of such a method, the genomic sequence has less than 1 false single nucleotide variant per megabase.

According to another aspect of the present invention, a method for phasing genomic sequence variants of an individual organism comprising multiple chromosomes is provided, the method comprising: providing a sample comprising fragments of the multiple chromosomes; sequencing the fragments without cloning the fragments in vectors to generate a whole genome sequence, wherein the whole genome sequence comprises multiple sequence variants; and phasing the sequence variants. According to one embodiment of such a method, phasing the sequence variants occurs during assembly of the whole genome sequence.

Example

Example 1: Comparison of DNA amplification methods

Preimplantation genetic diagnosis (PGD) is a form of prenatal diagnostics that consists of embryos (typically 10 per cycle on average) produced by genetic screening in vitro fertilization (IVF) and then transferring them to future mothers. It is generally suitable for women in late maternal age (greater than 34 years old) or couples with the risk of transmitting genetic diseases. The technology currently used for genetic screening is fluorescence in situ hybridization (FISH), comparative genomic hybridization (CGH), SNP arrays and array CGH for detecting chromosomal abnormalities, and SNP arrays and PCR for detecting gene defects. PGD for single-gene defects currently consists of a unique custom-designed assay for each patient, which often combines specific mutation detection with linkage analysis as a backup and controls and monitors contamination. Typically, 1 cell is obtained from each embryo biopsy on the 3rd day of development, and results are given on the 5th day (which is the latest day that embryos can be transferred). Blastocyst biopsy is initially used, which consists of a biopsy of 3-15 cells from the trophectoderm of the blastocyst (5th day embryo), followed by embryo freezing. Embryos can be kept frozen indefinitely without significant loss of potency and are suitable for whole genome sequencing, allowing a biopsy to be obtained at one site and then transferred to another site for whole genome sequencing. Whole genome sequencing of blastocyst biopsies will enable "universal" PGD testing for single gene defects and other genetic abnormalities that can be identified by this technology.

After conventional ovarian stimulation and egg retrieval, ovum is fertilized to avoid the sperm contamination in the PGD test by intracytoplasmic sperm injection (ICSI).After growing to the 3rd day, fine glass needle biopsy is used to obtain embryo, and one cell is taken out from each embryo.Each blastomere is added to clean pipe separately, covered with molecular grade oil, and transported to PGD laboratory on ice.After arrival, use the test treatment sample that is designed for the sudden change that CTG repeats expansion and two linkage markers in the amplification gene DMPK immediately.

After clinical PGD testing and embryo transfer, unused embryos were donated to an IVF clinic and used in the development of new PGD test formats. Eight blastocysts were donated and used in these experiments.

A blastocyst biopsy provides approximately 6.6 picograms (pg) of genomic DNA per cell. Amplification provides enough DNA for whole-genome sequencing. Figure 15 shows the results of MDA amplification of 1.031 pg, 8.25 pg, and 66 pg of purified genomic DNA standards and 1 or 10 PVP40 cells using our protocol (described below). The MDA reaction can be run for as long as necessary to obtain the amount of DNA required for a particular sequencing method (e.g., 30 to 120 minutes). It is expected that greater amplification will produce more GC bias.

Two DNA amplification methods were compared to identify a method that generates template DNA of sufficient quality for whole genome sequence analysis while minimizing the introduction of GC bias.We compared our protocol with SurePlex amplification (Rubicon Genomics Inc., Ann Arbor, Michigan) and a modified MDA commonly used for array CGH.

A biopsy of 10-20 cells was obtained from an embryo affected by the R-1MT mutation of myotonic dystrophy. The sample was lysed and the DNA denatured in a single tube, then amplified by MDA using our protocol and the SurePlex kit according to the manufacturer's instructions. Approximately 2 ug of DNA was generated by these two amplification methods. Before whole-genome sequence analysis, the amplified samples were screened with 96 independent qPCR markers scattered across the genome to select samples with the lowest amount of preference. Figure 16 shows the results. In brief, we determined the average number of cycles across the entire plate and deducted this number from each individual marker to calculate the "Δ cycle" number. The Δ cycle was plotted relative to the GC content of 1000 base pairs around each marker to indicate the relative GC preference of each sample. In order to clarify the overall "noise" of the sample, the absolute value of each Δ cycle was summed to produce a "Δ sum" measure. A lower Δ sum and a relatively flat data plot relative to GC content produce well-presented whole-genome sequences in our experience. The sum of the deltas is 61 (for our MDA method) and 287 (for SurePlex amplified DNA), indicating that our protocol produces much less GC bias than the SurePlex protocol.

Example 2: Complete genome sequencing of blastocyst biopsies for preimplantation genetic diagnosis (PGD)

Adopt the multiple displacement amplification (MDA) (Dean et al. (2002) Proc Natl Acad Sci U S A99,5261-5266) of modification to generate enough template DNA (about 1 μ g) for whole genome sequence analysis, as described herein. In brief, 5-20 cells from each 5-day-old embryonic cell are separated, frozen, and on dry ice, transported from the laboratory that separates them. The sample is melted, and cracked to release genomic DNA. In the situation that genomic DNA is not purified and is kept away from cell impurities, DNA is alkaline denatured by adding 1 μ l 400 mM KOH/10 mM EDTA. Use the multiple displacement amplification (MDA) reaction based on phi29 polymerase to carry out whole genome amplification to generate enough DNA (about 1 μ g) to carry out sequencing to embryonic genomic DNA. After 1 minute of alkaline denaturation, the denatured DNA is added with random 8 polymers protected by thiol. After 2 minutes, the mixture was neutralized and a master mix containing a final concentration of 50 mM Tris-HCl (pH 7.5), 10 mM MgCl2, 10 mM (NH4)2SO4, 4 mM DTT, 250 μM dNTPs (USB, Cleveland, OH) and 12 units of phi29 polymerase (Enzymatics, Beverly, MA) was added to yield a total reaction volume of 100 μl. The MDA reaction was incubated at 37°C for 45 minutes and inactivated at 65°C for 5 minutes. Approximately 2 μg of DNA was generated by the MDA reaction. This amplified DNA was then fragmented and used for library construction and sequencing as described above.

Myotonic dystrophy type 1 (DM1) is an autosomal dominant disease caused by a trinucleotide repeat expansion, i.e., cytosine-thymine-guanine (CTG)n, in the 3' untranslated region of the gene encoding myotonic dystrophy protein kinase (DMPK). We examined the clonal coverage of the DMPK CTG repeat interval. The sequencing technology described herein results in 35bp paired-end reads, which typically span approximately 400bp. For unaffected individuals and an unknown sample, 400bp is sufficient to span this CTG repeat region of both alleles, resulting in a copy number of approximately 2. In affected individuals and an unknown sample, a copy number of approximately 1 was observed, suggesting that the repeat expansion is too large for the 400bp paired-end to span; only the unaffected allele has coverage in this region.

Table 1 below provides summary information for positioning and assembling PGD embryo samples. All variation and positioning statistics are relative to the National Center for Biotechnology Information (NCBI) 37th edition human genome reference assembly. Samples 2A, 5B and 5C had poor amplification quality, resulting in less genomic response and a decrease in the total number of SNPs identified. Samples 5B and 5C are from different biopsies of the same embryo. Sample NA20502 was processed according to standard procedures and not amplified before library preparation.

Figure 17 shows the genome coverage of two samples (7C and 10C).Use the 10 megabase moving average of the 100 kilobase coverage window of standardization relative to the haploid genome to cover the coverage drawing.The dotted lines of copy number 1 and 3 represent haploid and triploid copy number respectively.These two embryos are males, and have haploid copy number for X and Y chromosome.Other loss or acquisition without full chromosome or chromosome large segment is obvious in these samples.

The worst performing sample achieved 85% genome coverage, while the best sample covered 95% of the genome, a level similar to that of a standard whole genome sequencing method ("standard sequencing") performed using the method described above using a few micrograms of purified, unamplified human genomic DNA. Generally, coverage is "noisy" compared to standard sequencing, but using a moving average of 10 megabases allows accurate detection of whole genome and chromosome arm amplifications and deletions. We also demonstrated that many polymorphisms can be detected, and that, in addition to DMPK mutations, the risk of developing certain diseases can be used for blastocyst implantation selection.

In this embodiment, the initial genomic DNA is extensively amplified (more than necessary approximately 10 times) to ensure that a sufficient amount of genomic DNA can be used for sequencing. Expectedly, reducing the degree of amplification improves sequence coverage and sequencing quality. Amplification can also be reduced by allowing the tissue obtained by biopsy (or other starting materials, such as cancer biopsy or needle aspirates, fetuses or cancer cells separated from blood flow, etc.) to grow in culture. This method slightly increases the overall turnaround time of the method. However, cultivating a small amount of available cells results in high-fidelity " amplification" of genomic DNA in the cell process of chromosome replication.

Because DMPK mutation is a trinucleotide repeat disease, it is difficult to analyze the mutation using current sequencing methods that use mate pair reads of approximately 400 bp in length. Longer mate pair reads (e.g., 1 kilobase or longer) can be used to span these regions and therefore sequence between them, which leads to accurate determination of repeat size.

Example 3: Clinically accurate genome sequencing and haplotype determination from 10-20 human cells

In this example, 65-130 pg (10-20 cells) of human genomic DNA (50% 60-500 kb in length) was divided into 384 aliquots, amplified, fragmented, and tagged in each aliquot. After sequencing, the diploid (phased) genome was assembled without DNA cloning or separation of metaphase chromosomes. Ten LFR libraries were used to generate approximately 3.3 terabases (Tb) of mapping reads from 7 unique genomes. Up to 97% of heterozygous single nucleotide variants (SNVs) were assembled into contigs, with 50% of the coverage bases (N50) being in contigs longer than about 500 kb (for European ethnic samples) and about 1 Mb (for African samples). In extensive comparisons between replicate libraries, the LFR unit type was found to be highly accurate, with 1 false positive SNV per 10 megabases (Mb). Despite starting with 100 picograms (pg) of DNA and 10,000-fold in vitro amplification, this 20-30-fold increase in accuracy compared to non-LFR genomes was achieved (Drmanac et al., Science 327:78, 2010; Roach et al., Am. J. Hum. Genet. 89:382-397, 2011), because most errors were inconsistent with the true haplotype, we have demonstrated cost-effective and clinically accurate genome sequencing and haplotype determination from 10-20 human cells.

LFR technology is a cost-effective DNA pre-processing step without cloning or whole metaphase chromosome isolation, which allows complete sequencing and assembly of different parental chromosomes at a clinically relevant cost and scale. LFR can be used as a pre-processing step before any sequencing method, although we used short-read sequencing technology, as described in detail above.

LFR can produce long range phased SNPs because it is conceptually similar to single molecule sequencing of fragments with a length of 10-1000kb. This is achieved by randomly dividing the corresponding parental DNA fragments into physically unique sets in the absence of any DNA cloning step, followed by fragmentation to generate shorter fragments (this is similar to the aliquot sampling of fosmid clones (Kitzman et al., Nat. Biotechnol. 29: 59-63, 2011; Suk et al., Genome Res. 21: 1672-1685, 2011)). Due to the fact that the fraction of genomes in each set is reduced to less than the haploid genome, the statistical probability of having corresponding fragments from the two parental chromosomes in the same set is significantly reduced. Similarly, the more single sets are inquired, the larger the number of times the fragments from maternal and paternal homologs are analyzed in different sets.

For example, a 384-well plate with 0.1 genome equivalents in each well produces a theoretical 19x coverage of both the maternal and paternal alleles of each fragment. This type of high initial DNA redundancy of approximately 19x produces a more complete genome coverage and higher variant response and phasing accuracy than using a strategy employing a fosmid collection (which results in a range of approximately 3x (Kitzman et al., Nat. Biotechnol 29:59-63, 2011) to approximately 6x (Suk et al., Genome Res. 21:1672-1685, 2011) coverage).

In order to prepare LFR libraries in a high-throughput manner, we developed an automated method that implements all LFR-specific steps in the same 384-well plate. The following is an overview of the method. First, a modified phi29-based multiple displacement amplification (MDA; Dean et al., Proc. Natl. Acad. Sci. U.S.A. 99: 5261, 2002) is used to implement highly consistent amplification to replicate each fragment approximately 10,000 times. Next, the DNA is fragmented and connected to the barcode adapters via an enzymatic step process in each well without an intervening purification step. In brief, long DNA molecules are processed into blunt-ended 300-1,500bp fragments by controlled random enzymatic fragmentation (CoRE). CoRE fragments the DNA by removing uridine bases, which are incorporated at a predetermined frequency during the MDA process by uracil DNA glycosylase and endonuclease IV. The fragments were resolved by nick translation from the resulting single-base gaps using E. coli polymerase 1, and blunt ends were produced. Then, unique 10-base Reed-Solomon error-corrected barcode adapters (PCT/US2010/023083, published as WO 2010/091107, incorporated herein by reference) (designed to reduce any bias (Figure 18) caused by sequence and concentration differences of each barcode) were connected to fragment the DNA in each well using a high-yield, low-chimera formation scheme (Drmanac et al., Science 327:78, 2010). Finally, all 384 wells were combined and unsaturated polymerase chain reaction was employed to generate templates sufficient for a short read sequencing platform using primers common to the connection adapters. More details about the LFR scheme we adopted are provided below.

High molecular weight DNA was purified from cell lines GM12877, GM12878, GM12885, GM12886, GM12891, GM12892 GM19240, and GM20431 (Coriell Institute for Medical Research, Camden, NJ) using the RecoverEase DNA Isolation Kit (Agilent, La Jolla, CA) following the manufacturer's protocol. The high molecular weight DNA was partially sheared to make it more suitable for manipulation using a Rainin P1000 pipette by pipetting 20-40 times. 200 ng of genomic DNA was analyzed on a 1% agarose gel with 0.5X TBE buffer using a BioRad CHEF-DR II using the following parameters: 6 V/cm, 50-90 seconds gradient transition time, and a 20-hour total run. The length of purified genomic DNA was determined using 500 ng of Yeast Chromosomal PFG Marker (New England Biolabs, Ipswich, MA) and Lambda Ladder PFG Marker (New England Biolabs, Ipswich, MA).

Separately, the immortalized cell line GM19240 (Coriell Institute for Medical Research, Camden, NJ) was cultured in RPMI supplemented with 10% FBS under standard cell culture conditions. Single cells were isolated at 200x magnification using a micromanipulator (Eppendorf, Hamburg, Germany) and placed into 1.5 ml microtubes containing 10 μl of dHO. The cells were denatured with 1 μl of 20 mM KOH and 0.5 mM EDTA. The denatured cells were then subjected to the LFR process.

DNA from each of the multiple cell lines was diluted and denatured at a concentration of 50 pg/ul in a 20 mM KOH and 0.5 mM EDTA solution. After incubation at room temperature for 1 minute, 120 pg of denatured DNA was removed and added to 32 ul of 1 mM 3'thiol-protected random octamer (IDT, Coralville, IA). After 2 minutes, the mixture was brought to a volume of 400 ul with dH2O, and 1 ul was dispensed into each well of a 384-well plate. 1 μl of 2X multiple displacement amplification (MDA) mix based on phi29 polymerase (Enzymatics Inc., Beverly, MA) was added to each well to generate approximately 3-10 ng of DNA (10,000 to 25,000-fold amplification). The MDA reaction consisted of 50 mM Tris-HCl (pH 7.5), 10 mM MgCl2, 10 mM (NH4)2SO4, 4 mM DTT, 250 uM dNTP (USB, Cleveland, OH), 10 uM 2'-deoxyuridine 5'-triphosphate (dUTP) (USB, Cleveland, OH), and 0.25 units of phi29 polymerase.

Then, controlled random enzymatic fragmentation (CoRE) was performed. Excess nucleotides were inactivated and uracil bases were removed by incubating the MDA reaction with a mixture of 0.031 units of shrimp alkaline phosphatase (SAP) (USB, Cleveland, OH), 0.039 units of uracil DNA glycosylase (New England Biolabs, Ipswich, MA), and 0.078 units of endonuclease IV (New England Biolabs, Ipswich, MA) at 37°C for 120 minutes. SAP was heat-inactivated at 65°C for 15 minutes. Gaps were resolved by 60 minutes of room temperature nick translation with 0.1 units of Escherichia coli DNA polymerase 1 (New England Biolabs, Ipswich, MA) in the same buffer with 0.1 nanomolar dNTPs (USB, Cleveland, OH), and the DNA was fragmented into 300-1,300 base pair fragments. E. coli DNA polymerase 1 was heat-inactivated at 65°C for 10 minutes. Remaining 5' phosphate groups were removed by incubation with 0.031 units of SAP (USB, Cleveland, OH) for 60 minutes at 37°C. SAP was heat inactivated at 65°C for 15 minutes.

Then, tagging adapter ligation and nick translation were performed. 10 base DNA barcode adapters (unique for each well) were attached to the fragmented DNA using a two-part directional ligation method. Approximately 0.03 pmol of fragmented MDA product was incubated for 4 hours at room temperature in a reaction in a total volume of 7 ul containing 50 mM Tris-HCl (pH 7.8), 2.5% PEG 8000, 10 mM MgCl2, 1 mM rATP, a 100-fold molar excess of 5'-phosphorylated (5'PO4) and 3' dideoxy-terminated (3'dd) common Ad1 (Figure 18) and 75 units of T4 DNA ligase (Enzymatics, Beverly, MA). Ad1 contains a common overhang region for ligation and hybridization with unique barcode adapters. After 4 hours, a 200-fold molar excess of unique 5' phosphorylated tagging adapters was added to each well and allowed to incubate for 16 hours. 384 wells were assembled into a total volume of approximately 2.5 ml and purified by adding 2.5 ml of AMPure beads (Beckman-Coulter, Brea, CA). One round of PCR was performed to create molecules with a 5' adapter and tag on one side and a 3' blunt end on the other side. The 3' adapter was added in a ligation reaction similar to the 5' adapter, as described above. To seal the nick created by ligation, the DNA was incubated at 60°C for 5 minutes in a reaction containing 0.33 uM Ad1 PCR1 primer, 10 mM Tris-HCl (pH 78.3), 50 mM KCl, 1.5 mM MgCl2, 1 mM rATP, 100 uM dNTPs to exchange the 3' dideoxy-terminated Ad1 oligomer with the 3'-OH-terminated Ad1 PCR1 primer. The reaction was then cooled to 37°C and incubated for an additional 30 minutes at 37°C after the addition of 90 units of Taq DNA polymerase (New England Biolabs, Ipswich, MA) and 21,600 units of T4 DNA ligase to create a functional 5'-PO4 gDNA end from the 3'-OH end of the Ad1 PCR1 primer by Taq-catalyzed nick translation and to seal the resulting repaired nick by T4 DNA ligation. At this point, the material was incorporated into a standard DNA nanoarray sequencing method.

Starting from total RNA, RNA-Seq data were obtained using the Ovation RNA-Seq kit (NuGen, San Carlos, CA) and SPRIWork (Beckman-Coulter, Brea, CA) to prepare sequencing libraries with an average insert size of 150-200 bp. 75 bp paired-end sequencing reactions were performed on a HiSeq 2000 (Illumina, San Diego, CA) at the Center for Personalized Genetic Medicine (Harvard Medical School, Boston, MA). Paired-end reads were assembled using bowtie v0.12.7 (Langmead et al., Genome Biol. 10: R25, 2009) with tophat v1.2.0 (Trapnell et al., Bioinformatics 25: 1105-1111, 2009), and single nucleotide variants (SNVs) were mapped using GATK UnifiedGenotyper v1.1 (http://www.broadinstitute.org/gsa/wiki/index.php/GATK_release_1.1) with reference hg19 and dbSNP version 132 for annotating known SNPs. SNVs were also mapped to genes from RefSeq and isoforms in transcriptomes identified using cufflinks v1.0.3 (http://cufflinks.cbcb.umd.edu/tutorial.html).

In order to identify the haplotypes of co-expressed alleles, data on heterozygous SNVs were filtered that occurred simultaneously on the same LFR contig and on the same gene with at least one other heterozygous SNV. In the case where transcripts exhibit allele-specific expression, the heterozygous alleles expressed on the LFR phased haplotype should have a higher or lower read count than their counterparts on the other haplotype. Here, we identify the haplotypes with higher expression as haplotypes in which most of their heterozygous alleles exhibit higher expression than their counterparts. If its expression is consistent with the haplotype it contains, the heterozygous calculation is "consistent". In the case of a split (where there is no haplotype majority), half of the heterozygous SNVs are calculated as consistent. In addition, in order to be fully considered, heterozygous SNVs are required to have at least 20 times RNA-Seq read coverage. The noise from the GATK genotyper was further filtered for heterozygous SNVs by randomly using a binomial test with the probability of selecting ASE and coverage.

For error correction purposes, each DNB was tagged with a 10-base Reed-Solomon code with either 1-base error correction capability for unknown error positions or 2-base error correction capability when the error position was known (U.S. patent application Ser. No. 12/697,995, published as US2010/0199155, incorporated herein by reference). These 384 codes were selected from a comprehensive set of 4096 Reed-Solomon codes with the properties described above (U.S. patent application Ser. No. 12/697,995, incorporated herein by reference). Each code from this set had a minimum Hamming distance of 3 to any other code in the set. For this study, the error position was assumed to be unknown.

Results . To demonstrate the ability of LFR to determine accurate diploid genome sequences, we generated three libraries of the Yoruba female HapMap sample NA19240. NA19240 was extensively interrogated as part of a triplicate study (NA19240 is a descendant of samples NA19238 and NA19239) in the HapMap Project (Consortium, Nature 437:1299-1320, 2005; Frazer et al., Nature 449:851-861, 2007), the 1,000 Genomes Project (Nature 467:1061-1073, 2010), and our own efforts (www.completegenomics.com/sequence-data/download-data/). Thus, based on the redundant sequence data of the parental samples NA19238 and NA19239, highly accurate haplotype information for 1.7 million heterozygous SNPs could be generated. Starting with 10 cells (65pg DNA) of the corresponding immortalized B cell line, one NA19240 LFR library was generated. Based on the total effective read coverage of 60x and the use of 384 unique fragment aliquots or collections, we estimated that if the DNA was denatured before being distributed to the wells (20 cell equivalents of dsDNA; Table 1 below), the optimal number of starting cells would be 10. Two replicate libraries were generated from an estimated 100-130pg (15-20 cell equivalents) of denatured high molecular weight genomic DNA. It was determined that the optimal amount of each library would be about 100pg when starting from denatured isolated DNA. This amount was selected to achieve more consistent genome coverage by minimizing random sampling of the sample.

All three libraries were analyzed using DNA nanoarray sequencing (Drmanac et al., Science 327:78-81, 2010). Custom alignment algorithms (Drmanac et al., Science 327:78-81, 2010; Carnevali et al., J. Computational Biol., 19, 2011) were used to map the mate pair reads of 35 bases to the reference genome, generating an average of more than 230Gb of positioning data (Table 1 below) with an average genome coverage greater than 80x. Analysis of the positioning LFR data showed two unique features attributable to MDA: a slight underrepresentation of GC-rich sequences (Figure 19) and an increase in chimeric sequences. Additionally, the variability of the standardized coverage between 100kb windows was approximately 2 times greater. However, nearly all genomic regions were covered with enough reads (5 or more), indicating that the 10,000 times of MDA amplification performed by our optimization scheme can be used for comprehensive genome sequencing.

Barcodes were used to graphically group the reads based on their physical well positions within each library (which showed pulses of coverage, i.e., sparse regions of coverage interspersed between long spans with little to no read coverage). On average, each well contained 10-20% of the haploid genome (300-600 Mb) in fragments ranging in length from 10 kb to more than 300 kb, with an N50 of approximately 60 kb (Figure 20). Initial fragment coverage was very consistent across chromosomes. As assessed from all detected fragments, the total amount of DNA actually used to generate the two libraries from the extracted DNA was approximately 62 pg and 84 pg (9.4 and 12.7 cell equivalents, Figure 20). This is less than the expected 100-130 pg, indicating some loss or undetectable DNA or imprecision in DNA quantification. Interestingly, the library of 10 cells appeared to be generated from approximately 90 pg (13.6 cells) of DNA, most likely due to some cells being in the S phase during separation (Figure 20).

Using a two-step custom genotyping algorithm designed to query low coverage read result data (less than 2x coverage) from about 40 single wells, overlapping heterozygous SNPs from fragments of the same parent chromosome located in different wells are assembled into unit type contigs (Figure 21). Unlike other experimental methods (Kitzman et al., Nat. Biotechnol. 29: 59-63, 2011; Suk et al., Genome Res. 21: 1672-1685, 2011; Duitama et al., Nucl. Acids Res. 40: 2041-2053, 2012), LFR does not limit the unit type of each initial fragment. Instead, LFR ensures complete presentation of the genome by maximizing DNA fragment input in terms of the number of aliquots and given read result coverage.

In the first step, the heterozygous SNPs from the unphased NA19240 genome assembly (www.completegenomics.com/sequence-data/download-data/) are combined with each LFR library to create a comprehensive SNP group for phasing. Then, a network is constructed for each chromosome, in which nodes correspond to heterozygous SNP responses, and connections relate to the connectivity score between each pair of SNPs. Together with the connection score, a direction is also obtained as part of the search for the best hypothesis about each pair of heterozygous SNPs. This highly redundant sparse connection network is then trimmed using domain knowledge and subsequently optimized using Kruskal's minimum spanning tree (MST) algorithm. This produces longer overlapping groups, with N50 from 950-1200kb available from these libraries (Figure 20).

By LFR phasing in each library a total of about 2.4 million heterozygous SNPs (Figure 20). LFR phasing is expected to gradually adopt about 90% of the heterozygous SNPs of these libraries. The library phasing of 10 cells is more than 98% of the variants of the two library phasings generated by self-isolated DNA, proving the potentiality that LFR works by a small amount of isolated cells. The number of read results is doubled to about 160x coverage and the number of phased heterozygous SNPs is further increased to more than 2.58 million, thus increasing the phasing rate to 96% (Figure 20). Combination repeats 1 and 2 (768 independent holes in total) (each with 80x coverage) produce more than 2.65 million phased heterozygous SNPs, and produce 97% phasing rate. Only using the SNP loci (omitting step 1 of the LFR algorithm) that responds in the LFR library for phasing usually causes the phased SNP sum to reduce by 5-15% (Figure 20).

Importantly, the number of phased SNPs obtained by LFR alone (starting with DNA from only 10-20 cells) is slightly higher than the number of SNPs phased by current fosmid methods (Kitzman et al., Nat. Biotechnol. 29:59-63, 2011; Suk et al., Genome Res. 21:1672-1685, 2011; Duitama et al., Nucl. Acids Res. 40:2041-2053, 2012). Because parents share a large fraction of variants in their children, this is substantially more than 81% of the heterozygous SNPs that can be phased using standard parental sequences (Roach et al., Am. J. Hum. Genet. 89:382-397, 2011). Adding parent-derived haplotype data to the 768-well library improved the phasing rate to 98%. Approximately 115,000 (approximately 4%) phased heterozygous SNPs were from the high-coverage LFR library and were not called in the standard library, indicating that MDA amplification and 160x coverage helped some regions get enough correctly called reads (5 or more). The high-coverage LFR phasing rate can be adjusted to balance haplotype completeness against phasing error.

Unit typing of European pedigrees . To further our understanding of the performance of LFR, we generated additional libraries from European ancestral pedigrees. CEPH family 1463 was selected because it has three generations of individuals, allowing for a comprehensive study of heritability. This family has previously been studied as part of a public data release (www.completegenomics.com/sequence-data/download-data/). Libraries were generated from individuals of each generation. A total of more than 1.6 Tb of sequence data was generated for NA12877, NA12885, NA12886, NA12891, and NA12892. In general, phasing was very high across all samples with approximately 92% of the attempted SNPs phased into the contig (Figure 20). Combining two LFR libraries (Figure 20) or LFR with parent-based phasing improved the overall ratio of phased SNPs to 97%. The N50 contig length across all analyzed family members was 500-600 kb. This length was limited to less than the length of NA19240. Investigation of the distribution of SNPs across the genomes of several different ethnic groups explains this difference.

The origin and impact of low heterozygosity regions in non-African populations . There are about twice as many 30kb-3Mb low heterozygosity regions (RLHs, defined as 30kb genomic regions with less than 1.4 heterozygous SNPs per 10kb, about 7 times lower than the planting density) in European pedigree samples than in NA19240, clarifying the previously reported relative excess of homozygotes in non-Africans (Gibson et al., Hum. Mol. Genet. 15: 789-795, 2006; Lohmueller et al., Nature 451: 994-997, 2008) and further supported by analysis of 52 complete genomes (Nicholas Schork, personal communication). These regions are obstacles to phasing, resulting in twice as small N50 contig lengths. More than 90% of contigs in European genomes end with these RLHs, which vary between unrelated individuals.

About 3% of all heterozygous SNPs in non-African genomes (30-60% of all unphased heterozygous SNPs) belong to these RLHs, which cover a very large fraction (30-40%) of these genomes. In Chinese and European genomes, longer RLHs cluster around 45 heterozygous SNPs per Mb (genomic coverage is about 1000 per Mb outside the RLHs), indicating that they shared a common ancestor around 37,000-43,000 years ago (based on a mutation rate of 60-70 SNPs per 20-year generation; Roach et al., Science 328:636-639, 2010; Conrad et al., Nat. Genet. 43:712-714, 2011). This may be due to a strong bottleneck at or after the time humans left Africa and within the previously determined range of 10,000-65,000 years ago (Li and Durbin, Nature 475:493-496, 2011). In addition, an excess of RLHs was observed on the X chromosome in European and Indian women (NA12885, NA12892, and NA20847) when compared to African women (NA19240), covering approximately 50% versus 17% of this chromosome, respectively (30% versus 14% for the entire genome in these same individuals). This indicates an even stronger out-of-Africa bottleneck on the X chromosome. A possible explanation is that substantially fewer women remained in Africa and had offspring with multiple men.

These observations suggest that whole-genome variation analysis in thousands of diverse genomes, including haplotype determination, will provide a deep understanding of human population genetics and the impact of these extensive "inbreeding" zones (which typically each contain more than 100 homozygous variants) on human disease and other extreme phenotypes. In addition, it shows that approximately 2,000 RLHs greater than 100 kb in length are present in all non-African individuals. Populations with a limited number of high-frequency haplotypes (which may be derived from recent bottlenecks or inbreeding (Gibson et al., Hum. Mol. Genet. 15: 789-795, 2006)) may also have long runs of the same heterozygous SNPs present in both parents, which limits the parents for phasing or allocating shorter LFR overlaps. Thus, population history and some reproductive methods can make phasing challenging, as shown by the X chromosomes of non-African women. Regardless of these factors, LFR phasing performance is roughly equivalent, phasing up to 97% of heterozygous SNPs in both European and African individuals, a result that should translate across all populations. In addition to combining LFR with standard genotyping of one parent as described below (a strategy that would be more limited to some families, as discussed above), using initial DNA fragments longer than 300 kb (e.g., by capturing cells in gel blocks or pre-purified DNA (Cook, EMBO J. 3: 1837-1842, 1984)) would span approximately 95% of all RLHs and determine haplotypes for most de novo mutations occurring in these regions. This would not be feasible with current fosmid cloning strategies (Kitzman et al., Nat. Biotechnol. 29: 59-63, 2011; Suk et al., Genome Res. 21: 1672-1685, 2011), which are limited to 40 kb fragments.

LFR reproducibility and phasing error rate analysis . In an effort to understand the reproducibility of LFR, we compared the haplotype data between the two NA19240 replicate libraries. In general, the libraries were very consistent, with the two libraries phasing only 64 differences in each library out of approximately 2.2 million heterozygous SNPs (Figure 22). This represents a phasing error rate of 0.003% or 1 error in 44Mb. LFR was also highly accurate when compared to the conservative but accurate whole chromosome phasing generated from the parental genomes NA19238 and NA19239, which had been previously sequenced by multiple methods. Only about 60 examples of 1.57 million equivalent single loci were found in which LFR phased variants that were inconsistent with the variants determined by the parental haplotype (if half of the inconsistencies were due to sequencing errors in the parental genomes, the phase rate was assumed to be 0.002%). The LFR data also contained approximately 135 overlapping groups (2.2%) per library that had one or more flipped haplotype blocks (Figure 22). Extending these analyses to the European descent repeat library of sample NA12877 ( FIG. 22 ) and comparing them to a recent high-quality family-based analysis (Roach et al., Am. J. Hum. Genet. 89:382-397, 2011) using 4 children of NA12877 and their mother NA12878 yielded similar results, assuming that each method contributes half of the observed inconsistencies. In both the NA19240 and NA12877 libraries, several contigs had numerous flipped segments. Most of these contigs tended to be located in regions of low heterozygosity (RLHs), low read coverage, or in repeat regions (e.g., subtelomeric or centromeric regions) observed in an unexpectedly large number of wells.

The unit type overlap group is classified into the parent chromosome . Most of the flip errors can be corrected by imposing the LFR phasing algorithm on the terminal overlap groups in these regions. Alternatively, these errors can be removed by simply and inexpensively adding the standard high-density array genotype data (about 1 million or more SNPs) from at least one parent to the LFR assembly. In addition, we found that the parental genotype can connect 98% of the LFR phased heterozygous SNPs between the whole chromosomes. In addition, this data allows the unit type to be classified into the maternal and paternal pedigrees, which can be used to incorporate parental imprint information in genetic diagnosis. If parental data are not available, population genotype data can also be used to connect the LFR overlap groups between the whole chromosomes, although this method can increase phasing errors (Browning and Browning, Nat. Rev. Genet. 12: 703-714, 2011). Even technically challenging approaches such as metaphase chromosome segregation (which has demonstrated whole chromosome haplotype determination) cannot assign parental origins without some form of parental genotype data (Fan et al., Nat. Biotechnol. 29:51-57, 2011). This combination of two simple techniques (i.e., LFR and parental genotyping) provides accurate, complete, and annotated haplotypes at a low cost.

Phased de novo mutations . As a demonstration of the completeness and accuracy of our diploid genome sequencing, we evaluated the phasing of 35 de novo mutations recently reported in the NA19240 genome (Conrad et al., Nat. Genet. 43: 712-714, 2011). 34 of these mutations responded in one of the standard genome or LFR libraries. Among those, 32 de novo mutations (16 from each parent) were phased in at least one of the two replicate LFR libraries. Not surprisingly, two non-phased variants reside in the RLH. Of these 32 variants, 21 were phased by Conrad et al. (supra), and 18 were consistent with the LFR phasing results. Three inconsistencies may be due to errors in previous studies (Matthew Hurles personal communication), confirming LFR accuracy without affecting the essential conclusions of the report.

Only LFR libraries were used to sequence the genome and determine the unit type from 100pg DNA . The analysis described above incorporates heterozygous SNPs from both standard and LFR libraries. However, in view of the complete presentation of the expected genome due to starting with the DNA amount equivalent to the DNA amount present in 10-20 cells, it is possible to use only LFR libraries. We have demonstrated that MDA provides sufficiently consistent amplification, and with high (80x) overall read coverage, the LFR library used alone allows detection of up to 93% of heterozygous SNPs without any modification to our standard library variation-response algorithm. In order to demonstrate the potential of using only LFR libraries, we phased NA19240 repeat 1 and an additional 250Gb of reads from the same library (500Gb in total). We observed that the total number of phased SNPs decreased by 15% and 5% (Figure 20), respectively. This result is not surprising in view of the fact that this library was generated from 60 pg DNA, instead of the optimal amount of 200 pg (Table 1 below) and also in view of the previously mentioned GC bias incorporated during in vitro amplification by MDA. Another 285 Gb LFR library responded and phased only 90% of all variants from the combined standard and LFR libraries (Figure 20). Despite the reduction in the total number of SNPs phased, the contig length was largely unaffected (N50>1 Mb).

Error reduction achieved by LFR for accurate genome sequencing from 10 cells . Substantial error correction (about 1 SNV in 100-1,000 kilobases of response) is a common attribute of all current massively parallel sequencing technologies. These rates may be too high for diagnostic purposes, and they complicate many studies searching for new mutations. The vast majority of false-positive variants are no longer likely to occur on maternal or paternal chromosomes. LFR can exploit this lack of consistent connectivity with surrounding true variants to eliminate these errors from the final assembled unit type. Both the Yoruba trio and the European pedigree provide an excellent platform for demonstrating the error reduction capabilities of LFR. We defined a set of heterozygous SNPs in NA19240 and NA12877 (greater than 85% of all heterozygous SNPs) that were reported with high confidence in each of the individual's parents as matching the human reference genome on both alleles. There are approximately 44,000 heterozygous SNPs in NA19240 and 30,000 in NA12877 that meet this criterion. The variants of the present invention are shown in Figure 23. The variants of the present invention are shown in Figure 24. The variants of the present invention are shown in Figure 25. The variants of the present invention are shown in Figure 25. The variants of the present invention are shown in Figure 25. The variants of the present invention are shown in Figure 25. The variants of the present invention are shown in Figure 25. The variants of the present invention are shown in Figure 25. The variants of the present invention are shown in Figure 25. The variants of the present invention are shown in Figure 25. Compared to genomes sequenced by standard methods, each LFR library exhibited a 15-fold increase in library-specific false positive responses before phasing. Most of these false-positive SNVs may have been introduced by MDA; sampling of rare cell line variants can cause a smaller percentage. Although LFR libraries were generated from 100 pg DNA and amplified via MDA, the application of the LFR phasing algorithm reduced the overall sequencing error rate to 99.99999% (approximately 600 false heterozygous SNVs/6 Gb), which is about 10 times lower than the error rate observed using the same ligation-based sequencing chemistry (Roach et al., Am. J. Human Genet. 89: 382-397, 2011).

Improve base calls with LFR information . In addition to phasing and eliminating false positive heterozygous SNVs, LFR can "rescue""no-response" positions or verify other responses (such as homozygous references or homozygous variants) by evaluating the origin of the holes that support the read results of each base call. As a demonstration, we found positions in the genome of NA19240 repeat 1 that had no response but were adjacent to adjacent phased heterozygous SNPs. In these examples, the positions were able to be "re-responded" because the phased heterozygous SNPs were indeed targeted at the presence of shared holes between the adjacent phased SNPs and the no-response position (Figure 24). Although LFR may not be able to rescue all no-response positions, this simple demonstration highlights the usefulness of LFR in more accurately responding to all genomic positions to reduce no-responses.

Highly divergent haplotypes exist in African and non-African genomes . Haplotype analysis, achieved through large-scale genotyping studies such as the HapMap project, is very important for understanding population genetics. However, the resolution of complete haplotypes for individuals is largely intractable or prohibitively expensive. Highly accurate haplotypes (filtering out clustered pseudoheterozygotes accumulated due to false positioning of repeat regions) (Li and Durbin, Nature 475:493-496, 2011; Roach et al., Science 328:636-639, 2010) will help understand many population phenomena found within individual genomes. As a demonstration, we scanned the LFR contigs of NA19240 for highly divergent regions between the maternal and paternal copies. 7,000 10-kb regions containing greater than 33 SNVs were identified; a 3-fold increase over the expected 10 SNVs. Assuming 0.1% standing variation and 0.15% base divergence per million years (based on 1% divergence between the human and chimpanzee genomes, which evolved approximately 6 million years from a common ancestor), our calculations suggest that approximately 50 Mb of these regions found in this African genome (approximately 2.0% of the "non-inbreeding" genome) may have evolved separately for more than 1.5 million years. If the chimpanzee-human split was less than 5 million years ago, this estimate is closer to 1 Myr (Hobolth et al., Genome Res. 21:349-356, 2011). This genome-wide analysis is consistent with a current study by Hammer et al. (Proc. Natl. Acad. Sci. USA 108:15123-15128, 2011) of several targeted genomic regions in African populations, assuming possible interbreeding between different hominin species in Africa. Our analysis revealed that 2.1% of non-inbred genomes of European descent also harbor similarly divergent sequences, often at different genomic locations. Most of these were likely introduced before humans left Africa.

A single genome contains multiple genes with inactivating variants in both alleles . Highly accurate diploid genomes are essential for making human genome sequencing valuable in clinical settings. To demonstrate how LFR can be used in a diagnostic/prognostic setting, we analyzed nonsense and splice site disrupting variants in the coding SNP data of NA19240. We further analyzed all missense variants using PolyPhen2 (Adzhubei et al., Nat. Methods 7: 248-249, 2010) to select only those variants that encode unfavorable changes. Both "likely to damage" and "probably to damage" are considered unfavorable for protein function because they are both nonsense mutations. 3485 variants matched these criteria. After phasing and removing false positives, only 1252 variants remained; that is, a significant reduction in potentially misleading information. We further reduced the list to examine only those 316 heterozygous variants, at least two of which co-occurred in the same gene. Using phased data, we were able to identify 189 variants present in the same allele within 79 genes. The remaining 127 SNPs were found to be dispersed among 47 genes with at least one adverse variation in each allele ( FIG. 25 ). This number was increased to 65 genes by combining two LFR libraries to determine the haplotypes for NA19240. Extending this analysis to European ancestry demonstrated that a similar number of genes (32-49 with coding mutations in both alleles) were potentially altered to the point of expressing little to no effective protein product ( FIG. 25 ). Extending this analysis to variants that disrupt transcription factor binding sites (TFBS) introduced an additional approximately 100 genes per individual. Many of these are likely to be partial loss of function changes or no loss of function changes. Due to the high accuracy of LFR, it is unlikely that these variants are the result of sequencing errors. Many of the adverse mutations found may have been introduced during the proliferation of these cell lines. A small number of these genes were found in unrelated individuals, suggesting that they may be incorrectly annotated or the result of systematic positioning or reference errors. The genome of NA19240 contains approximately 10 additional genes in the complete loss of function category; this is most likely due to a bias introduced by annotating the African genome using a European reference genome. However, these numbers are consistent with those found in several current studies of phased individual genomes (Suk et al., Genome Res. 21: 1672-1685, 2011; Lohmueller et al., Nature 451: 994-997, 2008) and suggest that most generally healthy individuals may have a small number of genes that are not absolutely required for normal life and encode ineffective protein products. We have demonstrated that LFR can place SNPs into haplotypes at large genomic distances where the phase of those SNPs can cause potential complete loss of function to occur. Such information will be crucial for effective clinical interpretation of patient genomes and for carrier screening.

TFBS destruction associated with allelic expression differences . Long unit types covering both cis-regulatory regions and coding sequences are crucial for understanding and predicting the expression level of each allele of a gene. By analyzing the 5.6Gb non-exhaustive expression data from RNA sequencing of lymphocytes from NA20431, we identified a small number of genes with significant differences in allelic expression. In each of these genes, a 5kb regulatory region upstream of the transcription start stable point and a 1kb scan downstream were found to have significantly altered binding sites for more than 300 different transcription factors (Sandelin et al., 32: D91-D94, 2004). In six examples (Figure 26), it was found that the 1-3 bases between the two alleles were different in each gene, causing a significant impact on one or more putative binding sites and potentially explaining the differential expression observed between the alleles. Although this is only one data set and it is not yet clear how much impact these changes have on transcription factor binding, these results demonstrate that with the type of large-scale studies (Rozowsky et al., Mol. Syst. Biol. 7:522, 2011) that become feasible using LFR unit typing, it is possible to elucidate the consequences of sequence changes on transcription factor binding sites.

Discussion . We have demonstrated the ability of LFR to accurately phase up to 97% of all detected heterozygous SNPs in the genome into long contiguous segments of DNA (N50 of 400-1500 kb in length). Even in the absence of candidate heterozygous SNPs from a standard library and thus using only 10-20 human cells, LFR libraries phased were able to phase 85-94% of the available SNPs, although current implementations have limitations. In several examples, the LFR libraries used in this article had less than optimal starting input DNA (e.g., NA20431). The improvements in phasing rates seen by combining two replicate libraries (samples NA19240 and NA12877) or starting with more DNA (NA12892) are consistent with this conclusion. Additionally, insufficient representation of GC-rich sequences resulted in fewer genomes being responded to (90-93% versus greater than 96% for the standard library). Improvements to the MDA method (e.g., by adding region-specific primers or using less amplification by improving yields in other steps) or the way we implement base and variant calling in LFR libraries (perhaps by using the assignment of reads to wells) could help improve coverage in these regions. Furthermore, as the cost of whole-genome sequencing continues to decline, higher-coverage libraries (which significantly improve call rates and phasing) will become more affordable.

The consensus haploid sequence is sufficient for many applications; however, it lacks two very important parts of the data for personalized genomes: phasing heterozygous variants and identification of false positive and negative variant responses. One of the goals of a personal genome is to detect disease-causing variants and to determine with extreme confidence whether an individual carries such variants or has one or two unaffected alleles. By independently providing sequence information from both maternal and paternal chromosomes, LFR is able to detect regions in the genome assembly that have only covered one allele. Similarly, false positive responses are avoided because LFR independently sequences both maternal and paternal chromosomes 10-20 times in different aliquots. The result is a statistically low probability that a random sequence error will recur in several aliquots at the same base position on a parental allele. In this way, LFR allows for the first time both accurate and cost-effective sequencing of genomes from a small number (preferably 10-20) of human cells, despite the use of in vitro DNA amplification and the resulting large number of unavoidable polymerase errors. In addition, by phasing SNPs at hundreds of kilobases to multiple megabases (or by integrating LFR with conventional genotyping of one or two parents in the entire chromosome), LFR can more accurately predict the impact of complex regulatory variants and parental imprinting on allele-specific gene expression and function in multiple tissue types. In short, this provides a highly accurate report on potential genomic changes that can cause protein function gain or loss. The clinical use of such information obtained cheaply for each patient is less than that of genomic data and will be crucial. In addition, the successful and affordable diploid sequencing of the human genome starting from 10 cells opens the possibility of comprehensive and accurate genetic screening from a variety of tissue sources, such as circulating tumor cells or microbiopsies of pre-implantation embryos generated via in vitro fertilization.

While many different forms of embodiments satisfy the present invention, as described in detail in conjunction with the preferred embodiments of the present invention, it should be understood that the present disclosure should be considered illustrative of the principles of the invention and is not intended to limit the invention to the specific embodiments illustrated and described herein. Many changes can be made by those skilled in the art without departing from the spirit of the invention. The scope of the invention will be measured by the appended claims and their equivalents. The abstract and the invention title should not be construed as limiting the scope of the invention, as their purpose is to enable the appropriate authorities and the general public to quickly ascertain the general nature of the invention. In the appended claims, unless the term "means" is used, none of the features or elements described therein should be construed as falling within the means-plus-function limitation of 35 U.S.C. §112,6.

The present invention provides the following:

1. A method for determining the sequence of a complex nucleic acid of one or more organisms, the method comprising:

(a) receiving, on one or more computing devices, a plurality of reads of the complex nucleic acid; and

(b) generating an assembled sequence of the complex nucleic acid from the reads using the one or more computing devices, the assembled sequence comprising less than one spurious single nucleotide variant per megabase at a call rate of 70% or greater.

2. The method of claim 1 , further comprising identifying a plurality of sequence variants in the assembled sequence, and phasing the plurality of sequence variants to produce a phased sequence.

3. The method of claim 2, comprising phasing at least three of the sequence variants, and identifying sequence variants that are inconsistent with the phasing of at least two sequence variants as errors.

4. The method of claim 2, wherein the assembled sequence is a whole genome sequence, and the method comprises phasing at least 70% of the sequence variants.

5. The method of claim 2, wherein the assembled sequence is a whole genome sequence, and the method comprises phasing at least 80% of the sequence variants.

6. The method of claim 2, wherein the assembled sequence is a whole genome sequence, and the method comprises phasing at least 85% of the sequence variants.

7. The method of claim 2, wherein the assembled sequence is a whole genome sequence, and the method comprises phasing at least 90% of the sequence variants.

8. The method of claim 2, wherein the assembled sequence is a whole genome sequence, and the method comprises phasing at least 95% of the sequence variants.

9. The method of claim 1 , wherein the step of receiving multiple read results for the complex nucleic acid comprises receiving multiple read results for each of a plurality of aliquots, each aliquot comprising one or more fragments of the complex nucleic acid.

10. The method of claim 9, comprising calling a base at a position in the assembled sequence based on preliminary base calls at that position from two or more aliquots.

11. The method of claim 9, comprising identifying as true a base call that occurs three or more times in the read results of two or more aliquots.

12. The method of claim 9, wherein an aliquot-specific tag is attached to each of said fragments, said method further comprising determining which aliquot gave said read result by identifying said aliquot-specific tag.

13. The method of claim 12, wherein the aliquot-specific tags comprise an error correction code, and each read comprises tag sequence data and fragment sequence data, wherein the tag sequence data is correct tag sequence data or incorrect tag sequence data comprising one or more errors; the method further comprising:

(c) correcting the incorrect tag sequence data using the error correction code, thereby generating corrected tag sequence data and uncorrectable tag sequence data;

(d) using the read comprising the corrected tag sequence data and the corrected tag sequence data in a first computer method requiring tag sequence data and generating a first output; and

(e) using the reads comprising the uncorrectable tag sequence data in a second computer method that does not require tag sequence data, and generating a second output.

14. The method of claim 13, wherein the first computer method is selected from the group consisting of: sample multiplexing, library multiplexing, phasing, and error correction methods using tag sequence data.

15. The method of claim 13, wherein the second computer method comprises positioning, assembly, and set-based statistics.

16. The method of claim 13, wherein the error correction code is a Reed-Solomon code.

17. The method of claim 1, further comprising:

(c) providing a first phased sequence of a region of the complex nucleic acid, the region comprising short tandem repeats;

(d) comparing the reading result of the first phased sequence for the region with the reading result of the second phased sequence for the region; and

(e) identifying, based on the comparing, an expansion of the short tandem repeat in one of the first phased sequence or the second phased sequence.

18. The method of claim 1 , further comprising obtaining genotypic data from at least one parent of the organism, and generating an assembled sequence of the complex nucleic acid from the read results and the genotypic data of the at least one parent.

19. The method of claim 1 , further comprising adding population genotype data and generating an assembled sequence of the complex nucleic acid from the read results and the population genotype data.

20. The method of claim 1, further comprising:

(c) aligning the plurality of reads of the first region of the complex nucleic acid, thereby creating overlaps between the aligned reads;

(d) identifying N heterozygous candidates within the overlap, where N is an integer greater than 2;

(e) clustering the space of ^2N to ^4N possibilities of the N heterozygous candidates or a selected subspace of the space, thereby creating a plurality of clusters;

(f) identifying two clusters with the highest density, each identified cluster containing a substantially noise-free center; and

(g) repeating steps (a)-(d) for one or more additional regions of the complex nucleic acid.

21. The method of claim 1 , wherein the assembled sequence comprises fewer than 0.8 false single nucleotide variants per megabase.

22. The method of claim 1 , wherein the assembled sequence comprises fewer than 0.6 false single nucleotide variants per megabase.

23. The method of claim 1 , wherein the assembled sequence comprises fewer than 0.4 false single nucleotide variants per megabase.

24. The method of claim 1 , wherein the assembled sequence comprises less than 0.2 false single nucleotide variants per megabase.

25. The method of claim 1 , wherein the assembled sequence comprises less than 0.1 false single nucleotide variants per megabase.

26. The method of claim 1 , wherein the assembled sequence has a response rate of at least 80% of the complex nucleic acid.

27. The method of claim 1 , wherein the assembled sequence has a response rate of at least 85%.

28. The method of claim 1 , wherein the assembled sequence has a response rate of at least 90%.

29. The method of claim 1, further comprising: (a) providing a quantity of the complex nucleic acid, and (b) sequencing the quantity of the complex nucleic acid to generate the plurality of reads.

30. The method of claim 1 , wherein the complex nucleic acid is selected from the group consisting of a genome, an exome, a transcriptome, a methylome, a mixture of genomes from different organisms, a mixture of genomes from different cell types of one organism, and a subset thereof.

31. The method of claim 1 , wherein the organism is a mammal.

32. The method of claim 1 , wherein the organism is a human.

33. One or more computer-readable non-transitory storage media storing the assembled human genome sequence produced by the method of item 1.

34. A computer-readable non-transitory storage medium storing instructions that, when executed by one or more computing devices, cause the one or more computing devices to implement the method of item 1.

35. A method for determining a human genome sequence, the method comprising:

(a) receiving, on one or more computing devices, a plurality of reads of the genome; and

(b) generating, using the one or more computing devices, an assembled sequence of the genome from the reads, the assembled sequence comprising fewer than 600 false single nucleotide variants per gigabase at a genome call rate of 70% or greater.

36. The method of claim 34, wherein the assembled sequence of the human genome comprises a genome call rate of 70% and an exome call rate of 70% or greater.

37. A computer-readable non-transitory storage medium storing instructions that, when executed by one or more computing devices, cause the one or more computing devices to implement the method of item 35.

38. A method for determining a human genome sequence, the method comprising:

(a) receiving, on one or more computing devices, a plurality of reads from each of a plurality of aliquots, each aliquot comprising a fragment of the human genome; and

(b) generating, with the one or more computing devices, a phased assembly of the genome from the reads, the assembly comprising fewer than 1000 false single nucleotide variants per gigabase at a genome call rate of 70% or greater.

39. A computer-readable non-transitory storage medium storing instructions that, when executed by one or more computing devices, cause the one or more computing devices to implement the method of item 38.

Sequence Listing

<110> Complete Genomics, Inc.

Drmanac, Radoje

Peters, Brock A.

Kermani, Bahram G.

<120> Processing and analysis of complex nucleic acid sequence data

<130> 92171-836153 (5039-US)

<140> US 13/448,279

<141> 2012-04-16

<150> US 61/546,516

<151> 2011-10-12

<150> US 61/527,428

<151> 2011-08-25

<150> US 61/517,196

<151> 2011-04-14

<160> 10

<170> PatentIn version 3.5

<210> 1

<211> 20

<212> DNA

<213> Artificial sequence

<220>

<223> Synthetic polynucleotides

<400> 1

ccgcagtagc ttacgaatcg 20

<210> 2

<211> 20

<212> DNA

<213> Artificial sequence

<220>

<223> Synthetic polynucleotides

<400> 2

gatttaactg agcacttggc 20

<210> 3

<211> 10

<212> DNA

<213> Artificial sequence

<220>

<223> Synthetic polynucleotides

<400> 3

aacgagtatt 10

<210> 4

<211> 10

<212> DNA

<213> Artificial sequence

<220>

<223> Synthetic polynucleotides

<400> 4

tttggcgttc 10

<210> 5

<211> 10

<212> DNA

<213> Artificial sequence

<220>

<223> Synthetic polynucleotides

<400> 5

gtagtaccgg 10

<210> 6

<211> 10

<212> DNA

<213> Artificial sequence

<220>

<223> Synthetic polynucleotides

<400> 6

aactgagcgg 10

<210> 7

<211> 12

<212> DNA

<213> Artificial sequence

<220>

<223> Synthetic polynucleotides

<400> 7

cagtcaagtg at 12

<210> 8

<211> 12

<212> DNA

<213> Artificial sequence

<220>

<223> Synthetic polynucleotides

<400> 8

catgatgagg ac 12

<210> 9

<211> 12

<212> DNA

<213> Artificial sequence

<220>

<223> Synthetic polynucleotides

<400> 9

tcttagcatg ta 12

<210> 10

<211> 12

<212> DNA

<213> Artificial sequence

<220>

<223> Synthetic polynucleotides

<400> 10

gtaactattc ag 12

Claims (25)

1. A method for analyzing the genomic DNA of an organism, the method comprising: Multiple reads corresponding to fragments of genomic DNA from multiple aliquots are received at one or more computing devices, each fragment of genomic DNA being tagged with an aliquot-specific tag sequence, and each read containing a sequence from the fragment of genomic DNA and an aliquot-specific tag sequence, wherein the genomic DNA contained in each aliquot of the multiple aliquots is less than a haploid genome equivalent. The origin of the read segment was determined by identifying the specific tag sequence of the aliquot sample; The following uses one or more computing devices to generate a phase sequence from the read segment: Identify multiple heterozygous loci corresponding to at least a portion of the genome of the organism; and Phased sequences of the plurality of heterozygous loci to generate a first haplotype and a second haplotype, wherein the phasing uses aliquots of reads corresponding to the plurality of heterozygous loci to determine which alleles at the heterozygous loci are in the same haplotype, the phasing sequence corresponding to at least a portion of the genome of the organism, wherein phasing of the plurality of heterozygous loci includes: For each pair of multiple heterozygous loci, A matrix is used to determine the number of shared sample segments between alleles of the heterozygous loci of the pair on the read, wherein the heterozygous loci of the pair are located within a specified distance from each other. 2. The method of claim 1, wherein phasing of the plurality of heterozygous loci further comprises: The score and orientation of the corresponding pair of heterozygous loci are calculated using each matrix; and The first and second unit types are determined using the scores and directions. 3. The method of claim 2, wherein the direction specifies which allele of the first heterozygous locus of the corresponding pair is connected to the first allele of the second heterozygous locus of the corresponding pair, and wherein the forward direction specifies that the two alleles are connected as in the list, and the reverse direction specifies that the two alleles are connected in the reverse order of the list. 4. The method of claim 3, wherein a score is calculated for the connection of corresponding pairs of heterozygous loci, and wherein the calculation comprises: Determine the first value in the positive direction; and A second value for the opposite direction is determined, wherein the direction is determined based on the larger of the first value and the second value. 5. The method of claim 2, wherein a score is calculated for the connection of corresponding pairs of heterozygous loci, and wherein the calculation comprises: The impurity value is determined as the ratio of the sum of all matrix elements except the two connected matrix elements to the total sum of all matrix elements; and The score is calculated using the impurity values and the two matrix elements. 6. The method of claim 5, wherein the score is determined using a fuzzy inference engine based on the impurity value and the two matrix elements corresponding to the connection. 7. The method of claim 2, wherein determining the first and second unit types using the score and the direction comprises: Based on the scores and directions, the connection graph between pairs of heterozygous loci is optimized. 8. The method of claim 7, wherein the graph is optimized by generating a minimum span tree. 9. The method of claim 7, wherein optimizing the connection graph provides a plurality of subtrees, the method further comprising: Each of the plurality of subtrees is simplified into an overlap group, thereby forming a plurality of overlap groups; and Sequencing information from the parent organism of the said organism is used to phase multiple contigs to generate the first haplotype and the second haplotype. 10. The method of claim 7, further comprising: The first heterozygous locus is removed from the graph as a node when the first heterozygous locus has no at least two connections to another heterozygous locus in one direction and no at least one connection to another heterozygous locus in another direction. 11. The method of claim 1, wherein the matrix for determining the number of shared subsamples between alleles of a specific pair at a heterozygous locus comprises locating the reads to the heterozygous locus of the specific pair and calculating the location reads of the subsamples sharing the alleles. 12. The method of claim 1, further comprising: The assembly sequence of the first and second unit types is generated using the one or more computing devices. 13. A method for analyzing the genomic DNA of an organism, the method comprising: Multiple reads corresponding to fragments of genomic DNA from multiple aliquots are received at one or more computing devices, each fragment of genomic DNA being tagged with an aliquot-specific tag sequence, and each read containing a sequence from the fragment of genomic DNA and an aliquot-specific tag sequence, wherein the genomic DNA contained in each aliquot of the multiple aliquots is less than a haploid genomic equivalent. The origin of the read segment was determined by identifying the specific tag sequence of the aliquot sample; One or more computing devices are used to generate multiple assembled sequences aligned with overlapping regions of the genome of the organism, each assembled sequence in the overlapping regions corresponding to a different aliquot sample; The following uses one or more computing devices to generate a phase sequence from the read segment: Identify a plurality of heterozygous loci corresponding to at least a portion of the genome of the organism, wherein the plurality of heterozygous loci comprise N heterozygous loci, where N is an integer greater than 1; and Phased sequences of the plurality of heterozygous loci to generate a first haplotype and a second haplotype, wherein the phasing uses aliquots of reads corresponding to the plurality of heterozygous loci to determine which alleles at the heterozygous loci are in the same haplotype, the phasing sequence corresponding to at least a portion of the genome of the organism, wherein phasing of the plurality of heterozygous loci includes: Based on the alleles at N heterozygous loci for the corresponding assembled sequence, multiple clusters are created by clustering the assembled sequences in a space of ^2N to ^4N possibilities; and The two clusters with the highest density were identified as corresponding to the first and second modular types. 14. The method of claim 13, wherein phasing of the plurality of heterozygous loci comprises: Calculate an N-dimensional matrix, where each dimension corresponds to a heterozygous locus, and each matrix element corresponds to the number of assembled sequences having allelic combinations corresponding to the matrix element; Identify the first matrix element and the second matrix element, each of which is the center of one of the two clusters; Determine the first haplotype at N heterozygous loci from the first matrix element; and The second haplotype is determined at N heterozygous loci from the second matrix element. 15. The method of claim 14, further comprising: Weights are assigned to each of the matrix elements, wherein a first weight of a first combination of alleles observed in the population of interest is greater than a second weight of a second combination of alleles observed in the population of interest. 16. The method of claim 13, further comprising: Phase morphology is repeated for multiple regions, each corresponding to a different plurality of heterozygous loci, thereby forming two contigs from the two clusters with the highest density for each of the plurality of regions, to obtain multiple contigs; and Sequencing information from the parents of the organism is used to phase the multiple contigs to generate the first and second haplotypes. 17. A method for analyzing the genomic DNA of an organism, the method comprising: Multiple reads corresponding to fragments of genomic DNA from multiple aliquots are received at one or more computing devices, each fragment of genomic DNA being tagged with an aliquot-specific tag sequence, and each read containing a sequence from the fragment of genomic DNA and an aliquot-specific tag sequence, wherein the genomic DNA contained in each aliquot of the multiple aliquots is less than a haploid genomic equivalent. The origin of the read segment was determined by identifying the specific tag sequence of the aliquot sample; The following uses one or more computing devices to generate a phase sequence from the read segment: Identify multiple heterozygous loci corresponding to at least a portion of the genome of the organism; and The plurality of heterozygous loci are phased to produce a first haplotype and a second haplotype. The phasing uses aliquots of reads corresponding to the origins of the plurality of heterozygous loci to determine which alleles at the heterozygous loci are in the same haplotype. The phased sequence corresponds to at least a portion of the genome of the organism. Identify phased SNPs from multiple heterozygous loci, wherein the phased SNPs have a first allele and a second allele; Identify loci that are the nearest neighbors of the phased SNP, wherein the loci have no response and the loci have reads with third and fourth alleles. Calculate the first number of shared aliquots containing the first allele at the phase-determined SNP and the third allele at the locus; and The third allele at the locus is determined based on the first number of the shared equally divided samples. 18. The method of claim 17, further comprising: When the first number of the shared equally divided samples is higher than a threshold, the location of the third allele at the locus is determined, and the threshold is 2 or higher. 19. The method of claim 17, further comprising: Calculate the second number of shared aliquots containing the second allele at the phased SNP and the third allele at the locus; Calculate the third number of shared aliquots containing the first allele at the phase-determined SNP and the fourth allele at the locus; and When the first number and the second number are higher than the threshold, and the third number is lower than the threshold, the locus is determined to be homozygous for the third allele. 20. The method of claim 17, further comprising: Calculate the second number of shared aliquots containing the second allele at the phased SNP and the fourth allele at the locus; When all reads carrying the third allele share an equal sample with the first allele, and all reads carrying the fourth allele share an equal sample with the second allele, the locus is determined to be heterozygous for both the third and fourth alleles. 21. A computer-readable storage medium storing instructions that, when executed by one or more computing devices, cause the one or more computing devices to perform the method of any one of claims 1-20. 22. A computer system comprising the computer-readable storage medium of claim 21. 23. A computer system comprising one or more processors configured to perform the method of any one of claims 1-20. 24. A computer system comprising means for carrying out the method of any one of claims 1-20. 25. A computer system, comprising: One or more computer-readable storage media for storing assembled whole human genome sequences having no more than one pseudomononucleotide variant per megabase and a response rate of at least 70%, wherein the assembled whole human genome sequence is generated by sequencing between 1 pg and 10 ng of human genomic DNA.