Reconstruction by tracking reads with variable errors

Abstract

Polynucleotide sequencing generates multiple reads of polynucleotide molecules. Many or all of the reads contain errors. Tracking reconstruction requires multiple reads generated by a polynucleotide sequencer and uses these multiple reads to accurately reconstruct the nucleotide sequence of the polynucleotide molecule. Some reads may contain errors that cannot be corrected. Therefore, there may be reads that can be used over their entire length and other reads with uncertain errors that cannot be corrected. When an uncertain error is discovered, the portion of the read with the error is skipped, and the sequence of the read after the error is used to reconstruct the tracking instead of discarding the entire read. The amount of skipped reads is determined by the location of the subsequence that matches the consensus sequence of other reads after the error. The analysis continues at the location determined by the matched location.

Description

Background Art

[0001] Much of the world's data is stored on magnetic and optical media today. Tape technology has recently seen significant density improvements with a single tape cartridge storing 185TB and is the densest form of storage commercially available today at approximately 10GB / mm 3 Recent studies have reported the feasibility of optical disks capable of storing 1PB, yielding a density of about 100GB / mm 3 Despite this improvement, storing zettabytes (2 70 Data stored in 1000 terabytes or 100 terabytes will still take up millions of units and use a lot of physical space. But storage density is only one aspect of storage media; durability is also important. Spinning platters have a lifespan of 3 to 5 years, and magnetic tapes have a lifespan of 10 to 30 years. Long-term archival storage requires data refreshes to replace failed units and refresh the technology.

[0002] The demand for data storage is growing exponentially, but the capacity of existing storage media has not kept up. Polymers of deoxyribonucleic acid (DNA) can store information at high density. The theoretical density limit is 1 exabyte / mm 3 (10 9 GB / mm 3 ). Less than 100 grams of DNA could store all of the man-made data in the world today. DNA is also very durable, with an observed half-life of over 500 years under certain storage conditions. DNA is therefore attractive as an information storage technology due to its high information density and long life. Yet another advantage of DNA as a storage medium is its continued relevance. Operating systems and standards for storage media will change, potentially rendering data on older storage systems inaccessible. But DNA-based storage has the benefit of everlasting relevance: as long as DNA-based life exists, there will be a strong reason to maintain technology that can read and manipulate DNA.

[0003] Despite the advantages of the DNA storage system, it must overcome several challenges. For example, DNA synthesis, degradation during storage, and sequencing are all potential sources of error. Therefore, the DNA sequence output by the sequencer may be different from the DNA sequence originally provided to the oligonucleotide synthesizer. Summary of the invention

[0004] This Summary is provided to introduce a selection of concepts in a simplified form that are further described below in the Detailed Description. This Summary is not intended to identify key features or essential features of the claimed subject matter, nor is it intended to be used to limit the scope of the claimed subject matter.

[0005] Binary data currently used by computers to store text files, audio files, video files, software, etc. can be represented as a series of nucleic acids (i.e., DNA or ribonucleic acid (RNA)) in polynucleotides. There are a variety of techniques for representing the 0 and 1 of binary data as a series of nucleotides. A variety of techniques are known to those of ordinary skill in the art. Polynucleotide sequences are designed to store binary data and then synthesized using an oligonucleotide synthesizer. The synthesized polynucleotides are placed in a storage device and eventually read by a polynucleotide sequencer. The data generated by the polynucleotide sequencer is decoded to recover the stored binary data. The machines that write and read polynucleotide sequences are not 100% accurate and introduce errors. Some types of errors (such as insertions, deletions, or substitutions of nucleotides) can be identified and corrected. Other types of errors (especially "burst" errors, where multiple errors are present in local "bursts" adjacent to or close to each other) may be difficult or impossible to correct. Therefore, with some decoding techniques, sequence reads including burst errors may not be available. The present disclosure provides techniques for extracting useful information from reads of polynucleotide sequences that include burst errors or other types of indeterminant errors.

[0006] The total length of the reads that can be generated by polynucleotide sequencing of a batch of polynucleotides is many times longer than the total length of all polynucleotides in the sequencing batch. This is called "coverage depth". Coverage depth greater than 1 (e.g., 10×, 20×, 30×, etc.) indicates that on average, the number of times each polynucleotide is provided to the sequencer is indicated by the coverage depth. Because this is an average value, some polynucleotides may not be sequenced at all, while other polynucleotides are sequenced many times more than the coverage level. For a single sequencing run, the coverage depth can be calculated by dividing the total number of base pairs sequenced by the total length of the polynucleotides provided to the sequencing machine. Each read in multiple sequencing reads of the same polynucleotide may have a slightly different sequence. Reads that are similar to each other and are therefore likely to be generated by sequencing the same polynucleotide can be grouped together in a cluster for further analysis. On average, the number of reads in a cluster is the same as the coverage depth.

[0007] Analysis of the reads in the cluster is performed to identify a single common output sequence from the multiple reads in the cluster. This type of analysis may be referred to as "trace reconstruction". The common output sequence is more likely to represent the actual nucleotide sequence in the source polynucleotide molecule than any single read. The technology for creating a common sequence from multiple reads is known to those skilled in the art. Some technologies identify errors in a single position of a read and correct the errors based on the sequences of other reads at the same position. Errors that can be identified and corrected include insertions, deletions, and substitutions. However, when multiple errors are close together, such as an insertion followed by a substitution and then two deletions, the current technology for creating a common output sequence cannot correctly identify the error, and the entire read may be ignored when determining the common output sequence. The technology discussed in the present disclosure includes identifying the position of an indeterminate error, skipping the portion of the read containing the indeterminate error, and determining the common output sequence using subsequent base calls in the read, rather than ignoring or discarding data from reads with indeterminate errors.

[0008] Uncertain errors can be identified at positions where the consensus sequence derived from other reads in the same cluster does not match in the read and the type of error cannot be determined. For example, if a given base call in a read is different from the base calls of all other reads at the same position, but the difference cannot be identified as an insertion, deletion, or substitution error, then the base call may be characterized as an uncertain error. A set number of positions can be skipped, and the search window down the read is evaluated to determine whether there is a matching subsequence between the read with the uncertain error and the other reads in the cluster. Each subsequence of a given length within the search window can be evaluated to determine whether it matches the sequence from other reads. The match can be determined in part by the match between the base call in the read with the uncertain error and the base call of the consensus sequence. Additionally, the match can be determined by comparing the base call of the read with the uncertain error with the simple majority vote of the other reads. The simple majority vote actually considers the most common base call at a given position, without considering errors such as insertions, deletions, and substitutions.

[0009] If a match is found, a single position within the matching subsequence is used as a candidate position, and subsequent positions adjacent to the candidate position are used as positions for comparison, to continue using the reads for determining a common output sequence. Therefore, the position before the indefinite error in the read can be used together with some or all of the other reads in the cluster to determine a portion of the common output sequence. The positions of the part where the read has the indefinite error and some number of positions after the indefinite error are ignored, and therefore, the common output sequence of these positions is determined by the other reads in the cluster. After the candidate position is located, the next position will be used again together with the other reads in the cluster to determine the common output sequence. Therefore, instead of ignoring all data including the indefinite error provided by the read, only the part near the indefinite error of the read is ignored, and the part before and after the error of the read can be used for determining the common output sequence.

[0010] The alignment of multiple reads in a cluster can be based on the first position in the read or the last position in the read. As the generation of the common output sequence proceeds position by position (e.g., "front" to "back" or "back" to "front"), errors and uncertainties may accumulate due to insertions and deletions that cause reads to be out of phase with respect to each other. Known sequences in all reads are intentionally inserted as alignment anchors during synthesis, and reference locations other than the ends are provided to re-align multiple reads. Therefore, multiple reads in a cluster can be aligned relative to each other based on alignment anchors and based on the first and last positions of the reads. Reads can be designed to include one or more alignment anchors. Because the alignment anchors include alternative alignment points in addition to the beginning and end of the reads, the reads may be split into shorter segments at the alignment anchors, and each segment in the shorter segments can be evaluated separately to determine the common output sequence. The common output sequence of each segment in the shorter segment can be re-joined to create a single common output sequence for the full length of the original read.

[0011] Specific implementations and examples in the present disclosure may refer only to DNA; however, it is to be understood that the present disclosure is equally applicable to any polynucleotides including DNA, RNA, DNA-RNA hybrids, and polynucleotides including synthetic or non-natural bases. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] The detailed description is set forth with reference to the accompanying drawings. In the drawings, the left-most digit(s) of a reference number identifies the drawing in which the reference number first appears. The use of the same reference numbers in different drawings indicates similar or identical items.

[0013] Figure 1 An illustrative architecture for tracking the operation of a reconstruction system is shown.

[0014] Figure 2 is an illustrative schematic diagram showing the use of a tracking reconstruction system.

[0015] Figure 3 An illustrative representation of substitution errors identified in accordance with the techniques of this disclosure is shown.

[0016] Figure 4 An illustrative representation of missing errors identified in accordance with the techniques of this disclosure is shown.

[0017] Figure 5 An illustrative representation of insertion errors identified in accordance with the techniques of this disclosure is shown.

[0018] Figure 6 An illustrative representation of adventitious error and inactivity tracking in accordance with techniques of this disclosure is shown.

[0019] Figure 7 An illustrative representation of introducing delay into an inactivity trace and locating candidate positions within the inactivity trace is shown.

[0020] Figure 8 The placement of alignment anchors within the reads is shown.

[0021] Fig. 9 Generation of partial consensus sequences by segmentation of reads divided by alignment anchors is shown.

[0022] Fig.10 A block diagram of an illustrative tracking and reconstruction system is shown.

[0023] Fig.11A and Fig. 11B An illustrative process for determining a consensus output sequence from multiple reads is shown.

[0024] Fig. 12A and Fig. 12B An illustrative process for generating binary data from reads received from a polynucleotide sequencer is shown.

[0025] Fig.13 An illustrative process for generating a consensus output sequence by omitting a portion of reads that contain burst errors is shown.

[0026] Fig.14 An illustrative process for analyzing reads containing adventitious errors in determining whether there are located subsequences beyond the adventitious errors that match other reads is shown. DETAILED DESCRIPTION

[0027] As mentioned above, DNA has great potential as a storage medium for digital information. However, dealing with errors that may damage the data is one of the challenges of using DNA to store digital data. There are many steps involved in converting digital data to synthetic DNA molecules and then recovering digital data from synthetic DNA molecules. The technology described in this disclosure provides a technology for using some information in a DNA read when the read includes an indefinite error.

[0028] The term "DNA strand" or simply "strand" refers to a DNA molecule. As used herein, "DNA reads", "sequence reads" or simply "reads" refer to data strings generated by a polynucleotide sequencer when the polynucleotide sequencer reads the sequence of a DNA strand. Because reads are data strings, they can also be referred to as "strings". However, the reads generated by the polynucleotide sequencer often contain errors and therefore cannot represent the structure of the DNA strand with 100% accuracy. Fortunately, most DNA sequencing technologies are capable of generating multiple reads of a DNA strand. Reads are called "noise reads" because each read may contain one or more errors, and the distribution of these errors is roughly random. A given read may also be error-free, but unless referenced to each other, it may not be known which reads are error-free and which contain errors. The technology disclosed herein uses multiple different noise reads for a single DNA strand to create a common output sequence that may represent the true sequence of the DNA strand. The common output sequence is a data string similar to any read, but the common output sequence is generated by analyzing the reads rather than being output directly from the polynucleotide sequencer. The process of going from many noisy reads to a roughly accurate consensus output sequence is called "challenge reconstruction."

[0029] A naturally occurring DNA chain is composed of four types of nucleotides: adenine (A), cytosine (C), guanine (G), and thymine (T). A DNA chain or polynucleotide is a linear sequence of these nucleotides. The two ends of a DNA chain (called the 5' and 3' ends) are chemically different. DNA sequences are conventionally represented starting from the 5' nucleotide end. The interactions between different chains are predictable based on the sequence: two single chains can bind to each other and form a double helix (if they are complementary): A in one chain is aligned with T in the other chain, and C and G are aligned similarly. The two chains in the double helix have opposite directionality (the 5' end is attached to the 3' end of the other chain), so the two sequences are "reverse complements" of each other. Two chains do not need to be completely complementary to be able to bind to each other. Ribonucleic acid (RNA) has a structure similar to DNA, and naturally occurring RNA is composed of four nucleotides A, C, G and uracil (U), instead of T. For the sake of brevity and readability, the discussion in this disclosure only refers to DNA, but RNA can be used instead of DNA or in combination with DNA.

[0030] The trace reconstruction problem can be stated using the following mathematical notation. Let ∑ denote a finite alphabet, e.g., ∑ = {A, C, G, T}. Let X∈∑ n is the sequence of interest, which can be arbitrary or random. The goal is to accurately reconstruct the sequence of DNA strand X from a collection of noisy reads.

[0031] Some sequencing technologies, such as sequencing by synthesis (discussed more below), have error distributions where the noise is at least to some extent independently distributed. Therefore, the noise can be modeled as being independently and identically distributed (iid) throughout the chain. This type of noise can be created using synthetic test data. To do so, let Y 1 , Y 2 , …, Y m Let p be an iid sequence obtained from X in the following way. d 、p i and p s represent the probability of deletion, insertion and substitution respectively, so that p = p d +p i +p s ∈[0, 1]. To obtain the noisy read segment Y, start from an empty string and perform the following operations for the comparison position j=1, 2..., n:

[0032] (no error), with probability 1-p, copy X[j] to the end of Y and increase j by 1;

[0033] (missing), with probability p d , increase j by 1;

[0034] (Insert), with probability p i , copy X[j] to the end of Y, add a random symbol to the end of Y, and increase j by 1;

[0035] (replacement), with probability p s , add a random symbol to the end of Y and increase j by 1.

[0036] Other sequencing technologies, such as nanopore sequencing (discussed more below), have different error distributions, where errors are not independently distributed but tend to be clustered, so that if there is one error, other errors are likely nearby. Nanopore sequencing has more errors than synthetic sequencing. About 12% compared to about 1-2%. By creating localized high-error areas, "bursts" of errors can be created in the string. In addition to single position errors, nanopore sequencing results may also have larger errors, such as transpositions of adjacent base blocks.

[0037] Identify a single consensus sequence from multiple noisy reads - track reconstruction - output estimated reads It is an estimate of the true sequence X of the DNA strand. This estimate is created from the number m noisy reads Y. Therefore, The goal is to reconstruct X exactly, that is, to minimize Minimize the error in noisy reads while considering or ignoring the error in noisy reads. The probability of being different from X can be achieved by using as much of the noise read segment Y.

[0038] Figure 1 An illustrative architecture 100 for realizing tracking reconstruction system 102 is shown. In short, the digital information intended for DNA molecule storage is converted into information representing a nucleotide string. The information representing a nucleotide string (that is, a letter string representing a nucleotide base sequence) is used as a DNA synthesis template, which indicates that an oligonucleotide synthesizer 104 synthesizes DNA molecule nucleotides by nucleotide chemistry. Artificially synthesized DNA allows the creation of a synthetic DNA molecule with an arbitrary series of bases, in which each monomer of the base is assembled together as a polymer of nucleotides. Oligonucleotide synthesizer 104 can be any oligonucleotide synthesizer using any generally recognized DNA synthesis technology. The term "oligonucleotide" used herein is defined as a molecule comprising two or more nucleotides.

[0039] The coupling efficiency of the synthesis process is the probability that a nucleotide will bind to an existing partial chain at each step of the process. Although the coupling efficiency of each step can be higher than 99%, this small error still causes the product yield to decrease exponentially with increasing length and limits the size of oligonucleotides that can currently be effectively synthesized to about 200 nucleotides. Therefore, the length of the stored DNA chain is about 100 to 200 base pairs (bp). This length will increase with the development of oligonucleotide synthesis technology.

[0040] The synthetic DNA produced by oligonucleotide synthesizer 104 can be transferred to DNA storage repository 106. There are many possible ways to construct DNA storage repository 106. Except by attaching the identification sequence to the structure of DNA chain at molecular level, DNA storage repository 106 can be constructed by physically separating the DNA chain into one or more DNA pools 108. Here, DNA pool 108 is shown as a flip-top tube, which represents a physical container for multiple DNA chains. When DNA is stored in a liquid solution, DNA chains are usually most susceptible to being operated by biotechnology. Therefore, DNA pool 108 can be implemented as a chamber filled with liquid, which is water in many implementations, and thousands, millions or more individual DNA molecules can be present in DNA pool 108.

[0041] In addition to being in liquid suspension, the DNA chains in the DNA storage library 106 can also be in a glassy (or vitreous) state, as a lyophilized product, as part of a salt, adsorbed on the surface of a nanoparticle or in another format. The structure of the DNA pool 108 can be implemented as any type of mechanical, biological or chemical arrangement that maintains a certain volume of liquid including DNA in a physical location. The storage device can also be in a non-liquid form, such as a solid bead or by encapsulation. For example, a single flat surface on which droplets exist is an implementation of the DNA pool 108, and even if it is not completely enclosed in a container, the droplets are partially maintained by the surface tension of the liquid. The DNA pool 108 may include single-stranded DNA (ssDNA), double-stranded DNA (dsDNA), single-stranded RNA (ssRNA), double-stranded RNA (dsRNA), DNA-RNA hybrid chains, or any combination including the use of non-natural bases.

[0042] The DNA chains deleted from the DNA repository 106 can be sequenced with a polynucleotide sequencer 110. In some implementations, the DNA chains can be prepared to be sequenced by amplification using polymerase chain reaction (PCR) to create a large number of DNA chains, which are identical copies of each other. The need for PCR amplification before sequencing may depend on the specific sequencing technology used. Although the level of PCR is much lower than current sequencing technology, PCR itself may be a source of error. At present, PCR technology generally introduces about one error per 10,000 bases. Therefore, on average, for every 100 reads of 100 bases, there will be an error as a result of PCR. The errors introduced by PCR are usually randomly distributed, so the tracking and reconstruction system is able to correct some errors caused by PCR.

[0043] As mentioned above, the polynucleotide sequencer 110 reads the order of the nucleotide bases in the DNA chain and generates one or more reads from the chain. The polynucleotide sequencer 110 uses a variety of techniques to interpret molecular information, and errors can be introduced into the data in a systematic and random manner. Errors can generally be classified as substitution errors, in which the true nucleotide is replaced by incorrect base calls (e.g., A is exchanged with G), insertions or deletions, in which random base calls are inserted (e.g., AGT becomes AGCT) or missing (e.g., AGTA becomes ATA). The transposition of a longer DNA region is another type of error. Each position in the read is a single base call determined by the polynucleotide sequencer 110 based on the properties sensed by the components of the polynucleotide sequencer 110. The various properties sensed by the polynucleotide sequencer 110 vary depending on the specific sequencing technology used. Base call represents which of the four nucleotide bases (A, G, C and T (or U)) in the DNA (or RNA) chain is present at a given position in the chain. Sometimes, base calls are wrong, and this is a source of error introduced by sequencing. Polynucleotide sequencing includes any method or technique used to generate base calls from DNA or RNA strands.

[0044] The sequencing technology that can be used is synthesis sequencing ( Sequencing). Synthesis sequencing is based on amplifying DNA on a solid surface using reentry PCR and anchor primers. The DNA is fragmented and adapters are added to the 5' and 3' ends of the fragments. The DNA fragments attached to the surface of the flow cell channel are extended and bridged amplified. The fragments become double-stranded and the double-stranded molecules are denatured. Multiple cycles of denaturing solid-phase amplification can then create millions of clusters of about 1,000 copies of single-stranded DNA molecules of the same template in each channel of the flow cell. Primers, DNA polymerase, and four fluorophore-labeled reversible terminator nucleotides are used to perform sequential sequencing. After nucleotide incorporation, a laser is used to excite the fluorophore, and an image is captured and the identification of the first base is recorded. The 3' terminator and fluorophore from each incorporated base are deleted, and the incorporation, detection, and identification steps are repeated.

[0045] Another example of a sequencing technique that can be used is nanopore sequencing. A nanopore is a small hole with a diameter of about 1 nanometer. Immersing the nanopore in a conductive fluid and applying an electric potential across the nanopore results in a tiny current due to the conduction of ions through the nanopore. The amount of current flowing through the nanopore is sensitive to the size of the nanopore. When a DNA molecule passes through the nanopore, each nucleotide on the DNA molecule blocks the nanopore to varying degrees. Therefore, as the DNA molecule passes through the nanopore, the change in the current passing through the nanopore represents the reading of the DNA sequence.

[0046] Other currently known sequencing technologies include single-molecule real-time (SMRT) sequencing from Pacific Biosciences. TM ) technology, Helicos true single molecule sequencing (tSMS), SOLiD available from Applied Biosystems TM technology, the use of chemically sensitive field effect transistor (chemFET) arrays and the use of electron microscopes to read the nucleotide base sequence.

[0047] All technologies used to sequence DNA are associated with some degree of error, and the type and frequency of errors vary by sequencing technology. For example, sequencing by synthesis creates errors in approximately 2% of base calls, and the errors tend to be independently and evenly distributed (iid). Most of these errors are substitution errors. Nanopore sequencing has a higher error rate of approximately 15% to 40%, and most errors caused by this sequencing technology are deletions. Errors in sequences generated by nanopore sequencing tend to appear in clusters, so if there is an error in one position, there is a high probability that an adjacent position will also have an error. Therefore, errors may be considered "bursty". The error distribution of a particular sequencing technology can describe the overall frequency of errors, the distribution characteristics of the errors, and the relative frequencies of various types of errors.

[0048] In some implementations, the polynucleotide sequencer 110 provides quality information indicating the confidence level of the accuracy of a given base call. The quality information may indicate that there is a high confidence level or a low confidence level in a particular base call. For example, the quality information may be expressed as a percentage of the accuracy of the base call, such as 80% confidence. Additionally, the quality information may be expressed as the confidence level that each of the four bases is the correct base call for a given position in the DNA chain. For example, the quality information may indicate an 80% confidence that the base call is T, an 18% confidence that the base call is A, a 1% confidence that the base call is G, and a 1% confidence that the base call is C. Therefore, the result of the base call will be T because the confidence that the nucleotide is the correct base call is higher than that of any other nucleotide. The quality information cannot identify the source of error, but only suggests which base calls are more or less likely to be accurate.

[0049] The polynucleotide sequencer 110 provides an output, a plurality of noisy reads (typically a plurality of DNA strands), in an electronic format to the tracking reconstruction system 102. The output may include quality information as metadata otherwise associated with the reads produced by the polynucleotide sequencer 110.

[0050] The tracking and reconstruction system 102 can be implemented as an integral part of the polynucleotide sequencer 110. The polynucleotide sequencer 110 can include an onboard computer that implements the tracking and reconstruction system 102. Alternatively, the tracking and reconstruction system 102 can be implemented as part of a separate computing device 112 that is directly connected to the polynucleotide sequencer 110 via a wired or wireless connection that does not span a network. For example, the computing device 112 can be a desktop computer or a notebook computer that is used to receive data from the polynucleotide sequencer 110 and / or control the polynucleotide sequencer 110. The wired connection can include one or more wires or cables that physically connect the computing device 112 to the polynucleotide sequencer 110. The wired connection can be created by a headphone cable, a telephone cable, a SCSI cable, a USB cable, an Ethernet cable, FireWire, etc. A wireless connection can be created by radio waves (e.g., any version of Bluetooth, ANT, Wi-Fi IEEE 802.11, etc.), infrared light, etc. The tracking and reconstruction system 102 may also be implemented as part of a cloud-based or network-based system using one or more servers 114 that communicate with the polynucleotide sequencer 110 via a network 116. The network 116 may be implemented as any type of communication network, such as a local area network, a wide area network, a mesh network, an ad hoc network, a peer-to-peer network, the Internet, a cable network, a telephone network, etc. Additionally, the tracking and reconstruction system 102 may be implemented in part by any combination of the polynucleotide sequencer 110, the computing device 112, and the server 114.

[0051] Figure 2 The use of the tracking reconstruction system 102 is shown as part of the process of decoding information stored in a synthetic DNA strand 200. The synthetic DNA strand 200 is a molecule having a specific sequence of nucleotide bases. Figure 1 As shown, the synthetic DNA chain 200 can be stored in the DNA pool 108. The synthetic DNA chain 200 can be present in the DNA pool 108 as a single-stranded molecule, or can be hybridized with a complementary ssDNA molecule to form a dsDNA. The polynucleotide sequencer 110 generates the output of multiple noise reads 202 from a single synthetic DNA chain 200. The number of reads is related to the coverage depth of sequencing. The greater the depth of sequencing, the more average sequences generated from the DNA chain. Due to sequencing errors, many reads in multiple reads of the same DNA chain may have different sequences. Each read in the read has a length (n), which is nine (corresponding to nine bases in the synthetic DNA chain 200) in this example. In actual sequencer data, the noise read can have an arbitrary length that is not all equal to each other. Deletion and insertion are a reason for the change in read length. For a given read, the length of the read can be expressed as n, but for all reads, n is not necessarily the same. In actual implementation, due to the current limitation on the maximum length of a DNA chain that can be artificially synthesized, the length of the read segment may be between 100 and 200. The positioning on the read segment can be referred to as a "position", such as from position one to position nine in this example. As used herein, "base" refers to the positioning of a given monomer in a DNA molecule, and "position" refers to the positioning along a data string such as a read segment. Therefore, assuming there are no errors, the third base in the synthetic DNA chain 200 corresponds to the third position in the read segment generated by the polynucleotide sequencer 110.

[0052] In this example, the number (m) of noise reads 202 provided to the tracking reconstruction system 102 is 5. However, any number can be used. In some implementations, the number of noise reads 202 provided to the tracking reconstruction system 102 can be 10, 20, or 100. The number of noise reads 202 provided to the tracking reconstruction system 102 can be less than the total number of reads generated by the polynucleotide sequencer 110. A subset of the total number of reads generated by the polynucleotide sequencer 110 can be randomly selected or heuristically analyzed by the tracking reconstruction system 102. In addition to random selection, other techniques can also be used to select which subset of reads is passed to the tracking reconstruction system 102. For example, quality information can be used to identify the m reads with the highest confidence in base calls from all reads generated by the polynucleotide sequencer 110. In some implementations, only reads of certain lengths are selected.

[0053] Tracking reconstruction system 102 analyzes noisy reads 202 according to the techniques of the present disclosure and generates a consensus output sequence 204. Consensus output sequence 204 represents the sequence of nucleotides in synthetic DNA strand 200 with less error than any of the individual noisy reads 202, and ideally no error.

[0054] The converter 206 converts the common output sequence 204 into binary data 208, thereby retrieving the digital information stored in the DNA storage library 106. The converter 206 may use additional error correction techniques to correct any errors that may remain in the common output sequence 204. Therefore, the tracking reconstruction system 102 does not have to correct all types of errors because there are other error correction techniques that can be used to recover the binary data 208.

[0055] Although the implementation discussed herein involves obtaining binary data 208 from reads of synthetic DNA strands 200, the tracking reconstruction system 102 operates equally well on reads of natural DNA strands. The output from the polynucleotide sequencer 110 is a plurality of noisy reads 202 for both synthetic DNA and natural DNA. Therefore, in implementations that do not involve the use of synthetic DNA to store binary data 208, the tracking reconstruction system 102 can be used to remove errors from reads generated by the polynucleotide sequencer 110.

[0056] Figure 3 Techniques for identifying substitution errors are shown. Reads can be aligned at the starting position or any other position (such as the alignment anchor discussed in more detail below). The starting position may correspond to the 5' end of the DNA chain from which the read was generated. In the figures of the present disclosure, the 5' end is oriented to the left. The position 300 of the comparison across the read is represented by a solid rectangular box. As each position in the read is analyzed in turn, the position 300 of the comparison can move from left to right along the read. Immediately following the position 300 of the comparison is a look-ahead window 302 represented by two dashed rectangular boxes. The look-ahead window 302 "looks forward" to the right of the position 300 of the comparison or toward the 3' end. That is, if the read is represented as Y j , and the position 300 of the comparison is denoted as p[j], then the look-ahead window 302 of length w is composed of Y j [p[j]+1], ...Y j [p[j]+w]. When the position 300 of comparison moves, the look-ahead window 302 can move along the read segment. In this example, the length of the look-ahead window 302 is two positions, but it can be longer, such as three, four or more.

[0057] The majority common base 304 is the most frequent base call at the position 300 of the comparison. Here, four of the five reads have G at this position, and one read has T. Because G is the most numerous base call, the majority common base 304 is G. In some implementations, the majority common base 304 can be determined by considering the quality information of the corresponding base calls at the position 300 of the comparison. Each base call at the position 300 of the comparison can be weighted based on the associated quality information. For example, if there is an 80% confidence that a given base call is G, then it can be counted as 0.8G when determining the majority common base 304, while a 30% confidence that a given base call is C will be counted as 0.3C when determining the majority common base 304. Therefore, when identifying the majority common base 304 for a given compared position 300, the confidence of a single base call can be considered. Additionally or alternatively, all base calls having quality information indicating that the confidence of the base call is less than a threshold level (e.g., 15%) may be omitted from determination of majority base call 304. If two or more different base calls are present with equal frequency, majority consensus base 304 may be randomly selected therefrom.

[0058] A read having a base call different from the majority consensus base 304 at the position 300 of the comparison is referred to as a variant read. k , base call Y k [p[k]] does not match the majority consensus base 304. In this example, the variant read is the third strand 308. When analyzed at a given compared position 300, in any read grouping, there may be zero, one, or more than one variant read.

[0059] The look-ahead window has 306 in total and is determined from the look-ahead window 302 in a manner similar to the majority of the common bases 304. The determination that the look-ahead window has 306 in total may also be affected by quality information. The look-ahead window has 306 in total and can be based on base calls weighted by their corresponding confidence levels and / or by omitting base calls with confidence levels below a threshold. The look-ahead window has 306 in total and is determined by considering the reads of the positions 300 that are not variant reads for comparison. Therefore, here, the look-ahead window 302 is shown as covering non-variant reads but not covering variant reads 308. In this example, the most common base call in the first position of the look-ahead window 302 is C, and the most common base call in the second position of the look-ahead window 302 is T. Therefore, the look-ahead window has 306 in total, which is a string of two positions of base calls: CT.

[0060] Next, the base calls in the look-ahead window 310 of the variant read (CT) are compared to the look-ahead window total 306 (CT). Because they match, the mismatch at the second position in the third read 308 is classified as a substitution. In mathematical notation, if Y k [p[k]+1], ...Y k [p[k]+w] matches the look-ahead window, then Y k The mismatch in is classified as a substitution. The majority consensus base 304 is used as the base call for that position in the consensus output sequence 204.

[0061] After the error type of the variant read is classified, the position 300 of the comparison is moved to continue analyzing the reads. For each read that is not a variant read, the position 300 of the comparison is moved one position to the right. In this example, these are the first read, the second read, the fourth read, and the fifth read. For the variant read whose error type is classified as a substitution (here, the third read 308), the position 300 of the comparison is also moved one position to the right. Therefore, if Figure 3 As shown in the lower portion of FIG. 3 , the compared position 300 is shifted one position to the right for all reads. The analysis is repeated at this new compared position 312, and in this iteration, the second read 314 is identified as a variant read.

[0062] Figure 4 Techniques for identifying deletion errors are shown. The compared position 400 is analyzed again to determine the most frequent base call at that position. In this example, three of the five strands have base call T, one strand has base call G, and one strand has base call C. Therefore, the most common base call is T, and the majority consensus base 402 for that position in the reads is T. The first strand 410 and the fourth strand are identified as variant reads.

[0063] The base calls in the look-ahead window 404 of the strand of non-variant reads are compared to determine the look-ahead window total 406. In this example, the values ​​of the two base calls in the look-ahead window 404 for the three non-variant reads are GA, GA, and TG. Therefore, the most common series of base calls is GA, and this becomes the look-ahead window total 406.

[0064] The value of the base call (AG) in the look-ahead window 408 of the first strand 410 is different from the look-ahead window total 406 (GA). Therefore, the error type responsible for the mismatch in the first strand 410 is not classified as a substitution. However, the base call in the compared position 400 and all positions except the final position of the look-ahead window 404 (GA) match the look-ahead window total 406 (GA). Therefore, the error type for this position of the first strand 410 is classified as a deletion. In this example, the length (w) of the look-ahead window 404 is 2, and therefore, all positions except the final position of the look-ahead window 404 are w-1 or the first base of the look-ahead window 404. If the length (w) of the look-ahead window 404 is 3, then when determining whether the error type in the variant read is a deletion, the first two bases (3-1) of the look-ahead window 404 will be considered. In mathematical notation, if Y k [p[k]], ...Y k [p[k]+w-1] matches the look-ahead window, then Y k Mismatches in are classified as missing.

[0065] After the error type of the variant reads is classified, the comparison position 400 is shifted one position to the right for each read in the reads that are not variant reads and for each read in the reads where the error is classified as a substitution. For the first strand 410 classified as a deletion, the comparison position 400 is not shifted. It remains at the same G located in the fifth position of the first strand 410. After realigning the strands after the differential shift to the new comparison position 412 for the first strand 410 (i.e., zero) and the other strands (i.e., 1), as shown in FIG. Figure 4 , the absence becomes apparent. This realignment, due to the shifting of the newly compared position 412 by different amounts based on the classification of the error type, keeps the strands in phase, further improving the analysis along the strands. After the newly compared position 412 is shifted (or not depending on the error type), the analysis can be repeated, here to identify a mismatch in the fifth strand 414.

[0066] Figure 5The technology for identifying insertion errors is shown. As discussed above, three possible error types are substitution, deletion and insertion. As identification substitution and deletion errors, the identification of insertion errors begins with analyzing base calls in the position 500 of comparison to determine most common bases 502, and analyzing base calls in the look-ahead window 504 to identify that the look-ahead window has 506. In this example, most of the total bases 502 are T. The 5th read 510 is a variant read, because it has A instead of T at the position 500 of comparison. The look-ahead window has 506 base calls that are GA. The base calls in the look-ahead window 508 of the 5th read 510 have 506 mismatches with the look-ahead window, so the error type is not classified as substitution. The first base call for the variant read (AT) in the base call in the position 500 of comparison and the look-ahead window 508 has 506 (GA) mismatches with the look-ahead window, so the error type is not classified as deletion.

[0067] However, the base calls in the look-ahead window 508 of the fifth read 510 match the base calls of the majority consensus 502 and all base calls except the last base call of the look-ahead window consensus 506 (i.e., position w-1) (i.e., both are TG). Therefore, the error is classified as an insertion of an A at position 5 of the fifth read 510. In mathematical notation, if Y k [p[k]+1] matches the majority common base 502, so Y k [p[k]+2], ...Y k [p[k]+2], ... matches the first w-1 coordinates of the look-ahead window, which has 506, so Y k The mismatch in is classified as an insertion. After the chain after the difference of the comparison position 500 is moved back to the new comparison position 512, as shown in FIG. Figure 5 As shown in the lower part of , the insertion error becomes obvious. The position 500 of the comparison is advanced by two positions for the read with insertion (here, the fifth read 510). For the other chains, the position 500 of the comparison is advanced by one position for the read of the non-variant read, by one position for the read with substitution, and by zero position for the read with deletion error.

[0068] Figures 3 to 5 The examples shown illustrate the analysis of only one error type, respectively. However, the techniques of the present disclosure are equally applicable to read groups with multiple errors in a compared position. There may also be multiple error types across a single compared position, e.g., in a total of 20 reads (m=20), three reads may have substitutions, one read may have deletions, and one read may have insertions.

[0069] Figure 6The situation where the error type cannot be identified is illustrated. A read may have a base call in the compared position that does not match the majority consensus base call, so this is a variant read, but the base call in the compared position and in the look-ahead window may not show any relationship to be classified as a substitution, deletion, or insertion. This is an error that cannot be classified according to the above techniques and is therefore referred to as an "indeterminate error" 604. Even if other techniques other than those described in the present disclosure are used to identify the error type, errors that fail to match the majority consensus base call that cannot be resolved as a specific type of error can be classified as indeterminate errors.

[0070] use Figures 3 to 5 The illustrated technique is unable to classify the base call at the location of the uncertain error 604 in the fifth read 602 as a substitution, insertion, or deletion. The compared position 600 immediately before the location of the uncertain error 604 is the last position before the uncertain error 604.

[0071] One way to handle reads with undetermined errors is to discard the reads from further processing. Thus, if a read has an error, and it cannot be resolved into a single error type, the read will be omitted from further analysis. If read 602 is discarded, it may be indicated as an "inactive track" because read 602 does not contribute any information to the track reconstruction. Once read 602 is identified as an inactive track, a consensus output sequence is generated from the remaining active tracks 606. Active tracks 606 may be all other reads from the same cluster, or may be smaller than all other reads if some are discarded due to also including undetermined errors, low confidence levels, or other reasons.

[0072] However, discarding reads completely also discards useful information that may be contained in the portion of the read that does not have an adventitious error. For some sequencing technologies, such as nanopore sequencing, discarding every read with an adventitious error can significantly reduce the number of clusters recovered, as this may even leave very few reads that cannot perform trace reconstruction.

[0073] Another way to handle ambiguous errors is to use a bias or tiebreaker to promote classification. The bias can be based on the error distribution of the polynucleotide sequencer used to generate the read. For example, if the frequency of the polynucleotide sequencer known to generate substitution errors is much higher than the frequency of deletion or insertion errors, all ambiguous errors can be classified as substitutions. If the error can be identified as one type in two possible error types, the relative frequency of those error types of the polynucleotide sequencer technology can be used to select between them. For example, if the error has been identified as a deletion or insertion (but not a substitution), and the polynucleotide sequencer has 80% substitution errors, 15% deletion errors, and 5% insertion errors, the error can be classified as a deletion error, because in this example, a deletion error is more likely to occur than an insertion error. However, this may reduce the reliability of the total output sequence, because base calls that may be incorrect may be included in the calculation of the total base call.

[0074] Additionally or alternatively, the quality information of each base call can be used to classify ambiguous errors.In one implementation, when error cannot be resolved as a single error type, the position of comparison and all base calls in the look-ahead window with quality information can be omitted from the majority of total bases and the look-ahead window total, and this quality information indicates that the base call confidence is less than the threshold level.Therefore, the total base call of the relevant position is determined on the most reliable base call of multiple reads.Ignoring the base call of low confidence may cause the above-mentioned technology to be able to resolve error as a single error type.However, if the confidence level of most base calls is very low, even if the confidence level is consistent, selecting specific base calls based on the confidence level may also not produce accurate results.

[0075] Figure 7 A technique is shown for skipping a portion of an inactive trace 700 that contains an indeterminate error to search for a match with an active trace 702 at a later location further along the inactive trace 700. The inactive trace 700 may be marked as an inactive trace due to determining an indeterminate error, such as Figure 6 The string of base calls represented by each read in the cluster can be represented as Y 1 , Y 2 , ..., Y m .. The value of m can vary based on the sequencing technology and coverage depth. For example, using nanopore sequencing, m may be about 10 to 100. The inactive trace 700 can be represented as Y j .

[0076] The consensus output sequence 704 may be generated from the active trace 702 and the portion of the inactive trace 700. The consensus output sequence 704 may be represented as the best estimate of the sequence X of the polynucleotide provided to the oligonucleotide sequencer. As previously discussed, the consensus output sequence 704 identifies a plurality of consensus base calls while taking into account insertions, deletions, and substitutions to maintain the alignment of the reads relative to each other.

[0077] Position 706 is the position in the inactivity track 700 that is estimated when the last read segment became inactive. This can be Figure 6 The position of the comparison is shown as 600. This position can be represented as i * Then p * [j] is Y j The position in the inactive track 700 is the position of comparison when the last read became inactive. A delay 708 of several positions (e.g., 5 to 10) is introduced in the inactive track 700 before continuing the evaluation of the base calls in the inactive track 700. The delay 708 is part of the inactive track 700 and it will be maintained as inactive regardless of whether there are any potential matches in the region. The last "good" position in the inactive track 700 is identified at position 706, but it is unknown how many positions are affected by adventitious errors. In sequences that include bursty errors (such as sequences from nanopore sequencing), it is likely that there are several positions of mismatches that are adjacent or close to each other that cannot be classified as insertions, deletions, or substitutions. The introduction of delay 708 reduces the likelihood that the inactive track 700 will return to an active state at a position where a mismatch exists.

[0078] Figure 7 The schematic representation of the inactive trace 700 at the bottom of shows an illustrative arrangement of how the inactive trace 700 can be evaluated to identify matching sequences. When an indefinite error is identified, some portion of the inactive trace 700 may have been over-analyzed and used to develop a common output sequence 704. This portion of the inactive trace 700 is indicated as the analysis sequence 710. This is followed by a delay 708, and then past the delay 708 is a portion of the inactive trace 700 that has a sequence 712 to be analyzed. Once it is identified how far past the delay 708, the inactive trace 700 can be reactivated and trace reconstruction can continue.

[0079] After the delay 708 there is a location 714 for comparison, which may be denoted as i. To determine whether the location 714 for comparison is within the delay 708 region, i is compared to i. * For comparison: If ii * ≤ delay, then the inactive trace 700 is kept inactive, but if ii *>delay, then the inactivity trace 700 is evaluated at position i as follows.

[0080] For the alignment of the read segments that have not been calculated for the common output sequence 704, the common look-ahead sequence 724 can be created from the activity tracking 702. By performing a majority common vote on the activity tracking 702 at the position of the comparison, the common look-ahead sequence 724 can be created in the same manner as the look-ahead window described earlier. This technology only uses the majority of common base calls as base calls for the common sequence at each position. If there are equal numbers of base calls (e.g., 5 Ts and 5 Cs), the connection will be arbitrarily broken. However, unlike the look-ahead window, the common look-ahead sequence 724 is not limited to a window of two or three positions.

[0081] Once the delay 708 region is passed, a candidate position 716 is sought. Candidate position 716 is a single position in the inactive trace 700 that can be compared with the active trace 702 and used to continue determining the common output sequence 704. Candidate position 716 may be denoted as k. Candidate position 716 may be immediately following the delay 708 at the compared position 714, or may be located further after the end of the delay 708.

[0082] If insertions and deletions occur with equal frequency from position 706 to candidate position 716, then the position at which the inactive trace 700 can be aligned again with the active trace 702 will be at or near the comparison position 714. The comparison position 714 can be denoted as k 0 . This can be done by combining candidate position i and position i * The difference between * It is calculated by adding the position compared to when the last read became inactive. Therefore, k 0 :=p * [j]+(ii * ). However, the number of insertions and deletions from position 706 to candidate position 716 may not be equal. Therefore, search window 718 is used to evaluate multiple possible locations to achieve the best result in k 0 Matches between subsequences of the inactive trace 700 and the active trace 702 in the surrounding region.

[0083] Search window 718 may include 0 Previous backward search window (sbw), exceeding k 0 The forward search window (sfw) and k 0 Therefore, the search window 718 is adjacent to k 0 The possible position range of k. The multiple possibilities of k can be expressed as k∈{k 0 -sbw,k 0-sbw+1,...,k 0 -1, k 0 , k 0 +1, ..., k 0 +sfw-1,k 0 +sfw}. Each possibility of k can be evaluated to determine whether there is a subsequence of the inactive trace 700 that matches the active trace 702. The larger search window 718 resulting from larger values ​​of sbw and sfw makes it more likely to find a match, but also increases the likelihood of returning a false positive result. The size of the search window 718 can be determined experimentally. The total length of the search window 718 is sbw+sfw+1. Therefore, if sbw and sfw are both 5 positions long, then the length of the search window 718 will be 11 positions. sbw and sfw can be the same length as each other or different lengths.

[0084] The subsequence of the inactive trace 700 within the search window 718 is compared with the comparison sequence Z to determine if there is a match. The comparison sequence Z includes The base call at position i in , and may include one or both of a backward match (mb) 720 region and a forward match (mf) 722 region.

[0085] The backmatch 720 region is located before position i and is aligned with a portion of the common output sequence 704. The length of the backmatch 720 sequence can be determined experimentally and can range, for example, from 0 to 10 positions. Therefore, the backmatch 720 sequence can be omitted. The backmatch 720 sequence can also overlap with the delay 708. At the position represented by the delay 708, the common output sequence 704 is determined by the active tracking 702, and the base calls from the inactive tracking 700 are not used at the location of the delay 708.

[0086] The first mb+1 position of Z is Therefore, this represents the portion of the comparison sequence Z that is identical to the common output sequence 704 at the position being evaluated. Even though the length of mb is 0, Also included in Z.

[0087] The forward match 722 sequence extends past position i, so there is no consensus output sequence 704 at the corresponding position. Therefore, for the purpose of comparison with the sequence in the inactive trace 700, the consensus lookahead sequence 724 sequence from the active trace 702 is used. The length of the forward match 722 sequence can be, for example, 0 to 10 positions. Therefore, if mf=0, then the forward match 722 sequence is not used. Although the backward match 720 sequence and the forward match 722 sequence are compared to the consensus sequence generated from the active trace 702, each is compared to a different consensus sequence. The backward match 720 sequence is compared to the consensus sequence generated from the active trace 702. The consensus output sequence 704 is compared, while the forward match 722 sequence is compared to the consensus look-ahead sequence 724 .

[0088] The length of the subsequence of the inactive trace 700 compared with Z is mb+mf+1, which is the same length as Z. Therefore, the match is performed between two strings of equal length. Recall that the inactive trace 700 can be represented as Y j , therefore, the subsequence of the inactive trace 700 compared with Z can be represented as Because its positioning is based on the positioning of k. Specifically, can be defined as Y j [k–mb], Y j [k–mb+1], …, Y j [k-1], Y j [k], Y j [k+1], …, Y j [k+mf–1], Y j [k+mf].

[0089] Using Z and These definitions can be used for every k∈{k 0 -sbw,k 0 -sbw+1,...,k 0 -1, k 0 , k 0 +1, ..., k 0 +sfw-1,k 0 +sfw} to determine if there are any candidate positions 716 that represent a match. A match may not exceed a distance threshold (dt) in the edit distance between the two strings. If dt=0, then only exact matches will not be considered a match. The number of comparisons may be the same as the length of the search window 718. Thus, if the search window is 12 positions long, a sliding window of 12 positions is moved along the search window 718 and compared to the corresponding base call sequences in Z.

[0090] If none of the candidate positions 716 within the search window 718 is associated with a position that matches the corresponding Z subsequence subsequence, then the inactive track 700 can be kept as inactive. The entire read can be discarded, or only a portion of the read before the uncertain error (ie, the analysis sequence 710) can be used.

[0091] If there is a sequence associated with a Z subsequence that matches the corresponding Z subsequence within the search window If the candidate position 716 of the subsequence is a subsequence, the position of comparison for the trace reconstruction can be set to the position immediately after the candidate position 716. The location of the compared position can be calculated as p[j]:=k+1. The inactive trace 700 can be designated as an active trace again, and the trace reconstruction including the read segment can proceed from the compared position forward.

[0092] There may be multiple candidate positions 716 (multiple locations of k) at which a match is found. If so, one of the candidate positions 716 is selected to be used as the location for restarting the tracking reconstruction. One approach is to select the one closest to k. 0 Alternatively, the distance k 0 The farthest position may be chosen, or the selection may be random.

[0093] A given read may have more than one undetermined error, and the process of identifying undetermined errors and then finding candidate positions 716 to restart the tracking reconstruction may be repeated multiple times in a single read. In addition, a read cluster that may include dozens of sequences may have multiple reads that include undetermined errors. Therefore, the technique can be used for more than one read in a cluster (including all reads). Because at multiple locations, i different reads in multiple reads may be active or inactive, the set of reads that constitute the active tracking 702 may be different at different alignment locations.

[0094] Figure 8 is a schematic diagram of several DNA strands 800, 802, 804 indicating the location of alignment anchors. Alignment anchors provide a reference location in addition to the start and end of a read, which can align multiple noisy reads, thereby assisting in constructing a consensus output sequence.

[0095] Typically, the accuracy of the total output sequence is more accurate toward the beginning or end of the read, where there is a higher confidence in the correct alignment with each other in the chain. Without being bound by theory, any error may be accumulated and cause other errors when analyzing subsequent positions along the read. For example, if the missing error is mistakenly identified as a substitution error, the remaining base calls in the read may be out of phase, and have a negative impact on the accuracy of subsequent total base calls. One way to minimize this impact is to use the first half of the "forward" analysis and the first half of the "reverse" analysis. For example, it is assumed that the length of the read to be analyzed is 200 positions. Analysis from left to right can be performed, which provides a total output sequence of positions 1 to 100. Analysis from right to left can also be performed, which provides a total output of positions 101 to 200. Two analyses can be performed in parallel. The resulting total output sequence provides base calls for all positions 1 to 200 by combining the analysis from left to right and the analysis from right to left. This is referred to as a combined total output sequence.

[0096] Because the use of the combined consensus output sequence is least reliable towards the end of the read, one arrangement is to add an alignment anchor 806 at the center of the DNA strand 800. This provides a location to re-align the reads at the position where the alignment is least reliable relative to each other. The marker provides an "anchor" for the consensus alignment to reset and start over without being affected by previously accumulated errors.

[0097] Comparison anchor is used as the mark of positioning, and even in the presence of noise and error created by sequencing, it should be identifiable in each read. Comparison anchor may be a relatively short sequence of, for example, 4 to 8 positions. The selection of the base sequence in the comparison anchor and the length of the comparison anchor can be arbitrary. Any sequence consistent on a plurality of DNA chains clustered together can be used as a comparison anchor. Because the comparison anchor uses some number of bases, this is at the expense of reducing the number of bases in the DNA chain that can be used to store digital information. However, including the comparison anchor in the DNA chain may make it possible to be used for tracking reconstruction to generate an accurate common output sequence with a lower coverage level, thereby reducing the total amount of sequencing required for recovering the information stored in the DNA.

[0098] In implementation, the sequence of the alignment anchor is created to avoid repeated bases and avoid internal alignments. Because the alignment anchor is used to align multiple reads, it may be beneficial to reduce the possibility of misalignment between the alignment anchor sequences in two different reads. For example, ACGT is a sequence without any repeated bases and without any internal alignments. However, if the first ACG in one read is aligned with the last ACG in a different read, ACGACG may form a partial alignment with itself. For example, if there is an error in the base call of a position within the alignment anchor, or the read sequence adjacent to the alignment anchor is the same as a portion of the alignment anchor, misalignment may occur.

[0099] One or more alignment anchors may be included in the polynucleotide sequence of the DNA strand. There are a variety of possible arrangements of alignment anchors within the DNA strand. One arrangement is to place alignment anchor 806 centered in the middle of DNA strand 800. Placing alignment anchor 806 in the center of DNA strand 800 will of course result in alignment anchor sequences being located in the middle or near the corresponding reads. For example, alignment anchor 806 may be located between 40% to 60%, 45% to 55%, 49% to 51%, or 50% of the distance from the beginning of DNA strand 800 to the end of DNA strand 800.

[0100] DNA strands 802 and 804 show examples of multiple alignment anchors. If multiple alignment anchors are used, they can be arranged in various ways. DNA strand 802 illustrates an example of alignment anchors 808, 810, and 812 being equidistant. Therefore, the number of positions between the start of DNA strand 802 and alignment anchor 808, between alignment anchor 808 and alignment anchor 810, between alignment anchor 810 and alignment anchor 812, and between alignment anchor 812 and the DNA strand end are all substantially the same. In order to align multiple reads, the alignment anchor can serve as an additional "end" of DNA strand 802. Therefore, if DNA strand 802 is 200bp long, instead of creating a common output sequence from right to left for 100 positions and a common output sequence from left to right for 100 positions, alignment anchors 806, 810, and 812 will divide DNA strand 802 into four regions, each region being approximately 50bp. Therefore, the sequencing length before the alignment point is reached is shortened from 200 bp to approximately 50 bp, and the chance of erroneous base calls that place one of the reads out of phase with the other reads is also reduced.

[0101] In the different arrangements illustrated by DNA strand 804, the spacing of the alignment anchors can be biased toward the center of the DNA strand. Therefore, alignment anchors 814, 816, and 818 near the center of DNA strand 804 are positioned closer together than alignment anchors 820 and 822 far from the center. Moreover, alignment anchors 820 and 822 can be farther from the corresponding ends of DNA strand 804. Each alignment marker used together should be unique so that it can be distinguished from other markers. Therefore, anchors 814 to 822 will have different sequences respectively. The lengths of the individual anchors can also be different.

[0102] Fig. 9 A schematic diagram of the reads is used to illustrate how partial consensus sequences generated from portions of the reads partitioned by alignment anchors can be joined together to create a complete consensus output sequence. Fig. 9 The illustrated read segment 900 may be obtained, for example, from Figure 8 The alignment anchor 902 may be located in the middle of the read segment 900, approximately 50% of the way between the left and right ends of the read segment 900. Even if the alignment anchor 902 is located exactly 50% along the length of the DNA strand, insertions and / or deletions introduced during sequencing may shift the position of the alignment anchor 902 in the read segment so that it is no longer exactly in the middle. Therefore, the alignment anchor 902 divides the read segment 900 into two parts of approximately equal length.

[0103] Because the sequence of the alignment anchor 902 is known (assuming there are no errors in that portion of the read 900), a search for that known sequence (e.g., AGCT) can identify one or more matches in the entire read 900. The placement of the alignment anchor 902 is also known, so any matches that are not within a threshold distance from the middle of the read 900 can be ignored. Thus, the location of the alignment anchor 902 can be identified in the read 900 and other reads in the same cluster.

[0104] If read 900 is 180 positions long and alignment anchor 902 is four positions long and located exactly in the middle, the left subsequence of read 900 located to the left of alignment anchor 902 will be 88 positions long, and the right subsequence located to the right of alignment anchor 902 will also be 88 positions long. The left subsequence can be analyzed by forward and reverse analysis to generate a left-to-right consensus output sequence 904 and a right-to-left consensus output sequence 906. Similarly, the right subsequence can be analyzed to generate consensus output sequences 908 and 910.

[0105] As discussed above, only the first half of each sequence in these four total output sequences 904 to 910 can be used to minimize the error introduced by the cumulative phase shift caused by the insertion and deletion caused by order-checking. In this implementation, before switching to the total output sequence running in the opposite direction, only 44 positions will be analyzed. Alternatively, more than half of each sequence in the corresponding total output sequence can be used so that there is a certain degree of overlap near the center of the subsequence. In the overlapping region, the base call from the total output sequence from left to right and from right to left can be used to determine the total output base call that will be assigned to a certain position.

[0106] The combination of the left-to-right consensus output sequence 904 and the right-to-left consensus output sequence 906 can generate a partial consensus sequence 912 for the left half of the read 900. Similarly, the consensus output sequences 908 and 910 can be combined in any suitable manner to generate a partial consensus sequence 914 for the right half of the read 900. The partial consensus sequence 912 can be aligned to the partial consensus sequence 914 by reference to the alignment anchor 902.

[0107] Once aligned, partial consensus sequence 912 and partial consensus sequence 914 can be joined using alignment anchors 902 to create a single joined consensus sequence 916. Joined consensus sequence 916 represents the consensus output sequence for read 900, which was developed from forward and reverse consensus output sequences 904 to 910. This is in contrast to implementations that do not use alignment anchors, where the consensus output sequence for read 900 is created from only one forward and one reverse consensus output sequence.

[0108] Because alignment anchor 902 does not encode digital information like other parts of read 900, the base call corresponding to alignment anchor 902 can be deleted to create a consensus output sequence without alignment anchor 918. The consensus output sequence without alignment anchor 918 represents a base call string, which can be decoded to recover the digital information contained in read 900.

[0109] Fig.10 Shows Figure 1 1000 of an illustrative block diagram of a tracking and reconstruction system 102 is shown. To review, the tracking and reconstruction system 102 can be implemented in whole or in part in any one of the computing device 112, the polynucleotide sequencer 110, and the server 114. Thus, the tracking and reconstruction system 102 can be implemented in a system comprising (multiple) one or more processing units 1002 and a memory 1004, both of which can be distributed across one or more physical or logical locations. The (multiple) processing units 1002 can include any combination of a central processing unit (CPU), a graphics processing unit (GPU), a single-core processor, a multi-core processor, an application-specific integrated circuit (ASIC), a programmable circuit such as a field programmable gate array (FPGA), etc. In addition to hardware implementations, one or more of the (multiple) processing units 1002 can be implemented in software and / or firmware. The software or firmware implementation of the (multiple) processing unit 1002 can include computer or machine executable instructions written in any suitable programming language to perform the various functions described. The software implementation of the (multiple) processing unit 1002 can be stored in whole or in part in the memory 1004.

[0110] Alternatively or additionally, the functionality of the tracking and reconstruction system 102 may be at least partially performed by one or more hardware logic components. For example, and not limitation, illustrative types of hardware logic components that may be used include field programmable gate arrays (FPGAs), application specific integrated circuits (ASICs), application specific standard products (ASSPs), systems on chips (SOCs), complex programmable logic devices (CPLDs), etc. Implementation as hardware logic components may be particularly suitable for portions of the tracking and reconstruction system 102 that are included as onboard portions of the polynucleotide sequencer 110.

[0111] The tracking reconstruction system 102 may include one or more input / output devices 1006 , such as a keyboard, pointing device, touch screen, microphone, camera, display, speaker, printer, and the like.

[0112] The memory 1004 of the tracking reconstruction system 102 may include removable storage devices, non-removable storage devices, local storage devices, and / or remote storage devices to provide storage for computer-readable instructions, data structures, program modules, and other data. The memory 1004 may be implemented as a computer-readable medium. Computer-readable media include at least two types of media: computer-readable storage media and communication media. Computer-readable storage media include volatile and non-volatile media, removable and non-removable media implemented in any method or technology for storing information (such as computer-readable instructions, data structures, program modules, or other data). Computer-readable storage media include but are not limited to RAM, ROM, EEPROM, flash memory or other memory technology, CD-ROM, digital versatile disk (DVD) or other optical storage devices, magnetic cassettes, magnetic tapes, magnetic disk storage devices or other magnetic storage devices, or any other non-transmission media that can be used to store information for access by a computing device.

[0113] In contrast, communication media may embody computer readable instructions, data structures, program modules or other data in a modulated data signal, such as a carrier wave or other transmission mechanism. As defined herein, computer readable storage media and communication media are mutually exclusive.

[0114] The tracking reconstruction system 102 can be connected to one or more polynucleotide sequencers 110 through a direct connection and / or a network connection via a sequence data interface 1008. The direct connection can be implemented as a wired connection, a wireless connection, or both. The network connection can traverse a network 116. The sequence data interface 1008 receives one or more reads from the polynucleotide sequencer 110.

[0115] The trace reconstruction system 102 includes a plurality of modules that may be implemented as instructions stored in the memory 1004 for execution by the processing unit(s) 1002 and / or implemented in whole or in part by one or more hardware logic components or firmware.

[0116] The randomization module 1010 randomizes the input digital data before encoding it into the DNA with the oligonucleotide synthesizer 104. The randomization module 1010 can create a random, more accurate pseudo-random string from the input digital data by performing an XOR of the input string and the random string. The random string can be generated using a seeded pseudo-random generator based on a function and a seed. This randomization of the input digital data increases the randomness of the synthetic DNA chain 200, resulting in the noise reads 202 from the polynucleotide sequencer 110 itself having A, G, C, and T pseudo-random sequences. Randomness supports decoding (i.e., clustering and tracking reconstruction).

[0117] Clustering module 1012 clusters subsets of multiple reads based on the possibility that the subsets of multiple reads are derived from the same DNA chain. The data received from the polynucleotide sequencer 110 at the sequence data interface 1008 may include a set of reads generated from each DNA chain provided to the polynucleotide sequencer 110. Although there may be errors in many reads, the reads from the same DNA chain are generally more similar to each other than the reads from another DNA chain. If the set of reads to be analyzed includes reads of different DNA chains, further analysis will be hindered. Therefore, clustering can be performed so that the data for further analysis is limited to a subset of reads that are considered to represent the same DNA chain. Poorly formed clusters may be "poorly" formed due to including too much or not enough. A cluster that includes too much poor formation is a cluster that groups reads of more than one chain in a single cluster together. A cluster that includes insufficient poor formation is a cluster that should be grouped into a single large cluster but is divided into multiple smaller clusters. The clustering module 1012 can use any suitable clustering technique, such as connectivity-based clustering (e.g., hierarchical clustering), centroid-based clustering (e.g., k-means clustering), distribution-based clustering (e.g., Gaussian mixture models), density-based clustering (e.g., density-based spatial clustering with noise applications (DBSCAN)), etc. The tracking reconstruction system 102 can analyze one or more, including all clusters derived from the data output by the polynucleotide sequencer 110. As mentioned above, the number of reads in a cluster may depend on the depth of sequencing. Using Nanopore sequencing, each cluster may have approximately 10 to 100 reads.

[0118] The read alignment module 1014 aligns multiple reads at positions that span the comparison of multiple reads. Initially, the left end of the read (corresponding to the 5' end of the DNA strand) can be aligned. This first position in the read can be used as the position of the initial comparison. As the analysis proceeds, the read alignment module 1014 moves the position of the comparison multiple positions along each read based on the identified error type.

[0119] The read alignment module 1014 advances the compared position by one position for the reads with the majority of common base calls at the compared position, advances the variant read by one position if the error type is classified as a substitution, advances the variant read by zero positions if the error type is classified as a deletion, and advances the variant read by two positions if the error type is classified as an insertion.

[0120] The read alignment module 1014 can also generate a "reverse" alignment that starts with the reads aligned at the right end (corresponding to the 3' end of the DNA strand). The analysis is then performed in the same manner, except that the movement to the "right" is changed to a movement to the left. The common output sequence 204 may be different for the same set of reads when analyzed from left to right compared to from right to left.

[0121] Typically, the accuracy of the total output sequence 204 is more accurate towards the beginning of the read, where higher confidence is correct in the comparison. Without being bound by theory, any error may be accumulated and cause other errors when analyzing subsequent positions along the read. For example, if the missing error is mistakenly identified as a substitution error, the remaining base calls in the read may be out of phase, and have a negative impact on the accuracy of subsequent error identification. A way to minimize this impact is to use the first half of the "forward" analysis and the first half of the "reverse" analysis. For example, it is assumed that the length of the read to be analyzed is 100 positions. Analysis from left to right can be performed, and this analysis provides a total output sequence 204 for the first 50 positions. Analysis from right to left can be performed, and this analysis provides a total output sequence 204 for the last 50 positions. Two analyses can be performed in parallel. The total output sequence 204 of gained is a combination of 50 base pairs identified by analyzing from left to right and 50 base pairs identified by analyzing from right to left. This is called a combined total output sequence.

[0122] The alignment anchor module 1016 functions by identifying alignment locations at locations other than the ends of the reads when alignment anchor sequences are present. When alignment anchors are present, the alignment anchor module 1016 can function with the read alignment module 1014 to create segmented alignments of the reads between one of the ends of the reads and the alignment anchor and between the two alignment anchors (if present). The alignment anchor module 1016 can identify alignment anchor sequences present in the reads. Alignment anchor sequences are predetermined sequences, such as Figure 8 and Fig. 9 As shown. The alignment anchor sequence can be a sequence of 3 to 8bp. In some implementations, the alignment anchor sequence can be 4 to 6bp long. The alignment anchor sequence can be any sequence included by the design in the DNA chain synthesized by the oligonucleotide synthesizer 104. Therefore, some nucleotide bases in the DNA chain are not used to encode digital information, but are used to insert a known sequence-alignment anchor sequence. In order to reduce the probability of not being aligned between two alignment anchor sequences, the sequence used as the alignment anchor sequence can be designed without repeating nucleotides so that it is not aligned with itself. For example, the sequence "AGCT" can be a suitable sequence for the alignment anchor, while the sequences "ACACAC" and "TTGG" will be more likely to be unaligned.

[0123] Once the alignment anchor is identified and the read is aligned at the alignment anchor, the alignment anchor module 1016 can also divide the read into multiple sub-reads, such as a first sub-read and a second sub-read. If there is only one alignment anchor sequence in the read, the first sub-read can extend from the beginning of the read (e.g., the 5' end) to the alignment anchor sequence, and the second sub-read can extend from the alignment anchor sequence to the end of the read (e.g., the 3' end). The alignment anchor module 1016 can do this for each read in the read being analyzed, thereby creating a third sub-read, a fourth sub-read, and so on.

[0124] The read alignment module 1014 can then align the sub-reads in a manner similar to aligning the entire read, except that the alignment is based not only on the ends of the reads but also on the alignment anchors that divide the reads. Thus, the read alignment module 1014 can align the first subsequence from the first half of the first read with the third subsequence from the first half of the second read. The second half of the first read represented by the second subsequence can also be aligned with the fourth subsequence representing the second half of the second read. Thus, instead of aligning the first read and the second read in their entirety, subsequences of the reads are aligned. Of course, in many implementations, more than two reads will be aligned.

[0125] Because subsequences are shorter than whole reads, the effect of accumulated errors in the generation of the consensus output sequence is reduced. There are fewer positions where errors may cause one read to be out of phase with the other reads. If errors do accumulate, or the read is out of phase with other reads in the same cluster, the alignment will be "reset" at the alignment anchor sequence. Subreads can be aligned in both forward and reverse orientations, as described above for whole reads and as described above for Fig. 9 Therefore, each alignment of subsequences (e.g., the first subsequence, the second subsequence, the third subsequence, and the fourth subsequence from the above example) can be analyzed in the same manner as the alignment of the entire sequence.

[0126] The variant read identification module 1018 determines the majority of common base calls at the compared position, and marks the reads with different base calls at the compared position as variant reads. In some implementations, the variant read identification module 1018 can use the error distribution associated with the polynucleotide sequencer 110 discussed above to determine the majority of common base calls. The variant read identification module 1018 can mark or otherwise identify each read that is a variant read for a given compared position. The marking then identifies the read as a read that should be identified for determining the error type. With each movement of the compared position, the identification of which reads are variant reads and which are not variant reads changes.

[0127] The error classification module 1020 will be used for the error type classification of variant reads as substitution, deletion or insertion. If the error cannot be uniquely classified, the error classification module 1020 can indicate that the type of error is uncertain, or the error type is limited to a possibility in two possibilities. The classification of the error by the error classification module 1020 can be based in part on the comparison of the common string of base calls in the look-ahead window of the subset of multiple reads with the majority of common base calls (that is, non-variant reads) and the base calls in the variant reads at the position of comparison. When determining the error type, different from the variant reads analyzed in a given iteration are not used. If the error classification module 1020 classifies the error in the read as an indefinite error, the variant read identification module 1018 can mark the read as inactive tracking.

[0128] As above Figure 3 The error type is classified as a substitution based on the consensus string of base calls in the look-ahead window that matches the base call string in the variant read following the position of comparison as described in . Figure 4 The error type is classified as a deletion based on a look-ahead window that matches the compared position in the variant read and base calls at one or more of the following positions in common as described in .

[0129] As above Figure 5 As described in , the error type is classified as an insertion. It is classified as an insertion because (1) the base call in the variant read after the compared position (i.e., the first base call in the look-ahead window of the variant read) matches the majority consensus base call, and (2) the consensus string of base calls in the look-ahead window that matches the base call string in the variant read is equal in length to the look-ahead window and starts two positions after the compared position. As described above in Figure 6 As described in , if the error type cannot be classified as a substitution, deletion, or insertion, it is classified as an indeterminate error.

[0130] The consensus output sequence generator 1022 determines the consensus output sequence 204 based at least in part on the plurality of consensus bases and the error types identified along the aligned reads. Each position in the consensus output sequence 204 is a plurality of consensus bases at that position in view of the adjusted alignment of the reads due to the error type classification. The error distribution of the polynucleotide sequencer 110 and / or the quality information of each base call may also be used to determine the consensus output sequence 204 by affecting the identification of the plurality of consensus bases and error types.

[0131] The consensus output sequence generator 1022 may also generate and assemble a consensus output sequence for alignment of subsequences defined by alignment anchor sequences. Assembling a consensus output sequence for two or more subsequences of a read cluster may include generating Fig. 9 Forward consensus sequence and reverse consensus sequence of each subsequence shown. Forward consensus sequence and reverse consensus output sequence can be combined as described above to create multiple combined consensus sequences for each read cluster. (Multiple) sequences of alignment anchors are used to attach the combined consensus sequence together. For example, the first combined consensus sequence representing the first half of the read and the second combined consensus sequence representing the second half of the read can be attached to each other by the alignment at the alignment anchor sequence. If there are multiple alignment anchors, multiple combined consensus sequences can be attached in the correct order by the alignment at the corresponding alignment anchor sequence. After two or more combined consensus sequences are attached together, (multiple) alignment anchor sequences can be deleted to leave a string of combined consensus sequences as reads without alignment anchor sequences. The combined consensus sequence assembled by multiple subsequences may be more accurate than the combined consensus sequence manufactured from the entire read length without alignment anchors, because the number of positions traversed by any forward consensus sequence or reverse consensus sequence (and the number of read opportunities that are out of phase with each other) is less. Reducing the number of cumulative errors can also reduce the coverage or sequencing depth required to form a cluster that can produce an accurate consensus sequence.

[0132] The error correction module 1026 may apply additional error correction techniques to decode the common output sequence 204. In some implementations, the error correction module 1026 uses a non-binary error correction code to decode the common output sequence 204 based on the redundant data encoded as a chain. An example of this type of error correction is Reed-Solomon error correction. In an example implementation, the Reed-Solomon outer code may be added to the starting binary data and eventually distributed over many DNA chains (e.g., 10,000-100,000 chains) when stored. If a threshold number of errors is exceeded, it is possible that the Reed-Solomon error correction may not be able to decode the common output sequence 204. If this occurs, the tracking reconstruction may be repeated with one of the parameters changed. Changing the parameters may result in a different common output sequence 204 that the Reed-Solomon error correction can decode. The length of the look-ahead window (w) is a parameter that can be changed. A look-ahead window of length 3 may be used instead of a look-ahead window of length 2 (or vice versa). The cutoff threshold for marking a read as inactive tracking, the length of the delay, the error type classification for accepting base calls based on quality information, and biasing ambiguous errors can be changed by making the threshold more lenient or more stringent. After changing one or more parameters, it can be determined whether the consensus output sequence 204 is different from the previous consensus output sequence 204, and if so, Reed-Solomon error correction can be applied to the new consensus output sequence 204 to see if it can decode the sequence.

[0133] The conversion module 1028 converts the consensus output sequence into binary data 208 representing at least a portion of the digital file. The conversion from a series of base calls to a binary data string 208 is performed by reversing the operations originally used to encode the binary data 208 into a series of base calls. These operations are known to the entity operating the DNA storage library 106. In some implementations, Figure 2 The converter 206 introduced in the embodiment of the present invention may include the same functionality as the conversion module 1028 and the error correction module 1026 and possibly other modules. The binary data 208 may be used in the same manner as any other type of binary data. If the various error correction techniques are sufficient, the binary data 208 will represent a perfect reproduction of the original binary data.

[0134] Illustrative Process

[0135] For ease of understanding, the processes discussed in the present disclosure are delimited to separate operations represented as independent blocks. However, these separately delimited operations should not be interpreted as necessarily being order-dependent in their performance. The order in which the processes are described is not intended to be interpreted as limiting, and any number of the described process blocks can be combined in any order to implement the process or an alternative process. Moreover, one or more of the operations provided may also be modified or omitted.

[0136] Fig.11A and Fig. 11B A process 1100 for correcting insertion, deletion, and substitution errors in sequence data generated by a polynucleotide sequencer is shown. The process 1100 may be performed by Figure 1 , Figure 2 and Fig.10 The tracking reconstruction system 102 is shown to implement it.

[0137] In 1102, binary data to be encoded into one or more DNA strands is reversibly randomized. This randomization occurs before DNA strand synthesis and, in some implementations, is present in all reads. The reads may be randomized via Fig.10 The read segment can be received by the sequence data interface 1008 shown. Fig.10 The randomization module 1010 is shown for randomization.

[0138] In 1104, sequence data generated by the polynucleotide sequencer is clustered using a clustering technique. Any suitable clustering technique may be used, and one of ordinary skill in the art will be able to identify a suitable clustering technique. Clustering creates groups of reads derived from the same source DNA strand. Clustering may be performed by Fig.10The clustering module 1012 shown is executed. In some implementations, clustering is performed on randomized data to improve the ability of clustering techniques to accurately separate multiple reads into different groups. A poorly formed cluster is a cluster containing reads derived from different DNA chains. Techniques such as discarding reads that deviate from the consensus sequence by more than a threshold amount can prevent poorly formed clusters from affecting the final consensus output sequence. However, discarding the entire read will lose any useful information that may be obtained from the read.

[0139] In 1106, the multiple reads classified as representing DNA strands are received for further analysis. Based on the clustering performed in 1104, the reads can be classified as representing the same DNA strand. Multiple reads can also be classified as representing the same DNA strand due to the use of sequencing technology in which the input to the polynucleotide sequencer is only a single DNA strand (or substantially identical copies produced by PCR). In some implementations, multiple reads can be classified as representing the same DNA strand via Fig.10 The sequence data interface 1008 is shown as being received. In other implementations, the plurality of reads may be received after clustering performed by the clustering module 1012.

[0140] At 1108, the positions of the comparison across the plurality of reads are identified. The positions of the comparison may be similar to Figure 3 The comparison position 300 is shown. Figure 4 The comparison position 400 is shown. Figure 5 The comparison position shown is 500 or Figure 6 The comparison position 600 is shown. In implementation, the comparison position can be determined by Fig.10 The read alignment module 1014 is shown as an identifier.

[0141] At 1110, the majority consensus base call at the compared position is determined by identifying the most common base call at the position. Ties may be arbitrarily broken. As described above, the most common base call may be identified in part by considering quality information of the base calls present at the compared position. In implementation, the majority consensus base call may be determined by Fig.10 The variant read identification module 1018 is shown to determine.

[0142] In 1112, it is determined whether the base call at the compared position is the same as the majority consensus base call. If the same, then the analyzed read has the expected base call at the position, is not a variant read relative to the position, and the process 1100 follows the "yes" path to 1114.

[0143] At 1114 , the process 1100 advances along the read to the next position. However, if the base call at the compared position does not match the majority consensus base call, the process 1100 follows a “no” path from 1112 to 1116 .

[0144] At 1116, reads from the plurality of reads whose base calls in the compared positions differ from the majority consensus base calls are identified as variant reads. In implementations, this identification may be by Fig.10 The variant read identification module 1018 is shown to execute.

[0145] Move to Fig. 11B , in 1118, the consensus string of base calls in the look-ahead window adjacent to the compared position is compared with the base calls in the variant reads. The consensus string of base calls in the look-ahead window can be limited to base calls from a subset of reads that have a majority consensus base call at the compared position. For example, in a set of 10 or 20 reads, there may be more than one variant read because the base call at the compared position does not match the majority consensus base call. When the comparison is made for one of the variant reads, the base calls in the look-ahead window of the other variant reads are not considered. This is because the other variant reads may have deletion or insertion errors, which will cause the base calls in the look-ahead window to be out of phase and possibly incorrect. Based at least in part on the comparison, the type of error can be determined for the variant read. In implementation, the comparison can be performed by Fig.10 The error classification module 1020 is shown to perform. In one implementation, the length of the look-ahead window may be 2 to 4 positions.

[0146] At 1120, based on the consensus string of base calls in the look-ahead window being identical to the string of base calls in the look-ahead window after the position of the comparison of the variant read, the error type of the variant read is determined to be a substitution. Thus, the look-ahead window of the variant read matches the look-ahead window of the non-variant read. This relationship is, for example, Figure 3 is shown in .

[0147] In 1122, the position of the comparison of the variant read is advanced by one position.

[0148] In an alternative, in 1124, based on the consensus string of base calls in the look-ahead window being the same as the base call string in the variant read that includes the base call in the compared position and the adjacent base call, the error type of the variant read is determined to be a deletion. The length of this base call string in the variant read is equal in length to the length of the look-ahead window. Thus, for example, if the look-ahead window is three positions long, the base call string in the variant read includes the base call in the compared position and the next two base calls. This relationship is, for example, Figure 4is shown in .

[0149] The compared position of the variant read is advanced by zero positions at 1126. Not advancing the compared position of the variant read will realign the strands because of the presence of the deletion so that the strands will be in phase for subsequent analysis.

[0150] As yet another alternative, at 1128, the error type for the variant read is determined to be an insertion based on base calls that match two specific patterns. First, the base calls in the variant read after the compared position are identical to the majority consensus base calls, and second, the consensus string of base calls in the look-ahead window is identical to the base call string in the variant read sequence. The base call string in the variant read sequence is equal in length to the look-ahead window, and the starting position of the base call string is two positions after the compared position. This relationship is, for example, Figure 5 is shown in .

[0151] In 1130, the position of the comparison of the variant read is advanced by two positions. The position of the comparison is advanced by one position to account for the insertion, and is advanced by a second position because the position of the comparison is advanced by one position for all non-variant chains. This maintains the alignment between the chains for subsequent analysis.

[0152] In each of 1120, 1124, and 1128, the error type of the variant read at the position of interest can be determined based at least in part on the error distribution associated with the polynucleotide sequencer. Considering the error distribution of the polynucleotide sequencer can change one or both of the determination of the majority consensus base call and the consensus string of base calls in the look-ahead window. In implementation, consideration of the error distribution can be performed by the consensus output sequence generator 1024.

[0153] At 1132, it is determined whether the variant read is less than a threshold reliability level. The threshold level may be the number of errors in the variant read; the number of errors in the variant that cannot be uniquely classified; the minimum, median, or mode of confidence levels of base calls in the variant reads; or (multiple) other factors. The threshold number may be a number of positions in the variant read from one position to the total number of positions. The threshold number may also be a percentage from 1% to 100%. If the variant read is less than the threshold reliability level, the process 1100 proceeds to 1134 along the "yes" path.

[0154] At 1134, the variant reads are omitted, and a single consensus output sequence from the plurality of reads is used without the variant reads. After omitting the variant reads, process 1100 proceeds to 1136. Alternatively, if the variant read is not less than a threshold reliability level (i.e., is considered reliable), the variant read is used for further analysis, and process 1100 proceeds along the "no" path to 1136. If the variant read cannot be classified as an insertion, deletion, or substitution due to the presence of an adventitious error, the variant read is classified as such, and rather than discarding the read, it can be further analyzed, as follows Fig.13 and / or Fig.14 shown.

[0155] At 1136, it is determined whether there are additional unanalyzed positions in the read. Thus, it is determined whether the "end" of the read has been reached and the majority of consensus base calls have been identified for the positions of the read. If used to break the read into multiple smaller segments, the end of the read can also be identified by the alignment anchor. If the analysis has not yet reached the end, then the process 1100 proceeds along the "yes" path to 1138.

[0156] In 1138, the positions of the subset of the plurality of reads having the majority of consensus base calls at the compared positions (i.e., the non-variant reads) are advanced by one. The newly compared positions of the variant reads are advanced by a certain amount based on the identified error type in 1122, 1126, or 1130. The newly compared positions may be similar to the newly compared positions 312, 412, 512, and 608, such as Figures 3 to 6 Process 1100 then returns to 1108 and the analysis continues.

[0157] At 1136 , if there are no unanalyzed positions in the read segment, process 1100 proceeds along the “no” path to 1140 .

[0158] At 1140, a single consensus output sequence is determined based in part on the majority consensus base calls and the error types. The single consensus output sequence may be determined by Fig.10 The common output sequence generator 1024 is shown to determine.

[0159] Fig. 12A and Fig. 12B A process 1200 is shown for recovering binary data encoded in a synthetic DNA strand by accounting for insertion, deletion, and / or substitution errors. The process 1200 may be performed by Figure 1 , Figure 2 and Fig.10 The tracking reconstruction system 102 is shown to implement it.

[0160] In 1202, the binary data to be encoded as DNA is reversibly randomized by taking the exclusive OR (XOR) of the binary data and a random sequence generated by a seed and a function. This operation affects DNA strands that, when read, produce reads that also have randomized properties. In implementation, randomization may be performed by Fig.10 The randomization module 1010 is shown to perform.

[0161] At 1204, a plurality of reads are received from a polynucleotide sequencer. In an implementation, the plurality of reads may be received by Fig.10 The sequence data interface 1008 is shown receiving.

[0162] In 1206, the plurality of reads are clustered into a plurality of clusters by sequence similarity. Similar sequences are likely to originate from sequencing of the same DNA strand (which is not identical due to errors introduced by the polynucleotide sequencer). Therefore, a cluster should represent all reads from the same DNA strand. Recall that a polynucleotide sequencer can sequence multiple different DNA strands simultaneously, and thus the raw output of sequence data from the polynucleotide sequencer may include reads corresponding to multiple different DNA strands. In implementation, clustering may be performed by Fig.10 The clustering module 1012 is shown executing.

[0163] At 1208, a cluster is selected from a plurality of clusters. A cluster contains a clustered set of reads. If the clustering is accurate, all reads in the clustered set of reads are from sequencing of the same DNA strand. At this point, before additional analysis is performed, the cluster is identified solely by the characteristics of the reads that are clustered together. Thus, in some implementations, each cluster is analyzed in turn, and the order in which the individual clusters are selected can be arbitrary. Multiple clusters in a cluster can also be analyzed in parallel. In an implementation, the selection of a cluster can be determined by Fig.10 The analysis may be performed by the clustering module 1012 as shown. The analysis may continue until trace reconstruction is performed on all clusters from the plurality of clusters.

[0164] In 1210, the clustered set of reads is aligned at a position of comparison across the clustered set of reads. In an implementation, the position of comparison may be a first position shared between the clustered set of reads. Thus, the original alignment may define a first position of comparison. The first position may be the leftmost position (corresponding to the 5' end), or alternatively may be the rightmost position (corresponding to the 3' end). In an implementation, the alignment may be by Fig.10 The read alignment module 1014 is shown to execute.

[0165] In 1212, the majority consensus base call at the first compared position is determined. The majority consensus base call is based at least in part on the most common base call across the clustered set of reads. The majority consensus base call can also be based in part on an error distribution associated with a polynucleotide sequencer. That is, the base calls can be weighted based on the associated error distribution (e.g., more certain base calls are more valuable, and fewer certain base calls have less impact on determining the majority consensus base call).

[0166] At 1214, variant reads from the clustered set of reads are identified. Variant reads have base calls that differ from the majority consensus base call at the compared position. In implementation, variant reads may be identified by Fig.10 The variant read identification module 1018 is shown to identify the variant reads.

[0167] Now move to Fig. 12B In 1216, the consensus string of base calls in the look-ahead window is identified. The consensus string is based on base calls of a subset of the clustered set of reads that have a majority of consensus base calls (i.e., non-variant reads) at the compared position. The look-ahead window is adjacent to the compared position. In some implementations, the look-ahead window can be two or three positions long.

[0168] At 1218, based at least in part on base calls in the look-ahead window of variant reads that match the consensus string of base calls, it is determined that the error type of the variant read at the compared position is a substitution. An example of this relationship is Figure 3 In implementation, the error type may be determined by the error classification module 1020 .

[0169] In 1220, the position of comparison of the variant chain is moved forward by one position.

[0170] At 1222, based at least in part on a series of base calls in the variant reads, including the base call at the compared position, and one or more base calls after the compared position that match the consensus string of base calls, it is determined that the error type of the variant read at the compared position is a deletion. An example of this relationship is Figure 4 In implementation, the error type may be determined by the error classification module 1020 .

[0171] In 1224, the position of comparison of the variant chain is shifted forward by zero positions.

[0172] At 1226, it is determined that the error type of the variant read at the compared position is an insertion. An insertion error is identified based at least in part on a base call in the variant read after the compared position that matches the majority consensus base call and a series of base calls in the variant read starting two positions after the compared position that matches the consensus string of base calls. An example of this relationship is described in Figure 5 In implementation, the error type may be determined by the error classification module 1020 .

[0173] In 1228, the position of comparison of the variant chain is moved forward two positions.

[0174] In 1230, the compared positions of reads in a subset of the clustered set of reads (ie, non-variant reads) are advanced forward by one position.

[0175] In 1232, a single consensus output sequence is determined by clustering the read segments. In an implementation, the single consensus output sequence is determined by Fig.10 The common output sequence generator 1024 is shown generating.

[0176] In 1234, the single consensus output sequence is converted to binary data. This may be the final manipulation of the information derived from the DNA strand before being used again as a digital computer file. In implementation, the change from sequence data to binary data may be accomplished by Fig.10 The conversion module 1028 is shown to perform.

[0177] Fig.13 A process 1300 is shown for identifying a portion of a read containing a burst error and generating a consensus output sequence using a portion of the read other than the portion containing the burst error. As an alternative to discarding the entire read, process 1300 can be performed with process 1100 or 1200. As described above, the read can be generated from a DNA strand storing digital data.

[0178] At 1302, the start of a portion of a read containing a burst error is identified. The start of the portion of a read containing a burst error can be identified by detecting the error itself. For example, in a read that does not match the common output sequence (variant read), a position that is not classified as an insertion, deletion, or substitution may be interpreted as the start of the portion of the read containing a burst error. An example of identifying an adventitious error is shown in FIG. Figure 6 is shown in .

[0179] In 1304, the end of the location containing the burst error is identified. Identifying the end of the location containing the burst error will find the position of the location beyond the burst error in the read segment, so that further analysis of the read segment will find a sequence that can be compared with other read segments in the cluster. The end of the burst error is identified by a candidate position in the read segment. The candidate position is connected to one or both sides of the backward matching region (mb) or the forward matching region (mf). As described above, the backward matching region matches the consensus sequence generated from other read segments in multiple read segments in the same cluster, and the forward matching sequence matches the sequence generated from the sequences of other read segments in the cluster by majority voting.

[0180] In 1306, the portion of the read containing the burst error is omitted from the generation of the common output sequence. Therefore, the corresponding positions of the read between the position identified as the start of the burst error and the position identified as the end of the burst error do not contribute to the determination of the common output sequence. For these positions in the common output sequence, base calls are determined based on some or all of the other reads in the cluster using techniques such as those illustrated in Figures 11 and 12.

[0181] In 1308, the common output sequence is generated using the portion of the read segment containing the burst error on either side. Therefore, even if a portion of the read segment is not used to generate the common output sequence, other portions of the read segment are also used. By comparing the values ​​of multiple read segments at the compared position, while aligning multiple read segments relative to each other based on insertions, deletions and substitutions, the portions on either side of the portion of the read segment with the burst error are used to generate the common output sequence. The common output sequence can be generated by the common output sequence generator 1024.

[0182] Fig.14 A process 1400 is shown for determining whether a read with an undetermined error can be "brought back" by finding a matching sequence below the read. As an alternative to discarding the entire read if the read includes an undetermined error, process 1400 can be performed with process 1100 or 1200. An undetermined error can be a burst error that is not classified as any of an insertion, deletion, or substitution.

[0183] In 1402, a wild error is identified in a read that is one of a plurality of reads of a polynucleotide sequencer. The plurality of reads may be a group of reads in the same cluster as formed, for example, by a clustering module 1012. The wild error is an error relative to a common output sequence of the plurality of reads. The common output sequence may be generated by a common output sequence generator 1024.

[0184] The uncertainty error can be Figure 6In addition, the error is identified as an undetermined error can be performed by the error classification module 1020. The undetermined error is located at the first position in the read segment. For example, the undetermined error may be located at position 100 in a read segment with a total length of 200 base pairs. Therefore, before determining that there is an undetermined error at position 100, the common output sequence can be identified for positions 1 to 99.

[0185] At 1404, a search window is defined that includes a second position in the read segment that is located at least a distance delayed beyond the first position. In an implementation, the second position is determined by k 0 In an implementation, the length of the delay can be between 4 and 10 positions. For example, if the delay distance is 8 positions, then the search window is a window that includes position 108 in the example read. The search window may include some positions where the base calls are error-free.

[0186] Limiting the search to the search window rather than the entire rest of the read segment focuses the search for matching sequences in the vicinity of where a matching sequence (if present) should be found. Shorter search windows make it less likely that a matching sequence will be identified, but longer search windows increase the probability of false positive matches. In implementations, the length of the search window can be between 9 and 20 positions. As described above, the search window can include two parts: a search window for k and a search window for k. 0 The backward search window (bsw) of the previous position and the 0 The forward search window (fsw) of the positions after bsw is defined. The lengths of bsw and fsw can be the same as each other. For example, both bsw and fsw may be about 5 to 8 positions long. Therefore, the search window includes bsw, k 0 and fsw.

[0187] In 1406, the edit distance is calculated for the subsequence in the search window determined from the comparison sequence. The comparison sequence may include a backward match sequence (mb), which is a common output sequence of corresponding positions in some or all of the other reads in the plurality of sequence reads. The length of the backward match sequence may be about 1 to 11 positions.

[0188] The consensus output sequence of multiple sequence reads can be determined by a majority vote identifying the considered positions, taking into account insertion errors, deletion errors, and substitution errors, and proceeding sequentially from the currently considered base call position to the adjacent base call position. Thus, all other reads from the same cluster (or all other reads that also do not have an adventitious error at that position) can be used to generate the consensus output sequence as described above.

[0189] Comparison sequence can also include forward matching sequence (mf).As mentioned above, forward matching sequence is adjacent to backward matching sequence, and is positioned to exceed the positioning of backward matching sequence (that is, away from the positioning of indefinite error).Forward matching sequence corresponds to a part for comparison reads, and its total output sequence has not yet been determined.Because total output sequence has not yet been determined for the position included in forward matching sequence, the majority of total votes of base calls from other reads in cluster are used to determine whether there is matching.The majority of total votes simply identify the most frequent base calls at each position, without attempting to consider errors such as insertion, deletion and substitution.If forward matching sequence is included, the length of forward matching sequence is 1 to 10 positions.

[0190] In implementation, the sliding window can be moved over the search window, and the sequence within the sliding window can be used as a subsequence to be compared with the common output sequence. The edit distance can be calculated for each of the positions of the sliding window. Depending on its position, the sliding window can include relatively more or relatively fewer backward matching sequences and forward matching sequences.

[0191] At 1408, it is determined whether a match is found between any of the subsequences at corresponding positions and the consensus output sequence. A match is identified by an edit distance less than a threshold. For example, if the threshold is 0, then only exact matches are considered matches. If the threshold for the edit distance is 1, then two sequences with a single base call may be classified as a match.

[0192] If a match is not found, then process 1400 proceeds along the "no" path to 1410 and the read is removed from consideration. Lack of a match indicates that there is no subsequence within the search window that matches the consensus output sequence / majority consensus vote. Therefore, the read may have a large amount of error, and the most accurate consensus output sequence can be obtained by ignoring the read.

[0193] However, if a match is found, then process 1400 proceeds along the "yes" path to 1412. In 1412, it is determined whether multiple matches are found. Within the search window, there may be more than one subsequence that matches the corresponding base call of the common output sequence. If this is the case, then process 1400 proceeds along the "yes" path to 1414.

[0194] At 1414, one of the matching subsequences is selected. Of the two or more matching subsequences, one subsequence can be selected by any of a variety of techniques, such as random selection, selecting one of the subsequences that is closest to the beginning of the read (e.g., the 5'-end), selecting one of the subsequences that is farthest from where the read was made, or by another criterion.

[0195] A single subsequence in the search window with an edit distance less than a threshold is identified in 1416. This is a matching subsequence, indicating that after localization of the adventitious error, there is a portion of the read that can be aligned again with other reads in the cluster.

[0196] At 1418, a candidate position within the matching subsequence is selected. The length of the matching subsequence may be approximately 5 to 15 positions, and one of these positions is selected to serve as a candidate position. The candidate position may be Figure 7 The candidate position (k) shown is the same. In implementation, the candidate position can be the position in the backward matching sequence that is farthest from the position where the indefinite error begins. Therefore, if there is no forward matching sequence, the candidate position can be a position in the matching subsequence at the end of the matching subsequence. If there is a forward matching sequence, the candidate position can be a position in the backward matching sequence that is directly adjacent to the forward matching sequence.

[0197] In 1420, the position of comparison is set to be one position beyond the candidate position. The position of comparison can be Figure 6 The position of comparison shown is the same as position 600. Therefore, the location of the matching subsequence and the candidate position are used to determine how to realign the read with other reads after the ad-hoc error. The position of comparison is at the third position in the read, which can be indicated as being at k+1 using the notation introduced above.

[0198] In 1422, a common output sequence from the plurality of reads is determined at a third position. The position of the comparison is used to compare the reads in the cluster with other reads to determine the common output sequence as described above. Analysis of the reads can continue until the end of the reads is reached, and any variations in the reads from the output sequence can be handled as described earlier in this disclosure.

[0199] Many parameters (such as delay, bsw, fsw, mb and mf) that are used to identify matching subsequences and finally determine the position of comparison have variable lengths. The length of a given set that is used to analyze sequence data can be experimentally determined by testing the different values ​​of each of delay, search window, mb and mf. Each can vary within the range of values, and different combinations can be tested.

[0200] For each combination of lengths, the number of clusters generated from the reads produced by the polynucleotide sequencer is determined. Generating a large number of clusters indicates that more reads can be used even with uncertain errors, and more data output by the polynucleotide sequencer helps to generate a common output sequence. Therefore, the combination of values ​​that results in the formation of the maximum number of read clusters from multiple sequence reads can be used as the length of the parameter.

[0201] Example

[0202] The technology described in the present disclosure is tested on three different files generated by nanopore sequencing, with sizes of 32kB, 115kB and 1500kB. Each of these files contains digital information encoded with the sequence of the DNA chain. This is not synthetic data, but an actual computer file encoded with DNA. Therefore, the qualitative test of the technology is partly a test of the ability to successfully recover files from the DNA chain. The quantitative measurement is the number of clusters correctly recovered from the output of the polynucleotide synthesizer and the coverage depth of sequencing.

[0203] This novel technique, which ignores portions of a read containing adventitious errors and uses portions that do not contain adventitious errors (the "new technique"), was compared to an earlier technique that attempted to classify errors as insertions, deletions, or substitutions and discarded the entire read from further consideration if the error could not be classified as one of these three types (the "old technique"). For each of the two smaller files, the new technique produced better results as measured by cluster recovery and coverage.

[0204] The number of correctly recovered clusters indicates the number of unique DNA chains for which the analysis can generate a consensus sequence. The total number of clusters formed includes clusters that do not produce a consensus sequence. Therefore, recovering a larger number of clusters indicates recovering sequences from a larger number of DNA chains. Comparison is presented in Table 1 below.

[0205]

[0206]

[0207] Table 1. Comparison of cluster recovery between two different techniques for trace reconstruction from multiple noisy reads.

[0208] Therefore, the new technology correctly recovers 2% to 4% more clusters. The largest file of 1500kB cannot be decoded using the old technology (therefore, there is no comparative cluster recovery value), but has been successfully decoded using the new technology. By trying various values ​​of delay, search window forward (swf), search window backward (swb), backward match (mb) and forward match (mf), about 5000 different parameter combinations are tested to see which combination will generate the maximum number of correct clusters. For all test results shown in Table 1, the distance threshold (dt) is set to zero, so an exact match is required. Swf and swb are set to the same value indicated in the sw(f / b) column. The parameters that produce the maximum number of correct clusters for 32kB files and 115kB files are shown in Table 2 below. Four different parameter combinations that successfully decoded 1500kB files are also shown in Table 2.

[0209] File size (kbytes) Delay Sw(f / b) mb mf 32 8 8 7 0 115 6 5 1 8 1500 5 5 5 5 1500 8 8 5 5 1500 8 8 8 0 1500 8 8 0 8

[0210] Table 2. Comparison of parameter values ​​for new techniques to successfully decode a file and maximize the number of correctly recovered clusters.

[0211] For the tested files, the delay length is 5 to 8, and the search window (forward and backward) is also 5 to 8. As indicated by the example of length 0, the forward match region or the backward match region can be omitted. However, for all the examples shown in Table 2, the combined length of the backward match region and the forward match region is 7 to 10.

[0212] In addition, the new technique was able to recover files with lower coverage levels. For a 32kB file, it was successfully decoded with 22× coverage using the new technique, but the old technique failed to decode it with 38× coverage. A 115kB file was successfully decoded with 27× coverage, but 32× coverage was insufficient using the old technique. For a large file of 1500kB, it was successfully decoded with 27.2× coverage. Testing with a random subsample of 1500kB files determined that it could be recovered with as low as 23.1× coverage. Using the new technique, all test samples could be recovered with 27× or less coverage, while the old technique required a coverage depth of well over 30× before the file could be successfully recovered. Reduced coverage allows the digital information stored in DNA to be decoded with less sequencing and at less expense than higher coverage levels.

[0213] Illustrative Embodiments

[0214] The following clauses describe a number of possible embodiments for implementing the features described in this disclosure. The various embodiments described herein are not restrictive, nor are every feature from any given embodiment required to be present in another embodiment. Unless the context clearly indicates otherwise, any two or more embodiments may be combined together. As used herein, in this document, "or" means and / or. For example, "A or B" means A without B, B without A, or A and B. As used herein, "comprising" means including all listed features and may include the addition of other features that are not listed. "Substantially consisting of..." means including the listed features and those additional features that do not substantially affect the basic characteristics and novel characteristics of the listed features. "Consisting of..." means only the listed features, and excludes any features that are not listed.

[0215] Item 1. A method for generating a consensus sequence from multiple reads of a deoxyribonucleic acid (DNA) chain storing digital data, the multiple reads being generated with a coverage of less than 30× by a sequencing technology that introduces burst errors into the reads among the multiple reads, the method comprising: omitting a portion of the read containing the burst error from generation of the consensus output sequence; and generating the consensus output sequence using portions of the read from either side of the portion of the read containing the burst error.

[0216] Clause 2. The method of clause 1, wherein the sequencing technology is nanopore sequencing.

[0217] Clause 3. The method of any one of clauses 1 to 2, wherein the plurality of reads are generated with a coverage of less than 25×.

[0218] Clause 4. The method of any one of clauses 1 to 3, further comprising: identifying the start of a portion of the read containing a burst error by identifying positions in the read that do not match the consensus sequence and are not classified as insertions, deletions, or substitutions.

[0219] Item 5. The method of any one of items 1 to 4 also includes: identifying the end of the portion of the read segment containing the burst error by identifying the position in the read segment that is flanked by at least one of the following: a backward matching region that matches a consensus sequence, or a forward matching region that matches a sequence generated by majority voting from sequences of at least two other read segments among multiple read segments.

[0220] Clause 6. The method of any one of clauses 1 to 5, wherein the consensus sequence is generated by comparing values ​​for multiple reads at compared positions while aligning the multiple reads to each other based on insertions, deletions, and substitutions.

[0221] Clause 7. A computer-readable medium encoding instructions that, when executed by a processing unit, cause a computing device to perform the method of any of clauses 1 to 6.

[0222] Clause 8. A system comprising a processing unit configured to implement the method of any one of clauses 1 to 6 and a memory.

[0223] Item 9. A method comprising: in a read among multiple sequence reads from a polynucleotide sequencer, identifying an uncertain error relative to a common output sequence of the multiple sequence reads, the uncertain error being located at a first position; defining a search window including a second position in the read, the second position being located at a distance at least delayed beyond the first position; calculating an edit distance for a subsequence in the search window and a comparison sequence, the comparison sequence including a backward matching sequence, the backward matching sequence being a common output sequence of corresponding positions in at least two of the multiple sequence reads; identifying a subsequence having an edit distance less than a threshold in the search window; selecting a candidate position within the subsequence based on the length of the backward matching sequence; setting a position for comparison with at least two other reads from the multiple sequence reads to a third position in the read, the third position being a position beyond the candidate position; and determining a common output sequence from the multiple sequence reads at the third position.

[0224] Clause 10. The method of any of Clause 9, wherein an adventitious error is an emergent error that is not classified as any of an insertion, a deletion, or a substitution.

[0225] Clause 11. The method of any one of clauses 9 to 10, wherein a common output sequence of multiple sequence reads is determined by identifying a majority vote for a considered position while taking into account insertion errors, deletion errors, and substitution errors, and the common output sequence of multiple sequence reads proceeds sequentially from a currently considered base call position to an adjacent base call position.

[0226] Clause 12. The method of any one of clauses 9 to 11, further comprising: identifying a second subsequence in the search window having an edit distance less than a threshold, and selecting the subsequence or one of the second subsequences that is closest to the candidate position.

[0227] Clause 13. The method of any one of clauses 9 to 12, wherein the threshold value of the edit distance is 0.

[0228] Clause 14. The method of any one of clauses 9 to 13, wherein comparing the sequences further comprises a forward match sequence that is a majority consensus vote for base calls from at least two other reads from the plurality of sequence reads.

[0229] Clause 15. The method of clause 14, wherein the length of the delay is between 4 and 10 positions, the length of the search window is between 9 and 20 positions, the length of the backward matching sequence is between 1 and 11 positions, and the length of the forward matching sequence is between 1 and 10 positions.

[0230] Clause 16. The method of Clause 15, wherein the sum of the length of the backward matching sequence and the length of the forward matching sequence is between 5 and 15 positions.

[0231] Clause 17. The method of clause 15 or 16, wherein the length of the delay, the length of the search window, the length of the backward matching sequence, and the length of the forward matching sequence are each experimentally determined by testing different values ​​for each length and selecting a combination of values ​​that results in the maximum number of read clusters recovered from multiple sequence reads.

[0232] Clause 18. A computer-readable medium encoding instructions that, when executed by a processing unit, cause a computing device to perform the method of any of clauses 9 to 17.

[0233] Clause 19. A system comprising a processing unit configured to implement the method of any of clauses 9 to 17 and a memory.

[0234] Clause 20. A system for alignment of sequence reads, comprising: one or more processing units; a memory coupled to the one or more processing units; an alignment anchor module, stored in the memory and implemented by the one or more processing units to: identify an alignment anchor sequence having a predetermined sequence, the predetermined sequence being present in a first sequence read and a second sequence read, divide the first sequence read into a first sub-read and a second sub-read, the first sub-read extending from the beginning of the first sequence read to the alignment anchor sequence, the second sub-read extending from the alignment anchor sequence to the end of the first sequence read, and divide the second sub-read into a first sub-read and a second sub-read, the first sub-read extending from the beginning of the first sequence read to the alignment anchor sequence, The reads are divided into a third sub-read and a fourth sub-read, the third sub-read extending from the beginning of the second sequence read to the alignment anchor sequence, and the fourth sub-read extending from the alignment anchor sequence to the end of the second sequence read; and a read alignment module, which is stored in a memory and implemented on one or more processing units to: align the first sub-read with the third sub-read based on the beginning of the first sequence read, the beginning of the second sequence read, and the alignment anchor sequence, and align the second sub-read with the fourth sub-read based on the alignment anchor sequence, the end of the first sequence read, and the end of the second sequence read.

[0235] Clause 21. The system of clause 20, wherein the alignment anchor sequence has a sequence that is not aligned to itself.

[0236] Clause 22. The system of any of clauses 20 to 21, wherein the alignment anchor sequence is located between 45% and 55% of the distance from the beginning of the first sequence read to the end of the first sequence read, and is located between 45% and 55% of the distance from the beginning of the second sequence read to the end of the second sequence read.

[0237] Clause 23. The system of any one of clauses 20 to 22, further comprising an oligonucleotide synthesizer configured to read a first polynucleotide molecule with a first sequence comprising only one instance of an alignment anchor sequence, and to synthesize a second polynucleotide sequence, wherein the second polynucleotide sequence comprises only one instance of the alignment anchor sequence.

[0238] Clause 24. The system of any one of clauses 20 to 23 also includes a consensus output sequence generator module, which is stored in a memory and implemented on one or more processing units to: generate a first forward consensus sequence from the first sub-segment and the third sub-segment; generate a first reverse consensus sequence from the first sub-segment and the third sub-segment; combine the first forward consensus sequence and the first reverse consensus sequence into a first combined consensus sequence; generate a second forward consensus sequence from the second sub-segment and the fourth sub-segment; generate a second reverse consensus sequence from the second sub-segment and the fourth sub-read; combine the second forward consensus sequence and the second reverse consensus sequence into a second combined consensus sequence; append the second combined consensus sequence to the end of the first combined consensus sequence by alignment at the alignment anchor sequence; and delete the alignment anchor sequence.

[0239] Item 25. A system for aligning sequence reads, comprising: one or more components for processing; a memory coupled to the one or more components for processing; a component for identifying an alignment anchor sequence having a predetermined sequence, the predetermined sequence being present in a first sequence read and a second sequence read, a component for dividing the first sequence read into a first sub-segment extending from the beginning of the first sequence read to the alignment anchor sequence, and a second sub-segment extending from the alignment anchor sequence to the end of a first polynucleotide sequence, a component for dividing the second sequence read into a third sub-segment extending from the beginning of the second sequence read to the alignment anchor sequence, and a fourth sub-segment extending from the alignment anchor sequence to the end of the second polynucleotide sequence; a component for aligning the first sub-segment with the third sub-segment based on the beginning of the first sequence read, the beginning of the second sequence read, and the alignment anchor sequence, and a component for aligning the second sub-segment with the fourth sub-read based on the alignment anchor sequence, the end of the first sequence read, and the end of the second sequence read.

[0240] in conclusion

[0241] Although the subject matter has been described in language specific to structural features and / or methodological acts, it is to be understood that the subject matter defined in the appended claims is not necessarily limited to the specific features or acts described above. Rather, the specific features and acts are disclosed as example forms of implementing the claims.

Claims

1. A method for improving the sequencing accuracy of a polynucleotide sequencer, include: Sequencing the polynucleotide molecule using the polynucleotide sequencer to generate a plurality of sequence reads; generating a first consensus output sequence for a first portion of the plurality of sequence reads from the polynucleotide sequencer, the first portion being before an adventitious error; identifying, in a read in the plurality of sequence reads, the variable error at a first position relative to a majority consensus base call of the plurality of sequence reads, wherein the variable error is a base call that does not match the majority consensus base call at the first position and the variable error cannot be identified as an insertion, deletion, or substitution; defining a search window including a second position in the read segment, the second position being located at at least a delay beyond the first position, wherein the delay has a predetermined length; Calculating an edit distance between a subsequence in the search window and a comparison sequence, the comparison sequence comprising a backward matching sequence, the backward matching sequence being a common output sequence for a predetermined number of positions before the second position; Identifying subsequences in the search window having an edit distance equal to or less than a threshold; selecting a candidate position outside the backward matching sequence and within the subsequence; Determining a second common output sequence for a second portion of the plurality of sequence reads after the candidate position; as well as A single consensus output sequence is generated from the first consensus output sequence and the second consensus output sequence, the single consensus output sequence representing a nucleotide sequence of the polynucleotide molecule. 2 . The method according to claim 1 , wherein the uncertain error is a burst error including a plurality of adjacent errors.

3. A method according to claim 1, wherein the first common output sequence of the multiple sequence reads is determined by: identifying a majority vote for a considered position, taking into account insertion errors, deletion errors, and substitution errors, and proceeding sequentially from a currently considered base call position to an adjacent base call position.

4. The method according to claim 1, further comprising: include: A second subsequence in the search window having an edit distance less than the threshold is identified, and the subsequence or the second subsequence closest to the first position is selected. The method of claim 1 , wherein the threshold for the edit distance is 0. 6 .

6. The method according to claim 1, further comprising: include: A forward match sequence is generated from a majority consensus vote for base calls from reads of the plurality of sequence reads other than the read with the adventitious error, the forward match sequence being located after the second position, and wherein the comparison sequence also includes at least a portion of the forward match sequence.

7. The method of claim 6, wherein the length of the delay is between 4 and 10 positions, the length of the search window is between 9 and 20 positions, the length of the backward matching sequence is between 1 and 11 positions, and the length of the forward matching sequence is between 1 and 10 positions.

8. The method of claim 7, wherein the sum of the length of the backward matching sequence and the length of the forward matching sequence is between 5 and 15 positions.

9. The method of claim 6, wherein the plurality of sequence read segments encode binary data, wherein the length of the delay, the length of the search window, the length of the backward matching sequence, and the length of the forward matching sequence are respectively determined experimentally by: testing different values ​​for each length, and selecting a combination of values ​​that results in recovery of the binary data with a minimum coverage depth.

10. A system for improving the sequencing accuracy of a polynucleotide sequencer, the system include: a polynucleotide sequencer configured to generate a plurality of sequence reads by sequencing the polynucleotides; one or more processing units; a memory coupled to the one or more processing units; a variant read identification module stored in the memory and implemented by the one or more processing units to identify an error in a read in the plurality of sequence reads, wherein the error is a base call that does not match a majority consensus base call at a first position in an alignment with the sequence reads; an error classification module, stored in the memory and implemented by the one or more processing units, to determine that the error is an indeterminate error that cannot be identified as an insertion, deletion, or substitution; a common output sequence generator, stored in the memory and implemented by the one or more processing units to: generating a first common output sequence for a first portion of the plurality of sequence reads preceding the first position including the uncertain error; defining a search window including a second position in the read segment, the second position being located at at least a delay beyond the first position, wherein the delay has a predetermined length; Calculating an edit distance between a subsequence in the search window and a comparison sequence, the comparison sequence comprising a backward matching sequence, the backward matching sequence being a common output sequence for a predetermined number of positions before the second position; Identifying subsequences in the search window having an edit distance equal to or less than a threshold; selecting a candidate position outside the backward matching sequence and within the subsequence; Determining a second common output sequence for a second portion of the plurality of sequence reads after the candidate position; as well as A single consensus output sequence is generated from the first consensus output sequence and the second consensus output sequence, the single consensus output sequence representing a nucleotide sequence of the polynucleotide molecule. The system according to claim 10 , wherein the uncertain error is a burst error including a plurality of adjacent errors.

12. A system according to claim 10, wherein the consensus output sequence generator determines the first consensus output sequence of the multiple sequence reads by identifying a majority vote for a considered position, taking into account insertion errors, deletion errors, and substitution errors, and proceeding sequentially from a currently considered base call position to an adjacent base call position.

13. The system of claim 10, wherein the consensus output sequence generator identifies a second subsequence in the search window having an edit distance less than the threshold, and selects the subsequence or the second subsequence closest to the first position. The system of claim 10 , wherein the threshold for the edit distance is 0.

15. A system according to claim 10, wherein the consensus output sequence generator generates a forward match sequence from a majority consensus vote for base calls from reads of the multiple sequence reads other than the read with the adventitious error, the forward match sequence is located after the second position, and wherein the comparison sequence also includes at least a portion of the forward match sequence.

16. The system of claim 15, wherein the length of the delay is between 4 and 10 positions, the length of the search window is between 9 and 20 positions, the length of the backward matching sequence is between 1 and 11 positions, and the length of the forward matching sequence is between 1 and 10 positions.

17. The system according to claim 10, further comprising a conversion module stored in the memory and implemented by the one or more processing units to: convert the single common output sequence into binary data according to an encoding.

18. The system of claim 17, wherein the polynucleotide sequencer is configured to generate the plurality of reads having a coverage of less than 25× and the conversion module is configured to convert the single consensus output sequence into binary data with an accuracy sufficient to restore the binary data to a digital file.

19. The system of claim 10, wherein the polynucleotide sequencer is a nanopore sequencer.

20. A system for improving the sequencing accuracy of a polynucleotide sequencer, the system include: means for sequencing a polynucleotide molecule to generate a plurality of sequence reads; means for generating a first common output sequence for a first portion of the plurality of sequence reads, the first portion being prior to an adventitious error; means for identifying, in a read of the plurality of sequence reads, the variable error at a first position relative to a majority consensus base call of the plurality of sequence reads, wherein the variable error is a base call that does not match the majority consensus base call at the first position and the variable error cannot be identified as an insertion, deletion, or substitution; means for defining a search window including a second position in the read segment, the second position being located at at least a delay beyond the first position, wherein the delay has a predetermined length; means for calculating an edit distance for a subsequence in the search window and a comparison sequence, the comparison sequence comprising a backward match sequence, the backward match sequence being a common output sequence for a predetermined number of positions before the second position; means for identifying subsequences in the search window having an edit distance equal to or less than a threshold; means for selecting a candidate position outside the backward matching sequence and within the subsequence; means for determining a second common output sequence for a second portion of the plurality of sequence reads after the candidate position; as well as Means for generating a single consensus output sequence from the first consensus output sequence and the second consensus output sequence, the single consensus output sequence representing a nucleotide sequence of the polynucleotide molecule.

Images