RU2799778C9 - Method for determining indicator correlated with probability that two mutated sequence readings are from the same sequence containing mutation - Google Patents

Abstract

FIELD: bioinformatics; biotechnology.

SUBSTANCE: computer-implemented method for determining the sequence of at least a portion of at least one template target nucleic acid by determining whether two mutated sequence reads originate from the same mutated sequence. Multiple mutated sequence reads are performed. Moreover, each reading of a mutated sequence corresponds to a subsequence of the sequence containing the mutations. A general minimizing function is applied for each mutated sequence read, thus determining one or more appropriate minimizers for each mutated sequence read. The positions of one or more corresponding minimizers in each mutated sequence read are determined. The positions of one or more mutations in each mutated sequence read are determined. At the same time, for at least two mutated reads of the sequence with a common minimizer, the number of mutations with the same position and/or with a mismatched position is counted, when the corresponding minimizers are aligned in order to determine the indicator correlated with the probability that the indicated at least two mutated sequence reads come from the same sequence containing the mutation. At least two mutated sequence reads are assembled based on the specified index. Sequences of at least a portion of at least one template target nucleic acid are determined based on said assembly. A method for determining at least a portion of the sequence of at least one template target nucleic acid molecule is also described, including using the above method after the step of sequencing regions of at least one mutated template target nucleic acid molecule to obtain a plurality of mutated sequence reads.

EFFECT: accurate sequencing of nucleic acids, as well as fast and accurate sequence assembly from short sequence reads.

29 cl, 5 dwg

Description

Scope of the invention

The invention relates to a computer-implemented method for determining an indicator correlated with the probability that two mutated sequence reads originate from the same sequence containing mutations, a method for obtaining at least a portion of the sequence, and a method for determining the sequence of at least one template nucleic acid molecule - targets.

Prerequisites for the creation of the invention

The ability to sequence nucleic acid molecules is a tool that is very useful in many different applications. However, it can be difficult to determine the exact sequences of nucleic acid molecules that contain problem structures, such as nucleic acid molecules that contain repetitive regions. It can also be difficult to resolve structural features such as the haplotype structure of diploid and polyploid organisms and structural variants in the genomes of these organisms.

Many of the more recent techniques (so-called next-generation sequencing techniques) are capable of accurately sequencing only short nucleic acid molecules. Next generation sequencing techniques can be used to sequence longer nucleic acid sequences, but this is often complex and expensive. Next generation sequencing techniques can be used to generate short sequence reads corresponding to the sequences of regions of the nucleic acid molecule, and the complete sequence can be assembled from these short sequence reads. If the nucleic acid molecule contains repetitive regions, it may not be clear to the user whether two sequence reads having similar sequences correspond to sequences of two repeats in a longer sequence, or two replications of the same sequence. Similarly, a user may wish to sequence two similar nucleic acid molecules simultaneously, and it may be difficult to determine whether two sequence reads having similar sequences correspond to sequences from the same parent nucleic acid molecule or two different parent nucleic acid molecules.

Sequencing-assisted mutagenesis (SAM) techniques can assist in assembling sequences from short sequence reads. In general, SAM involves introducing mutations into the template sequences of target nucleic acids. Insertion mutation schemes can assist the user of the method in assembling nucleic acid molecule sequences from short sequence reads.

For example, in cases where template nucleic acid molecules contain repeat regions, the repeats can be distinguished from each other by different mutation patterns, allowing the repeat regions to be resolved and assembled correctly.

In general, SAM techniques include introducing mutations into copies of a template target nucleic acid molecule to obtain a mutated template target nucleic acid molecule and/or one or more mutated sequences, sequencing one or more mutated sequences to obtain SAM data , including mutated sequence reads, and then assembling sequences from mutated sequence reads based on appropriate mutation patterns. Since different mutated copies will contain mutations at different positions, the assembled sequence can characterize the original template nucleic acid molecule.

However, there remains a need for more reliable and/or more computationally efficient methods for processing SAM data.

Statement of the Invention

The authors of the present invention have developed new improved methods for processing SAM data, including mutated sequence reads. Thus, in one aspect of the invention, a computer-implemented method is provided for determining a score correlated with the probability that two mutated sequence reads are from the same mutated sequence. The method includes receiving multiple mutated sequence reads. Each mutated sequence read corresponds to a subsequence of the mutated sequence. A sequence containing mutations contains mutations compared to a sequence that does not contain mutations. The method further includes applying a common minimizer function for each mutated sequence read, thereby determining one or more appropriate minimizers for each mutated sequence read. The method further includes determining the positions of one or more corresponding minimizers in each mutated sequence read. The method further includes determining the positions of one or more mutations in each mutated sequence read. For at least two mutated sequence reads with a common minimizer, the method further comprises counting the number of matched and/or mismatched mutations when the respective minimizers are aligned.

In another aspect of the present invention, a method is provided for obtaining at least a portion of the sequence of a template target nucleic acid molecule.

In another aspect of the present invention, a method is provided for determining at least a portion of the sequence of at least one template target nucleic acid molecule.

Additional aspects of the present invention are presented in the dependent claims and in the detailed description.

Brief description of graphic materials

On FIG. 1 shows an embodiment of a method for detecting at least a portion of at least one template target nucleic acid molecule in accordance with the present invention.

On FIG. 2 shows an embodiment of a method for determining a score correlated with the probability that two mutated sequence reads are from the same mutated sequence, in accordance with the present invention.

On FIG. 3 is an example of the step of determining the positions of one or more mutations in a mutated sequence read.

On FIG. 4A shows a comparative example of an assembly from short reads of the 2.3 Mb Arcobacter butzlerii genome. without using the method of the present invention.

On FIG. 4B shows an example assembly of the 2.3 Mbp genome of Arcobacter butzlerii. using the method of the present invention.

On FIG. 5 presents experimental data on the effect of coverage depth by short readings of a long template on the results of the method of the present invention.

Detailed description of the invention

General definitions

Unless otherwise indicated, all technical and scientific terms used in this document have a generally accepted meaning, understandable to a person skilled in the field to which the present invention relates.

In general, the term "comprising" means "including, without limitation". For example, the phrase "method comprising [certain steps]" should be interpreted to mean that the method includes the steps listed, but additional steps may be performed.

In some embodiments of the invention, the word "comprising" is replaced by the phrase "consisting of". The term "consisting of" implies a limitation. For example, the phrase "method consisting of [certain steps]" should be understood to mean that the method includes these steps, and additional steps are not performed.

In some aspects, the invention provides a method for determining or obtaining at least a portion of the sequence of at least one template target nucleic acid molecule. The method can be used to determine or obtain the complete sequence of at least one template target nucleic acid molecule. In an alternative embodiment, the method can be used to determine or obtain a partial sequence, ie. sequences of a portion of at least one template target nucleic acid molecule. For example, if it is not possible or difficult to determine the complete sequence, the user may decide that the sequence of a portion of at least one template target nucleic acid molecule is useful or even sufficient for his purpose.

For the purposes of the present invention, the term "nucleic acid molecule" (or "unmutated nucleic acid molecule") refers to the polymeric form of nucleotides having any length. Nucleotides can be deoxyribonucleotides, ribonucleotides, or analogs thereof. Preferably at least one nucleic acid molecule consists of deoxyribonucleotides or ribonucleotides. Even more preferably, at least one nucleic acid molecule consists of deoxyribonucleotides, ie. at least one nucleic acid molecule is a DNA molecule.

A "Template Target Nucleic Acid Molecule" can be any nucleic acid molecule that the user wishes to sequence.

At least one "template target nucleic acid molecule" may be single stranded or may be part of a double stranded complex. If at least one template target nucleic acid molecule consists of deoxyribonucleotides, it may form part of a double-stranded DNA complex. In this case, one strand (eg, the coding strand) will be considered to be at least one template target nucleic acid molecule, and the other strand is a nucleic acid molecule that is complementary to at least one template target nucleic acid molecule. The at least one template target nucleic acid molecule may be a DNA molecule corresponding to a gene, may contain introns, may be an intergenic region, may be an intragenic region, may be a genomic region spanning multiple genes, or may actually represent the entire genome of an organism.

For purposes of the present invention, the term "mutated nucleic acid molecule" or "mutated template target nucleic acid molecule" refers to a "nucleic acid molecule" or "template target nucleic acid molecule" that has been mutated. The mutations may be substitution mutations, optionally transitional mutations. For the purposes of the present invention, the term "mutation by type of substitution" should be interpreted as that the nucleotide is replaced by another nucleotide. For example, converting an ATCC sequence to an AGCC sequence introduces one substitution type mutation. For the purposes of the present invention, the term “transition-type mutation” should be interpreted as that nucleotide A is replaced by nucleotide G and vice versa (i.e. A⇔G mutations) or that nucleotide C is replaced by nucleotide T and vice versa (i.e. mutations C ⇔T).

The phrase "introducing mutations into at least one template target nucleic acid molecule" means subjecting at least one template target nucleic acid molecule to the second of a pair of samples of conditions in which at least one template target nucleic acid molecule mutates. This can be achieved using any suitable method. For example, mutations can be introduced by chemical mutagenesis and/or enzymatic mutagenesis.

For the purposes of the present invention, the "mutated sequence" corresponds to at least a portion of the nucleotide sequence in the "mutated nucleic acid molecule" or "mutated template target nucleic acid molecule". A "mutated sequence" may also be referred to as a "mutated sequence". A “mutated sequence” is referred to herein as μ ⁱ , and a set (i.e., several) of “mutated sequences” is referred to as M, where μ ¹ ...μ ⁿ ∈ M. A “mutated sequence” corresponds to at least parts of a nucleotide sequence in a "nucleic acid molecule" or "template target nucleic acid molecule". A "non-mutated sequence" may also be referred to as a "non-mutated sequence". A “mutation-free sequence” is referred to herein as S ⁱ , and a set (i.e., several) of “mutation-free sequences” is referred to as S, where S ¹ ...S ⁿ ∈ S. A “mutated sequence” and “not mutated sequence" thus may correspond to at least a portion of the sequence of the nucleic acid molecule of the nucleotides (nt) of adenine (A), thymine (T), guanine (G), and cytosine (C). Such a chromosomal sequence may be in the range of 10 ³ to 10 ⁹ nucleotides (nt) or more in length.

For the purposes of the present invention, a “mutated sequence read” corresponds to a subsequence of a “mutated sequence”, i.e., a “mutated sequence read” may be substantially identical to at least a subsequence of the “mutated sequence”, but it contains mutations compared to containing mutations. sequence and may contain additional small differences due to read errors. A "mutated sequence read" is denoted as ρ ⁱ , and a set (i.e. several) of "mutated sequence reads" are denoted as P, where ρ ¹ ...p ⁿ ∈ P. An "unmutated sequence read" corresponds to a subsequence of "a sequence that does not contain mutations" , i.e. A "non-mutated sequence read" may be substantially identical to a subsequence of the "non-mutated sequence" except for read errors during sequencing. An "unmutated sequence read" is denoted as ^ri and a set (i.e., several) of "unmutated sequence reads" is denoted as R, where r ¹ ...r ⁿ ∈ R. A "mutated sequence read" can be obtained by sequencing the "mutated sequence read" region. template target nucleic acid molecule" and a "non-mutated sequence read" can be obtained by sequencing the "template target nucleic acid molecule" region. The sequence read may have a length that is less than the sequence, such as about 150 nt in length.

Sequence Analysis Method 10

On FIG. 1 shows a method 10 for determining at least a portion of at least one template target nucleic acid molecule in accordance with the invention.

Method 10 for determining at least a portion of at least one template target nucleic acid molecule may include sample preparation step S110. Sample preparation step S110 may include providing a pair of template target nucleic acid molecules and introducing mutations into one of the pair of template target nucleic acid molecules to obtain a mutated template target nucleic acid molecule. Sample preparation step S110 may include any known techniques for providing a template target nucleic acid molecule and a mutated template target nucleic acid molecule.

The method 10 for determining at least a portion of at least one template target nucleic acid molecule may further include a sequencing step S120. Sequencing step S120 includes sequencing regions of at least one mutated template target nucleic acid molecule, thereby obtaining a plurality of mutated P sequence reads. In addition, the sequencing step S120 may include sequencing regions of at least one (unmutated) template target nucleic acid molecule (a template target nucleic acid molecule that corresponds to a mutated template target nucleic acid molecule) resulting in multiple unmutated reads of the R sequence . Step S120 may include any known techniques for obtaining multiple mutated P sequence reads.

Method 10 for determining at least a portion of at least one target nucleic acid template molecule includes step 200 or method 200 for determining whether two mutated sequence reads ρ ⁱ , ρ ^j are (or are) derived from the same mutated sequence μ ⁱ . Determining whether two mutated reads ρ ⁱ , ρ ^j of a sequence are (or are) derived from the same mutated sequence μ ⁱ includes determining whether two mutated reads ρ i , ρ ^j are obtained (or are occurring ⁾ from the same or similar or overlapping part containing the mutation sequence μ ⁱ , i. e. whether both mutated sequence reads ρ ⁱ , ρ ^j contain a subsequence that corresponds to the same part of the mutated sequence μ ⁱ . Method 200 is a computer-implemented method and may be implemented by a computer processor. The method 200 generates a score correlated with the probability that two mutated sequence reads ρ ⁱ , ρ ^j come from the same mutated sequence μ ⁱ .

The method 10 for determining at least a portion of at least one template target nucleic acid molecule may further comprise the step of assembling a sequence S300. Sequence assembly step S300 includes assembly or reconstruction of at least a portion of the sequence μ ⁱ , S ⁱ . A mutated sequence μ ⁱ can be obtained by assembling a set of mutated reads P of a sequence based on a score correlated with the probability that the respective two mutated sequence reads ρ ⁱ , ρ ^j come from the same mutated sequence μ ⁱ . This can be achieved, for example, by grouping a plurality of mutated P sequence reads into groups corresponding to sequences containing μ ⁱ mutations, and then assembling each group separately to reconstruct some or all of the individual μ ⁱ mutated sequences. The mutant-free sequence S ⁱ can be obtained by correcting the errors of the mutated sequence μ ⁱ , for example by deriving the most likely mutation-free sequence S ⁱ from the mutated sequence μ ⁱ using a plurality of unmutated reads R of the sequence. The sequence assembly step S300 may include any known methods of assembling a mutated sequence μ ⁱ from a plurality of mutated sequence reads P based on a score correlated with the probability that the respective two mutated sequence reads ρ ⁱ , ρ ^j are from the same sequence containing mutations of the sequence μ ⁱ .

On FIG. 2 shows a method 200 for determining whether two mutated sequence reads ρ ⁱ , ρ ^j originate from the same mutated sequence μ ⁱ in accordance with the present invention.

The method 200 includes the step S210 of receiving a set of mutated reads ρ ¹ ...p ⁿ ∈ P of the sequence. Each mutated reading of the ρ ⁱ sequence corresponds to a subsequence of the mutated sequence μ ⁱ . The mutated sequence μ ⁱ contains mutations, eg substitution type mutations, optionally transition type mutations, compared to a non-mutated sequence S ⁱ . The mutated sequence μ ⁱ may be at least a portion of the sequence of the mutated template target nucleic acid, and the non-mutated sequence may be at least a portion of a (unmutated) template target nucleic acid molecule, the mutated template target nucleic acid molecule being obtained by introducing mutations, eg substitution type mutations, optionally transition type mutations, into the template target nucleic acid molecule. Each subsequence of the mutated sequence μ ⁱ may be at least part of the sequence of a fragment of a mutated template target nucleic acid molecule. Each subsequence of the mutant-free sequence S ⁱ may be at least part of the sequence of a fragment of the template target nucleic acid molecule. Step S210 of receiving a plurality of mutated P sequence reads may include receiving a plurality of mutated P sequence reads directly from a sequencer used for sequencing a mutated target nucleic acid template molecule, or receiving a plurality of mutated P sequence reads from a data store that stores a plurality of mutated P sequence reads.

The method 200 further includes the step S220 of applying a common minimizing function for each mutated read ρ ⁱ of the sequence. The application of a common minimizing function determines one or more appropriate minimizers for each mutated read ρ ⁱ of the sequence. The method 200 further includes the step S222 of determining the positions of one or more appropriate minimizers in each mutated read ρ ⁱ of the sequence.

In a preferred embodiment, method 200 includes the step of S224 allocating mutated reads of the P sequence to respective minimizer groups. A mutated reading ρ ⁱ of a sequence for which more than one minimizer is defined can be placed in multiple corresponding minimizer groups.

The method 200 further includes the step S230 of determining the positions of one or more mutations in each mutated reading ρ ⁱ of the sequence. Step S230 of determining the positions of one or more mutations in each mutated reading of the ρ ⁱ sequence may be performed before, after or simultaneously with steps S220, S222 and S224 related to the overall minimizing function.

For at least two mutated sequence reads ρ ⁱ , ρ ^j with a common minimizer, method 200 further includes counting the number of matched and/or mismatched mutations when the respective minimizers are aligned, i.e. when the positions of the nucleotides of one mutated reading of the ρ ⁱ sequence are shifted relative to the positions of the nucleotides of another mutated reading of the ρ ^j sequence in such a way that the position of the minimizer of one mutated reading of the ρ ⁱ sequence is identical to the position of the minimizer of the other mutated reading of the ρ ^j sequence. The number of matched and/or mismatched mutations can be an indicator correlated with the probability that two mutated sequence reads are from the same mutation containing sequence. In an alternative embodiment, method 200 may include the additional step of S242 determining a score correlated with the probability that two mutated sequence reads are from the same mutated sequence based on the number of matched and/or mismatched mutations.

Step S210 of receiving a plurality of mutated sequence reads

Step S210 includes receiving a set of mutated reads ρ ¹ ...ρ ⁿ ∈ P of the sequence. Step S210 may further include receiving a plurality of unmutated reads r ¹ ... r ^rn ∈ R of the sequence. Each mutated reading of the ρ ⁱ sequence may correspond to a subsequence of the mutated sequence μ ⁱ . Each unmutated read of the r ⁱ sequence may correspond to a subsequence of the unmutated sequence S ⁱ .

A mutated sequence μ ⁱ can be obtained by introducing mutations into a non-mutated sequence S ⁱ . Each mutated reading of the ρ ⁱ sequence can thus contain mutations, i.e. correspond to the region of the mutated target nucleic acid template molecule that includes mutations, i.e. match a subsequence containing the mutation sequence. In one embodiment, each mutated read of the ρ ⁱ sequence contains substitution type mutations, i. e. corresponds to the region of the mutated template target nucleic acid molecule, which includes substitution type mutations. In a preferred embodiment, the substitution mutations are transitional mutations such that each mutated reading of the ρ ⁱ sequence contains a transitional mutation, i. e. corresponds to the region of the mutated template target nucleic acid molecule, which includes transitional mutations.

Each nucleotide of each reading of the ρ ⁱ , r ⁱ sequence is preferably encoded in binary using two bits. This is advantageous in particular when multiple mutated P sequence reads contain transition type mutations (A⇔G and C⇔T) such that one of the two bits (e.g. the first bit) determines whether the nucleotide is purine (A or G) or pyrimidine (T or C). For example, nucleotides can be encoded in binary using the following format: A: 00, G: 01, C: 10, and T: 11. This encoding will be used throughout this specification. However, it will be apparent that the present invention is not limited to this coding and that the present invention can be easily implemented using any other nucleotide encoding.

Each read of the ρ ⁱ , r ⁱ sequence can be encoded to account for homopolymer errors in the read of the ρ ⁱ , r ⁱ sequence. Homopolymer errors occur when a pass misreads the length of the same nucleotide, for example, the sequence TAAAAGC can be misread as TAAGC because it is difficult for the sequencer to determine the number of A's if multiple A's are encountered in a pass. To account for such homopolymer errors, runs with several identical nucleotides can be encoded as one instance of a nucleotide. In an alternative embodiment, homopolymer errors may be accounted for during post-processing (i.e., not during initial encoding) of the ρ ⁱ , r ⁱ reads of the sequence, e.g., by encoding any k-mers used in method 200 and/or any seeding patterns used in step S230, so that passages with multiple identical nucleotides are encoded as a single nucleotide instance.

Steps S220 and S222: common minimizing function

A minimizer is a k-mer from a set of k-mers that satisfies the general minimizing function min(⋅) on the set of k-mers.

For the purposes of this application, k-mer is a nucleotide subsequence of length k. k-measures starting from position i in the sequence S=[S ₁ , S ₂ , …, S _n-1 , S _n ] of length n is denoted as k(S _i ), and k(S _i )=[S _i , S _i+1 , …, S _i+k-1 ]. The set of k-measures in the sequence S with initial positions between i and j is denoted as k(S _i …S _j ). The minimizer of all k-mers with initial positions in the range i to j of the sequence S will be denoted as min(k(S _i …S _j )).

The general minimizing function min(⋅) is used to determine one or more minimizers (i.e., one or more representative k-mers) from a set of k-mers, preferably from all or substantially all of the k-mers generated by reading ρ ⁱ , r ⁱ sequences, i.e. k-mers, preferably from all or essentially all k-mers that exist in reading the ρ ⁱ , r ⁱ sequence. For the purposes of the present invention, the set of k-mers that exist in the ρ ⁱ , r ⁱ read of the sequence may contain the reverse complement k-mers of the ρ ⁱ , r ⁱ read of the sequence. Preferably each minimizer is a k-mer of length equal to or greater than 5 (i.e. 5-mer or more), preferably equal to or greater than 10 (i.e. 10-mer or more), even more preferably equal to or greater than 15 (i.e. 15 measures or more). Each minimizer may be a k-mer of length less than 50, optionally less than 30, optionally less than 25. If the general minimization function min(⋅) is used to define longer minimizers, then a particular minimizer is more likely to be representative of a particular part of the sequence, those. less likely that the minimizer will appear in several separate and unrelated sections of the sequence. Setting an upper limit on the size of minimizers reduces the risk that minimizers will contain sequencing errors.

Step S220 of applying the general minimizing function min (⋅) may include identifying in the corresponding mutated reading ρ ⁱ of the sequence of one or more k-mers that is(s) listed first(s) in the ordered list of possible k-mers. One or more minimizers defined for the corresponding mutated reading of the ρ ⁱ sequence may represent a failure of the identified one or more k-mers. An ordered list of possible k-measures may contain all or some of the possible k-measures in a predetermined order. Step S220 may include generating an ordered list of possible k-mers, or may not include generating an ordered list of possible k-mers (eg, in situations where a direct comparison to a list is not required to determine the minimizer, as in some examples below).

For example, the general minimizing function min() may define as the k-mer minimizer with the integer minimum value of the two-bit binary encodings of all k-mers in the mutated reading ρ ⁱ of the sequence. In other words, the general minimizing function min() can identify the k-mer that is listed first in a list of k-mers that are ordered by the integer value of their two-bit binary encodings. For example, based on the binary encoding of A: 00, G: 01, C: 10, and T: 11, the general minimizing function can identify in the mutated reading the sequence of 5-mers that is listed first in the ordered list example AAAAA, AAAAG, AAAAC, AAAAT, AAAGA , AAAGG,…, STTTTS, STTTTT, TTTTT. For example, an example of a mutated sequence reading:

contains 5-measures ACGGA, CGGAA, GGAAA, … AGAGC, GAGCC. The 5-mer AAAGC is listed first in the ordered list example above, and the general min() function will determine AAAGC as the minimizer for this mutated sequence read example. It should be understood that this general minimizing function min(⋅) is not required to actually generate an ordered list of possible k-mers in order to define a minimizer for a set of k-mers.

Determining the integer minimum value of the two-bit binary encodings of all k-mers in a mutated read ρ ⁱ of a sequence is just one example of a general minimization function min() that can be applied to a mutated read of a ρ ⁱ sequence to determine a minimizer. Any other general minimizing function min() can be used. For example, it is preferred that the general minimizing function min() randomize the ordering of the integer minimum function. One way to achieve this randomization is to first apply a bitwise logical OR with an arbitrary bit vector to each k-mer contained in the mutated read ρ ⁱ of the sequence, after which an integer minimum function can be used.

In an alternative embodiment, instead of an ordered list of possible k-mers, a predefined set of possible k-mers can be used, and the application of the general minimizing function min() includes identifying one or more k-mers that exist in the predefined set of possible k-mers. One or more minimizers defined for the corresponding mutated reading of the ρ ⁱ sequence may represent a failure of the identified one or more k-mers. The predefined set of possible k-mers may be ordered or unordered. The predefined set of possible k-measures may be a set of k-measures, including only k-measures that are suitable or intended to be used as minimizers. Step S220 of applying the general minimizing function min( ) may include generating a predefined set of possible k-measures.

In a preferred embodiment, in an ordered list of possible k-mers, the k-mers are ordered based on the probability that the k-mers occur in a mutated sequence μ ⁱ and not occur in a non-mutated sequence S ⁱ , i.e. k-mers that are relatively unlikely to occur in a mutated sequence but not in a mutant-free sequence can be listed higher in the ordered list, and k-mers that are relatively unlikely to occur in a mutated sequence but not in a mutated-free sequence , can be listed below in an ordered list. In an alternative preferred embodiment, the predetermined set of possible k-mers contains k-mers that are relatively likely to be present in the mutated sequence but not in the non-mutated sequence, and optionally the set is free of k-mers that are relatively unlikely to be present in the mutated containing sequence. sequence, but not in a mutation-free sequence. Step S220 may include determining which k-mers contained in the set of mutated reads P of the sequence are relatively likely to occur in the mutated sequence but not in the non-mutated sequence, for example, by comparing the number of occurrences (or observations) of the k-mer in the set mutated reads of the P sequence with the number of occurrences of the k-mer in the set of unmutated reads of the R sequence. The step may include counting the number of occurrences of the k-mer in the set of mutated reads P of the sequence and counting the number of occurrences of the k-mer in the set of unmutated reads R of the sequence.

In both preferred embodiments, the overall minimizer function min( ) is chosen to preferably define as one or more minimizers those k-measures that are more likely to occur in the mutated reading of the ρ ⁱ sequence than in the unmutated reading of the r ⁱ sequence. This makes it more likely that each minimizer contains a mutation.

In a more preferred embodiment, the ordered list of possible k-mers contains only those k-mers, i. e. consists only of those k-mers that occur more often in the set of mutated reads of the P sequence than in the set of unmutated reads of the R sequence (or more often in a sequence containing mutations than in a sequence not containing mutations), i.e. k-mers for which the number of occurrences in the set of mutated reads P of the sequence is greater than the number of occurrences in the set of unmutated reads R of the sequence. In an alternative more preferred embodiment, the predefined set of possible k-mers contains only those k-mers, i. e. consists only of those k-mers that occur more frequently in the set of mutated reads of the P sequence than in the set of unmutated reads of the R sequence (or more often in a sequence containing mutations than in a sequence not containing mutations), i.e. k-mers for which the number of occurrences in the set of mutated reads P of the sequence is greater than the number of occurrences in the set of unmutated reads R of the sequence. Preferably, the ordered list of possible k-measures or the predefined set of possible k-measures contains only those k-measures, i. e. consists of those k-mers that occur n or more times in the set of mutated reads of the sequence and occur less than n times in the set of unmutated reads of the sequence, i.e. those k-mers for which the number of occurrences in the set of mutated reads R of the sequence is equal to or greater than n, and the number of occurrences in the set of unmutated reads R of the sequence is less than n. N may be an integer greater than or equal to 1. N may be an integer greater than or equal to 2. Preferably, "n" is equal to 2. In addition, an ordered list of possible k-measures or a predefined set of possible k-measures contains only those k-measures, i.e. consists only of those k-mers that do not occur in the set of unmutated readings of the sequence, i.e. those k-mers for which the number of occurrences in the set of unmutated reads R of the sequence is equal to 0.

For example, an ordered list of possible k-mers or a predefined set of possible k-mers may contain only those k-mers that occur at least twice in the set of k-mers of the set of mutated readings P of the sequence, but do not occur (or are rare) in the set of k-mers of the set of unmutated reads of the R sequence. This ensures that, with high probability, an ordered list of possible k-mers or a predefined set of possible k-mers will include minimizers that contain a mutation present in two or more of the set of mutated P sequence reads. Optionally, k-mers that occur more frequently in the set of mutated sequence reads than in the set of non-mutated sequence reads are relatively likely to occur in the sequence containing the mutation. Optionally, those k-mers that occur n or more times in the set of mutated sequence reads and occur less than n times in the set of unmutated sequence reads are relatively likely to occur in the sequence containing the mutation.

A predefined set of possible k-measures can be generated by constructing a set of mutational minimizers U _M , where U _M contains k-measures, preferably all or substantially all k-measures, for which the number of occurrences or observations in the set of mutated reads P of the sequence is greater than or is n (preferably where n ≥ 2, more preferably where n is 2) and the number of occurrences or observations in the set of unmutated reads P of the sequence is less than n (preferably where n is 0 or 1, more preferably where n is 0) . A set of mutational minimizers U _M can be created by counting the frequency of occurrence of each k-mer in the set of unmutated reads R of the sequence and the set of mutated reads P of the sequence. A set of mutational minimizers U _M can be efficiently computed from a set of unmutated reads R of a sequence and a set of mutated reads P of a sequence using probabilistic data structures such as the counting Bloom filter or the associated cuckoo and coefficient filter methods. An ordered list of possible k-mers can be generated from the entire set of mutational minimizers U _M .

The set of mutational minimizers U _M can be used as a predefined set of possible k-measures. In an alternative embodiment, the set of mutational minimizers U _M can be further processed to obtain a predefined set of possible k-mers. In a preferred embodiment, a subset W _M of the set of mutational minimizers U _M is used as a predefined set of possible k-mers. A subset W _M can be constructed by dividing each mutated reading ρ ⁱ ∈ _P into two or more non-overlapping sections (not necessarily of essentially the same size), e.g. non-overlapping sets of starting positions of k-mers of size L _w , e.g. _w +1…2L _W }, etc. A typical value for _Lw could be 50 using mutated sequence reads of length 150, resulting in the start positions of possible k-mers being divided into 3 groups. Then, for each set of initial positions, the subset W _M can be denoted as follows:

As a result, each of the plurality of mutated P sequence reads can be divided into two or more sections (eg, 3 sections) and a minimizer representing each section can be found. A minimizer is determined by first identifying potential minimizers by the intersection of k-mers in that section of the corresponding mutated reading of the sequence with a set of mutational minimizers U _M , and then applying a general minimizing function to this set to identify one minimizer for each section.

Thus, in a preferred embodiment, step S220 of applying the general minimization function min() to each mutated sequence read comprises:

creation of a set of mutational minimizers U _M which consists of k-mers, preferably all or essentially all k-mers in the set of mutated reads P of the sequence that occur n or more times in the set of mutated reads P of the sequence and occur less than n times in the set of unmutated reads reads R of the sequence, where n is an integer greater than or equal to 2;

optionally creating a subset W _M from a set of mutational minimizers U _M by dividing each of the set of mutated reads P of the sequence into two or more sections, identifying k-mers, preferably all or substantially all of the k-mers in each section of each of the set of mutated reads P of the sequence , which occur in the set of mutational minimizers U _M , and adding to the subset W _M one of the identified k-mers for each section of each of the set of mutated reads P of the sequence, optionally one of the identified k-mers for each section of each of the set of mutated reads P sequences are selected by applying a common minimizing function min(⋅) (eg, an integer minimum detection function or any other known minimizing function) to the identified k-mers of each section of each of the plurality of mutated P sequence reads; And

using a set of mutational minimizers U _M or a subset of the set of mutational minimizers U _M (for example, a subset of W _M ) as a predefined set of possible k-measures, and for each of the set of mutated readings P of the sequence, identifying k-mers, preferably all or essentially all k-mers in the corresponding mutated sequence reading μ ⁱ that occur in a predefined set of possible k-mers, wherein one or more minimizers defined for the corresponding mutated sequence reading are the identified k-mers.

The method 200 further includes step S222 to determine the positions] of one or more appropriate minimizers in each mutated read ρ ^{i of} the sequence. The positions j of each of the minimizers in each respective mutated read ρ ⁱ of the sequence may be stored as an integer bit value in association (eg, in the same location or minimizer group) with the corresponding minimizer.

Step S224: distribution to minimizer groups

In a preferred embodiment, method 200 includes the step S224 of allocating the mutated P sequence reads into one or more minimizer groups. Assigning mutated reads P of the sequence to one or more minimizer groups includes placing an index i characterizing the mutated read ρ ^{i of} the sequence into one or more minimizer groups. Each minimizer group may contain mutated P sequence reads that have a common minimizer, and may not contain mutated P sequence reads that do not have a common minimizer. The step S240 of counting the number of matched and/or mismatched mutations can only be performed on mutated P sequence reads that are in the same minimizer group. This improves the computational efficiency of step S240.

In other words, one or more minimizers can be used as hash keys to collect the reads of the sequence containing the minimizer into a common hash segment (herein referred to as the minimizer group), for example, in preparation for some additional processing (for example, step S240) held on these readings of the sequence.

Each minimizer, which is determined by applying the general minimizing function min() to the mutated P sequence reads, can be used to allocate the mutated P sequence reads to one or more minimizer groups. In one embodiment, each minimizer in an ordered list of possible k-mers, or each minimizer in a predefined set of possible k-mers (for example, each minimizer in the set of mutational minimizers U _M or a subset thereof, such as a subset of W _M ), can be used for distribution purposes. mutated reads of the P sequence into one or more minimizer groups.

Step S224 of distributing the mutated P sequence reads into one or more minimizer groups may include creating one or more minimizer groups. This may include creating one minimizer group for each minimizer defined with the common min() function, or one minimizer group for each minimizer (or k-measure) in a predefined set of possible k-mers U _M , or one minimizer group for each k-measure in subset W _M . Each minimizer group can be implemented as a contiguous block of RAM. Preferably, the minimizer collections are implemented as a file on a computer storage medium (such as a computer disk, such as a spinning magnetic disk or solid state disk), allowing each group to store large amounts of data (which is appropriate in cases of sequence analysis).

Step S224 of allocating mutated sequence reads P to one or more minimizer groups may include storing the mutated sequence read ρ ⁱ or an indicator i characterizing the mutated sequence read ρ ⁱ in the corresponding minimizer group. Step S222 determining positions] of one or more corresponding minimizers in each mutated read ρ ⁱ of the sequence may include storing the position j of the corresponding minimizer in the corresponding minimizer group. In addition, the position α=morphomuts(ρ ⁱ ,V _R ) of one or more mutations in each mutated sequence read ρ ⁱ determined in step S230 of determining the positions α of one or more mutations in each mutated sequence read can be stored in the corresponding minimizer group. Optionally, arbitrary additional values can be stored in the minimizer group, such as the sequence of the mutated read ρ ⁱ of the sequence, quality information regarding sequence accuracy, or other information if these data are useful for subsequent processing. These values associated with each mutated read ρ ⁱ of the sequence can be stored as a tuple in each minimizer group. For convention, the elements of the tuple of the y-th element of the z-th group of the minimizer b _z,y are denoted as b _z,y .i, b _z,y .j, and b _z,y .α. Each mutated reading of the ρ ⁱ sequence can be added to several minimizer groups.

Step S230: mutation positions

The method 200 includes the step S230 of determining the a positions of one or more mutations in each mutated read ρ ⁱ of the sequence. Step S230 of determining the α positions of one or more mutations in each mutated sequence read ρ ⁱ can be performed using methods that do not involve alignments.

Step S230 of determining the α positions of one or more mutations in each mutated sequence read ρ ⁱ may include obtaining a set of seeded unmutated k-mers V _R , i. a set of k-mers of an unmutated reading R of a sequence to which one or more seed patterns ψ have been applied. Obtaining a set of seeded unmutated k-mers V _R may include creating or generating a set of seeded masked unmutated k-mers V _R . A set of seeded, unmutated k-mers V _R can be obtained or generated by applying each of the one or more seed patterns to each k-mer in a non-mutated sequence, eg, to each k-mer in unmutated reads of the sequence. Applying the seed pattern to the k-measure may include determining the result of a bitwise logical AND to the seed pattern and the (two-bit-encoded) k-measure. Applying a seed pattern to a k-mer results in a seed-masked k-mer. The set of seed-masked unmutated k-mers V _R can be denoted as

those. a set of seed-masked unmutated k-mers V _R is created by applying each of one or more seed patterns ψ of the seed family ψ to each k-measure k(r _j ⁱ ) for all (or substantially all) positions of the j k-mer (i.e., e. from 1 to in each unmutated read r ⁱ for all (or essentially all) unmutated reads r ⁱ in the set of unmutated reads R of the sequence.

The seed pattern can be used to modify the way k-mers are compared to each other. A seed pattern is defined as a set of positions (ie, nucleotides) within two k-mers that must be identical in both k-mers in order for the seeded k-mers to be considered the same. The seed pattern may contain mask positions and non-mask positions. Applying seed patterns to a k-mer creates a seeded k-mer, in which the seeded k-mer positions corresponding to the masked positions of the corresponding seed pattern are ignored in any additional processing (such as comparisons), while the seeded k-mer positions are corresponding to non-masking positions of the corresponding seed pattern are not ignored by any additional processing (such as comparisons). For example, the seed pattern {1, 2, 4, 6, 7} requires that the first, second, fourth, sixth, and seventh positions (or nucleotides) in the two compared k-mers k(S _i ) and k(S _j ) be are identical to be considered the same (for k=7). The third and fifth positions in the two k-mers can be arbitrary nucleotides. This means that the third and fifth positions in the two seeded k-mers are masked by the seeding pattern.

The one or more seed patterns may optionally be one or more transition seed patterns. This is advantageous in particular when the mutated sequence M contains transition-type mutations compared to the non-mutated sequence S, ie. each of the plurality of mutated P sequence reads contains one or more transition type mutations.

A transition seed pattern is a specialized type of seed pattern where positions are divided into three classes instead of just two: each position must (1) match exactly, or (2) both must be purine or pyrimidine, or (3) each position must be any of four nucleotides to match. Transition seed patterns are particularly preferred when the mutated sequence contains transition-type mutations. When implemented on a computer using the two-bit nucleotide encoding proposed above, the position where an exact match is needed can be implemented as a bitmask of 11, while the position where only transition-type mutations are allowed is designated as 10, and the position at which any nucleotide is allowed is denoted as 00. The seed pattern {1, 2, 4, 6, 7} can be written as the bitmask 11110011001111. The transition seed pattern {1, 2, 4, 6, 7} can be written as the bitmask 11111011101111. Two k-measures can be evaluated for a match by calculating for each of them the result of a bitwise logical AND for the bitmask and two-bit encoding of the k-measure, and then checking the identity of the two seeded masked k-measures obtained. For convenience, a function that applies a seed pattern to the k-measure k(S _i ) by a bitwise logical AND will be denoted as a function ψ(k(S _i )).

In one embodiment, one or more seed patterns are selected such that the probability of obtaining identical seed masked k-mers when at least one of the one or more seed patterns is applied to any k-mer of a plurality of mutated P sequence reads (or sequences containing mutations) is ) and the corresponding k-mer of the plurality of unmutated R sequence reads (or mutated sequences) is greater than 90%, preferably greater than 95%, further preferably greater than 98%, most preferably greater than 99%.

One or more seed patterns may constitute a family of seed patterns ψ. A seed pattern family ψ is a set of two or more seed patterns that, when used together, are capable of identifying matches among k-mers at a particular percentage of nucleotide identity with a high probability, e.g., greater than 90% probability, preferably greater than 95%, further preferably greater than 98 %, most preferably more than 99%. A family of seed patterns ψ is denoted as a set of n different functions for applying seed patterns ψ ₁ …ψ _n ∈ ψ. The seed pattern weight w(ψ) is defined as the number of seed positions that must be identical for two k-measures to be considered the same, where w(ψ)≤k.

For each mutated sequence read, step S230 of determining the positions α of one or more mutations in each mutated sequence read ρ ⁱ may include applying each of the one or more seed patterns ψ _i to k-mers (optionally to each k-mer) in the corresponding mutated read ρ ⁱ sequences to obtain a set of seeded mutated k-mers. The positions of one or more mutations can be determined by identifying one or more positions in the mutated read ρ ⁱ of the sequence masked by all seed patterns that correspond to seeded mutated k-mers from the set of seeded mutated k-mers that occur in the set of primed masked unmutated k-mers. k-mers V _R . This means that positions that are not mutations in the mutated read ρ ⁱ of the sequence can be identified as one or more positions in the mutated read ρ ^{i of} the sequence that are not masked by any of the seed patterns that correspond to the seed masked mutated k-mers from sets of seeded mutated k-mers that exist in the set of seeded unmutated k-mers V _R .

For example, the positions a of one or more mutations of each mutated read ρ ⁱ of the sequence can be determined as follows:

• creating a bit vector a of length 2L and setting the initial values of the bit vector α to 0;

• creating a bit vector b of length 2k and setting all initial values of the bit vector b to 1;

• for every ψ ∈ Ψ and for every position j in reading between 1 and - calculation of ψ(k(ρ _j ⁱ )). If ψ(k(p _j ⁱ )) ∈ V _R , then the assignment α←α | (ψ(b)>>2j), where the operator | stands for the bitwise logical OR operator, and the >> operator stands for the right shift operator. Querying whether a set belongs to a V _R can be implemented either exactly, using something like a hash table, or approximately, using a high-performance probabilistic data structure such as a Bloom filter, a coefficient filter, or a similar approach.

• Optionally, for ease of further processing, converting the bit vector a from the length of the binary two-bit encoding of the mutated sequence read to the length of the mutated sequence read itself, by removing odd positions, eg α←{α2,α4,α6,…α2L}.

• Optionally, for ease of further processing, apply a logical NOT operation to the bits, so that a value of 1 represents positions for which no match was found with the seed.

The result of the above procedure will be a bit vector α, where each position containing 1, with a high probability corresponds to the position of the mutation. For convention, the function that calculates the bit vector α for the mutated read ρ ⁱ of the sequence is denoted as α=morphomuts(ρ ⁱ ,Y _R ).

On FIG. 3 shows an example illustrating how a bit vector α can be obtained for an example of a mutated read ρ=ACGCAAAGCGCTACGAGCGACTGATATT sequence using a single seed pattern ψ=1110110011. The 4th, 8th, 11th, 12th, and 16th positions of the mutated reading of the ρ sequence correspond to mutations in the mutated reading of the ρ sequence, i.e. the nucleotides at these positions in a mutation-free sequence will be different. In practice, a mutated reading ρ of a sequence can be encoded in a two-bit binary format, and each position of the seed pattern ψ can cover two bits (i.e., each 1 in the seed pattern ψ will be realized as two binary 1s, and each 0 in the seed pattern will be implemented as binary 00 or binary 10). A set of seeded, unmutated V _R k-mers was generated earlier in this example.

As shown in the example shown in FIG. 3, seeding patterns are applied to each k-mer in the mutated reading of the ρ sequence, thus generating one seeded k-mer for each k-mer in the mutated reading of the ρ sequence. It is then checked whether the seeded k-mers exist in the set of seeded non-mutated k-mers V _R . In the example shown, all of the 1st, 5th, 13th, 17th, 18th, and 19th seeded k-mers occur in the set of seeded non-mutated k-mers V _R . These seeded k-mers do not contain mutation positions not masked by the seeding pattern.

The 1st, 5th, 13th, 17th, 18th, and 19th seeded k-measures are then used to identify positions that are masked by all seed patterns corresponding to those seeded k-measures. The 4th position of the mutated read of the ρ sequence is masked by all these seed patterns, indicating that the 4th position of the mutated read of the ρ sequence is ignored when processing the 13th, 17th, 18th, and 19th seeded k-mers , i.e. The 4th position of the mutated reading of the ρ sequence is masked by a seed pattern for the 13th, 17th, 18th, and 19th seeded k-mers. None of these seed patterns mask the 4th position of the mutated reading of the ρ sequence. Thus, the 4th position of the mutated reading of the ρ sequence is identified as the mutation position. On the contrary, although the 7th position of the mutated reading of the ρ sequence is masked by all seed patterns corresponding to the 1st, 13th, 17th, 18th, and 19th seeded k-mers, this is the 7th position of the mutated reading of ρ sequence is not masked by a seed pattern corresponding to the 5th primed k-mer. Thus, the 7th position of the mutated reading of the ρ sequence is not identified as a mutation position. Instead, the 7th position of the mutated reading of the ρ sequence is identified as a position that is not a mutation.

Essentially all seed patterns corresponding to the 1st, 5th, 13th, 17th, 18th and 19th seeded k-mers are combined using a logical OR. The bits of the resulting bit vector may be mirrored (eg, using a logical NOT operation) to obtain the positions of the mutations in the mutated read ρ of the sequence as the bit vector α.

Alternative implementation of step 230 using a reference assembly

In the embodiment described above, step 230 of determining the positions α of one or more mutations in each mutated sequence read ρ ⁱ is performed using a plurality of mutated sequence reads P and a plurality of unmutated R sequence reads based on the application of seed patterns to each mutated sequence read ρ ⁱ .

In large and complex genomes, such as the human genome, a significant portion of the genome is made up of repetitive sequences. For example, more than half of the human genome is thought to be part of repetitive sequences. These repeat sequences are classified into "families" of similar repeat sequences. The most abundant in the human genome is the Alu family of short dispersed nuclear elements (SINEs), which are about 300 nt long and present in approximately 1 million copies. Another common family is the L1 family of long dispersed nuclear elements (LINEs) with element sizes ranging from 1 to 6.5 kb. and with a number of copies of about 10,000.

Different copies of repeating sequences in the genome may not be identical, for example, they contain differences in single bases. Because of the biology of the mutation, these differences are often transition type differences. In some situations, these differences may look similar to differences due to the introduction of mutations between a set of mutated P sequence reads and a set of unmutated R sequence reads. This is especially true for some polymerase mutagenesis approaches used to introduce mutations in at least one template target nucleic acid molecule in the context of generating multiple mutated P sequence reads, as transitional type mutations are often introduced.

As a result, the set of unmutated reads of the R sequence may contain a large number of k-mers that differ from each other only by a certain number of differences in the type of transition. Therefore, the set of mutated P sequence reads may include one or more k-mers that are identical to the k-mers of the set of unmutated R sequence reads, despite the presence of mutations compared to the unmutated R sequence reads. In some situations, it is possible that natural differences between different copies of the repeat sequence in different unmutated reads of the r ⁱ sequence will partially "mask" mutations introduced into multiple mutated reads of the P sequence. This is especially pronounced for SINEs from the Alu families.

Thus, it would be preferable if, in such situations, an embodiment of the method were provided to better distinguish intentionally introduced mutations from natural differences between copies of repeating sequences.

A first approach to improving the ability of the method to distinguish intentionally introduced mutations from natural differences between copies of repeat sequences is to use seed patterns with a much higher mass, so that mutated primed-masked k-mers are more likely to include one or more positions containing a difference that distinguishes copies repeating sequence. In one embodiment that employs the first approach, the mass w(ψ) of each seed pattern ψ is in the range of 50 to 100, preferably in the range of 70 to 90. For the human genome, a mass of approximately 80 will suffice for the first approach.

However, the first approach may not be ideal in all cases. A seed pattern with a mass of 80 will be very long, probably longer than the typical length of a mutated read ρ ⁱ of the sequence. In addition, the family size Ψ of seed patterns required to provide high sensitivity can become very large, which requires significant additional computing resources to process all seed patterns. Finally, the probability of insertion-deletion error coverage by the seeding pattern will increase, and additional algorithmic complexity will be required to adapt to the possibility of insertion-deletion errors. Thus, in some circumstances this first approach may not be preferred.

A second approach to improving the ability of a method to distinguish intentionally introduced mutations from natural differences between copies of repeat sequences is to use an approach based on aligning (or matching) multiple mutated P sequence reads with a reference assembly (or reference genome). The reference assembly can either be independently generated, such as the human hg38 genome generated by the Reference Genome Consortium (GRC), or can be a de-novo assembly based on multiple unmutated R sequence reads. For one or more mutated sequence reads in the second approach, the step of determining the positions of one or more mutations in each mutated sequence read comprises aligning the corresponding mutated sequence reads with a reference assembly.

This approach may be particularly suitable when the mutated ρ ⁱ sequence reads are paired-terminal sequence reads. An advantage of aligning mutated paired-ended sequence reads with a reference assembly, particularly as applied to SINE repeats, is that the fragment size in the shotgun short-read library is generally larger than the length of repeated sequences. The typical fragment size for paired-end sequencing is 400,600 bp, with about 150 bp. sequenced from each end of the fragment. Thus, if one matched-end sequence read of a pair of matched-end sequence reads falls on a repeat sequence, then another of the matched-end sequence reads of a pair of matched-end sequence reads is likely to fall on a unique sequence outside of the repeat sequence. Thus, a standard paired-end alignment program (eg, a Burrows-Wheeler alignment tool such as BWA-MEM) is able to reliably match a pair of paired-ended sequence reads to the correct location in the reference assembly, including the correct copy of the repeating sequence. The positions of any differences between the aligned mutated sequence reads ρ ⁱ and the reference assembly can then be recorded and stored in a bitmap α similar to that obtained using the seed pattern approach for each mutated sequence read ρ ⁱ . Thus, determining the positions of one or more mutations in the corresponding mutated sequence reading is accomplished by identifying, in the corresponding mutated sequence reading, position differences between the corresponding mutated sequence reading and the reference assembly.

However, aligning multiple mutated P sequence reads with a reference assembly may not be ideal in some situations, as any given template target nucleic acid molecule will typically have regions that are not represented in the reference assembly. Therefore, it is not possible to align the mutated sequence reads ρ ⁱ with those regions that are not present in the reference assembly and obtain a bit vector α from the differences between the aligned mutated sequence reads ρ ⁱ and the reference assembly. In addition, regions that are not represented in the reference assembly are often of clinical interest because they represent structural insertion variants relative to the reference assembly. In addition to the large insert regions, any mutations occurring in small inserts relative to the reference assembly will also be missed by an approach based on aligning multiple mutated P sequence reads with the reference assembly.

Thus, a third hybrid approach to improving the ability of the method to distinguish intentionally introduced mutations from natural differences between copies of repeat sequences is to combine the approach based on aligning multiple mutated reads of the P sequence with a reference assembly and the approach based on applying seed patterns to each mutated read ρ ⁱ sequences. This third approach can be used as an alternative implementation of step 230 of the present method.

In a third approach, the position of one or more mutations in each mutated sequence read is determined using both approaches: by aligning a set of mutated sequence reads P with a reference assembly, and by applying seed patterns to each mutated sequence read ρ ⁱ . If the position in the corresponding mutated sequence reading is aligned with the reference assembly, then the position in the corresponding mutated sequence reading is defined as the position of the mutation in the corresponding mutated sequence reading if the position in the corresponding mutated sequence reading is the position at which the corresponding mutated sequence reading differs from the reference assembly . If the position in the corresponding mutated sequence reading is not aligned with the reference assembly, then the position in the corresponding mutated sequence reading is defined as the position of the mutation in the corresponding mutated sequence reading if the position in the corresponding mutated sequence reading is a position that is masked by all seed patterns that match the masked seeded mutated k-mers from the set of seeded mutated k-mers that occur in the set of seeded unmutated k-mers.

To achieve this, a bit vector a of the type described above is independently obtained by both equalization-based and seeding-based approaches. The bit vector from the seed pattern approach to each mutated read ρ ⁱ of the sequence is denoted α _mmd , and the bit vector based from the approach based on applying alignment of the set of mutated reads P of the sequence with the reference assembly is denoted α _map. An additional alignment mask bit vector, denoted α _amask , is also constructed, which records those positions of each mutated sequence read that successfully align with the reference assembly. The equalization mask bit vector α _amask will have 1 at each position that is successfully aligned and 0 at positions that were not successfully aligned with the reference assembly.

The resulting hybrid bit vector α _hybrid is then constructed, which combines the bit vector from the seeding approach to each mutated sequence read ρ ⁱ , α _mmd , and the bit vector from the approach based on aligning the set of mutated reads P of the sequence with the reference assembly, α _map , in the following way:

α _hybrid= α _map | ( _α.mmd & ~ _αamask )

Where | denotes the bitwise logical OR operator, & denotes the bitwise logical AND operator, and ~ denotes the bitwise NOT operator.

Thus, a third approach uses mutation positions determined by alignment with a reference assembly at positions in the mutated sequence read where the alignment was successful, and mutation positions determined by seeding patterns at all other positions. This provides the advantage of being able to include a high quality reference assembly in the analysis while processing all types of inserts relative to the reference assembly. Alignment to an independent, high quality reference assembly, such as the human genome reference, can be much more accurate than alignment to an assembly of de novo short reads. Using mutation positions determined by alignment with a master assembly can provide more accurate estimates of mutation positions, especially in regions of repetitive sequences, while a seed pattern-based method without alignment can identify mutation positions in regions that are not present in the master assembly. The latter can occur without the need to compute the assembly, which is a computationally demanding task. Thus, the hybrid approach provides improved accuracy in identifying mutation positions and computational efficiency relative to using either approach alone.

It is also possible to "augment" the reference assembly with variants and locally assembled regions from a particular template target nucleic acid molecule to obtain an assembly graph specific to that template target nucleic acid molecule. A bit vector from the approach based on aligning a set of mutated reads P of a sequence with a reference assembly (denoted α _map ) can be obtained by aligning the mutated reads of a sequence with an enlarged assembly graph, and then combining with the approach based on applying seed patterns to each mutated read ρ ⁱ sequences for any regions of the template target nucleic acid molecule that remain difficult to align for technical or other reasons.

Steps S240 and S240: determining a score correlated with the probability that two mutated sequence reads are from the same mutated sequence

The method 200 includes step S240, in which at least two mutated reads of a sequence with a common minimizer are counted by the number of matched and/or mismatched mutations when the respective minimizers are aligned.

This can be achieved by first determining the difference in position j of the minimizer determined in step S222 for each of the two mutated sequence reads. For example, the difference in position j of the minimizer for each of two mutated reads ρ ^a and ρ c ^of the sequence stored in the minimizer group as a=b _z,y .i and c=b _z,x .i can be defined as d=b _{z ,y} .jb _z , _x .j.

Counting the number of mutations with matching positions may include determining the intersection size of the sets for the mutation positions determined in step S230 when the mutation positions determined for one of the two mutated readings ρ^X and ρ^at sequences stored as b_z,y and b_z,x, have a right shift d. For example, for two mutated readings ρ^X and ρ^at sequences stored as b_zy and b_z,x, the number of mutations with matching positions can be determined as follows:

λ _x,y =|Ω(b _z,x .α) ∩ (Ω(b _z,y .α) - d)|, where Ω(α) is defined as the set of position indicators in α that are non-zero (i.e., i.e., the set of mutation positions in the corresponding mutated reading ρ ⁱ of the sequence) and where Ω(b _z,y .α) - d is understood as element-wise subtraction of d from Ω(b _z,y .α). Set intersection can be efficiently implemented on a computer using CPU instructions such as bitwise shift and popcount.

Counting the number of mutations with mismatched positions may include determining the size of the symmetrical set difference for the mutation positions determined in step S230, when the mutation positions determined for one of the two mutated sequence reads ρ ^x and ρ ^y , stored as b _{z, y} and b _{z, x} , have a right shift d. For example, for two mutated reads ρ ^x and ρ y ^of the sequence, stored as b _{z, y} and b _{z, x} , the number of mutations with mismatches can be determined as follows:

δ _x,y =|(Ω(b _z,x .α) \ (Ω(b _z,y .α) - d))∪ ((Ω(b _z,y .α) - d) \ Ω(b _z,x .α))|.

Step S242 of determining an indicator correlated with the probability that at least two mutated sequence reads are from the same mutated sequence may be based on the number of mutations with a matching λ _x,y position and/or a mismatching δ _{x position, U} . In one embodiment, a score correlated with the probability that at least two mutated sequence reads are from the same mutated sequence corresponds to the number of mutations with the same position λ _x,y . The higher the number of mutations with matching positions λ _x,y , the higher the probability that at least two mutated sequence reads come from the same sequence containing the mutation. In an alternative embodiment, the score correlated with the probability that at least two mutated sequence reads are from the same mutated sequence corresponds to the number of mismatched mutations δx _,y . The lower the number of mutations with mismatched positions δ _x,y , the higher the probability that at least two mutated sequence reads come from the same sequence containing the mutation.

In a preferred embodiment, said metric, correlated with the probability that at least two mutated sequence reads are from the same mutated sequence, is one of: i) a probability density that at least two mutated sequence reads are from the same mutated sequence, and ii) a scoring function that is correlated with a probability density that at least two mutated sequence reads are from the same mutated sequence.

For example, the number of mutations with matching positions λ _x,y and the number of mutations with mismatches δ _x,y can be used to calculate in the probability density model that two reads are from the same mutation-containing sequence M, or an estimator that is correlated with this probability density. One such evaluation function is ω _a,c =Δ(λ _x,y ) - Δ(δ _x,y ), where Δ(n)=(0.5n)(n+1) for a=b _{z, x} .i and c=b _z,y .i. Thus, ω _{a, c} is the score or weight of the link between two mutated sequence reads ρ ^a and ρ ^c . A collection of such links may be obtained for all pairs of mutated reads ρ ⁱ of the sequence in the corresponding minimizer group b _z , or if there are a large number of entries in the b _z minimizer group, calculation or reporting of the links may be limited to randomly selected pairs of records in the minimizer group b _z .

Step S300: Sequence Assembly or Sequence Reconstruction

The method 10 may further include the step S300 of assembling or reconstructing a sequence, or at least a portion of a sequence, such as a mutated sequence or a non-mutated sequence. The assembled or reconstructed sequence may be a mutated or non-mutated sequence.

Method 200, such as sequence reconstruction or assembly step S300, may include generating an undirected weighted graph from a plurality of mutated sequence reads. An undirected weighted graph contains nodes corresponding to a set of mutated sequence reads. For example, each node may correspond to a corresponding mutated sequence read in the sense that it is represented by the read pointer i of the corresponding mutated sequence read or by the sequence of the corresponding mutated sequence read. Edges between nodes are associated with respective edge weights, where the weight of each edge can be determined based on the number of match and/or mismatch mutations determined for the two mutated sequence reads corresponding to the two nodes associated with the respective edge. The weight of each edge may correspond to a score correlated with the probability that at least two mutated sequence reads (i.e., two mutated sequence reads corresponding to nodes associated with the edge associated with the edge weight) are from the same containing mutation. sequences. Thus, the weight of an edge connecting two mutated sequence reads (nodes) represents the probability that the two mutated sequence reads were derived from the same mutated sequence, or some other arbitrary function that is correlated with that probability.

An undirected weighted graph can be constructed by processing each of the minimizer groups sequentially or in parallel, thus computing the edges between the mutated sequence reads in each minimizer group. The edge weight may be an estimate function ω _a,c .

An undirected weighted graph including the edge weights ω _{a, c} can then be used to process SAM data (e.g., mutated sequence reads), for example using any known or unknown techniques for applying such an undirected weighted graph to assemble a sequence. Assembling a sequence from an undirected weighted graph may include, for example, generating clusters of mutated sequence reads and assembling the mutated sequence reads in each cluster to reconstruct a template corresponding to at least a portion of the mutated sequence.

For example, method 200 or step S300 of reconstructing or assembling at least a portion of a sequence may include performing a graph clustering operation on an undirected weighted graph, thereby obtaining clusters of mutated sequence reads that are expected to come from the same containing mutation. sequences. Graph clustering can be done using any standard streaming graph clustering algorithm, such as Markov Clustering (MCL) or Infomap. In an alternative embodiment, the edges of an undirected weighted graph can be filtered by some minimum threshold weight, and then to represent mutated sequence reads can be the connected components of the graph are taken, which come from the same sequence containing mutations.

Step S300 of reconstructing or assembling at least a portion of the sequence may further comprise reconstructing at least a portion of the mutated sequence by assembling the mutated sequence reads into clusters. For example, mutated sequence reads in clusters can be processed by de novo standard assembly techniques to reconstruct the mutated sequence. Such de novo assembly methods include, for example, the IDBA-UD algorithm from "IDBA-UD: a de novo assembler for single-cell and metagenomic sequencing data with highly uneven depth", Peng Y et al., Bioinformatics. 2012 Junl;28(11):1420-8. doi: 10.1093/bioinformatics/bts174. Epub 2012 Apr 11, or the SPAdes method from SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing", Benkevich A et al., J Comput Biol. 2012 May; 19(5): 455-477, or A5-miseq method from "A5-miseq: an updated pipeline to assemble microbial genomes from Illumina MiSeq data", Coil D et al., Bioinformatics 2015 Feb 15;31(4):587-9 doi: 10.1093/bioinformatics /btu661.Epub 2014 Oct 22.

Step S300 of reconstructing or assembling the sequence may further comprise reconstructing at least a portion of the mutated sequence using error correction on the reconstructed portion of the mutated sequence, i.e. by deriving the most likely mutant-free sequence from the reconstructed portion of the mutated sequence using a plurality of unmutated sequence reads. Techniques for such error correction include, for example, the FMLRC method from "FMLRC: Hybrid long read error correction using an FM-index", Jeremy R. Wang et al., BMC Bioinformatics volume 19, Article number: 50 (2018). For example, a mutated sequence can be subjected to error correction using unmutated sequence reads to remove introduced mutations, thereby reconstructing portions of the mutated sequence. Error correction may include, for example, determining the possible sets of edits containing mutations of the sequence that would be required to convert the containing mutations to a sequence that does not contain mutations, compatible with unmutated reads of the sequence, determining the set of edits that have the smallest size (i.e., containing the smallest edits ) of the possible edit sets, and applying the determined set of edits having the smallest size to the mutated sequence to obtain a probable estimate of the mutation-free sequence. Mutation-free sequence parts can then be assembled using standard de novo long read assembly tools such as Canu, Flye, or PEREGRINE, or in combination with short reads in R using a tool such as Unicycler or MaSuRCA, with assembly such sequence that does not contain mutations.

Sample Pool Processing

When processing batches of samples containing multiple samples, you can enter sample barcodes as predefined marker sequences for each sample. If a user of method 200 wishes to use the method on multiple samples, where each sample contains one or more mutated template target nucleic acid molecules, one option is to process each sample (e.g., mutant and/or fragment) separately in the laboratory, and then inserting sample barcodes only at the final stage before sequencing. Another alternative is to insert sample barcodes only at the ends of the template target nucleic acid molecules, in which case it is possible to combine all barcoded template target nucleic acid molecules at the beginning of the sample preparation process, thus significantly reducing reagent costs and labor costs (thus called the early-stage pooling approach). Thus, sample preparation may include inserting appropriate sample barcodes into the ends of template target nucleic acid molecules in each sample such that each sample contains template target nucleic acid molecules having a different sample barcode relative to template target nucleic acid molecules in other samples. . Sample preparation may further include pooling samples to form a sample pool, introducing mutations and optionally fragmenting template target nucleic acid molecules in the sample pool, and sequencing portions of the mutated template target nucleic acid molecules in the sample pool.

However, the early sample pooling approach creates an additional processing problem because the resulting set of mutated P sequence reads contains an unlabeled mixture of mutated P sequence reads from many different samples. Samples can be processed individually to construct unmutated R sequence reads, in which case the unmutated R sequence reads contain multiple sets of unmutated R reads¹…R^ζ sequence, where ζ is the number of samples processed in the batch. Each set of unmutated sequence reads may be associated with a corresponding pattern. Method 200 may include receiving an unmutated read R of a sequence in a plurality of sets of unmutated reads R¹…R^ζ sequences, with each set of unmutated reads R¹…R^ζ sequence is associated with the corresponding one or more samples.

Thus, each of the plurality of mutated sequence reads may be a subsequence of the mutated sequence associated with one of the plurality of patterns. Each of the plurality of unmutated sequence reads may correspond to a subsequence of the non-mutated sequence associated with one of the plurality of patterns. Each mutated sequence may contain mutations compared to the corresponding non-mutated sequence. Obtaining a set of unmutated primed masked k-mers may include obtaining a corresponding set of unmutated primed masked k-mers for each sample.

A simple approach to processing data from ζ samples would be to apply method 200 once for each of the ζ samples. An alternative approach is to extend the method 200 such that all ζ samples can be processed at the same time. This can be achieved as described below.

Method 200 (e.g., step S230) may include generating a set of unmutated sample bitvectors, where for a corresponding k-mer in the set of unmutated primed masked k-mers V _R , each wild-pattern bitvector determines in which of the set of samples the corresponding k-mer occurs. (or occurs at least x times, where x is an integer greater than or equal to 1) and in which of the set of samples the corresponding k-mer does not occur (or occurs less than x times). The set of unmutated primed k-mers V _R can be generated from the set of unmutated k-mers in the manner already described above. The set of unmutated k-mers can be defined as the union of all k-mers in each of the set of samples R ¹ …R ^ζ , i.e. the set of unmutated k-mers R can be defined as R=∪R ¹ …R ^ζ ⋅

For example, method 200 may include defining a surjection V _R on a collection of bit vectors containing binary indicators of the presence of seeded k-mers in each sample. Each bitvector can have a 1 at position i if the i-th sample in the sample set contains k-mers (or contains them at least X times), otherwise it has 0 at position i. In a software implementation, the surjection can be stored using an unordered map data structure, such as a hash map, or an approximate membership query structure, such as a counting filter with coefficients. The surjection can be denoted as Z: V _R →ν where ν is the space of bit vectors of length ζ⋅

Step S230 of determining the positions of one or more mutations in each mutated sequence read can be extended to construct a bit vector a simultaneously for multiple samples. For each mutated sequence read, and for each set and/or each combination of sets of unmutated primed masked k-mers, determining the positions of one or more mutations may include identifying one or more positions in the mutated sequence read that are masked by all seed patterns corresponding to the mutated primed masked k -mers from a set of mutated primed masked k-mers that occur in the corresponding set or combination of sets of unmutated primed masked k-mers, and associating the identified one or more positions with one or more patterns associated with the corresponding set or combination of sets of unmutated primed k-mers -mers. This can be achieved, for example, with the multi-sample oriented variant morphomutsMS(ρ ⁱ ,V _R ) of the function morphomuts(ρ ⁱ ,V _R ), which includes the following steps:

1. Initialization of the set A of bit vectors a with one initial bit vector a ₀ of length 2, containing only bits 0; initialization of bit vector b of length 2k containing only bits 1; initialization of the set C of bit vectors with one initial element with ₀ of length ζ, containing only bits 1; mapping initialization G:

2. For each position j in reading between

a. For each seed pattern determine ψ(k(ρ _j ⁱ ))

b. If , perform the following steps:

i. For each element with (i.e. for each ) calculate d←c∧Z(ψ(k(ρ _j ⁱ ));

ii. If d contains only 0 bits, then return to 2b.i to process the next C element (or, if there are no more, to process the next seed pattern, or the next position j), otherwise continue with 2b.iii;

iii. assign α←Г(с) | (ψ(b)>>2j), where | means bitwise logical OR and >> means bitwise right shift operator, and remove c from C;

iv. Add d to C and α to A and determine the juxtaposition d→a in D;

v. If not equal to zero, then add to C and G(c) to A and define the mapping

vi. return to 2b.i to process the next element c from C. If there is no more c in C, return to 2a to process the next seed pattern ψ⋅ If there is no more ψ in Ψ, return to 2 to process the next position j. Otherwise, continue with 3, otherwise:

3. Convert the bit vectors to A using the transformation applied to create α in the morphomuts(⋅) function; And

4. Return C, A, and match D as the result of the function.

Optionally, if there are too few matching positions in bitvector A (e.g., less than a given number y, where y is an integer greater than or equal to 1, preferably greater than or equal to 2, 3, 4, or 5), corresponding entries in C and A may be discarded. This is advantageous because such entries may appear due to random similarities between input patterns, and thus the resulting bit vectors are the result of false matches with a bad pattern. By discarding these positions before further processing, unnecessary calculations can be avoided. The method 200 may include comparing the number of identified positions to a predetermined number y, where y is an integer greater than or equal to 1, preferably greater than or equal to 2, and if the number of identified positions is less than the predetermined number y, discarding (or ignoring if further processing) of the identified one or more positions and the association of the identified one or more positions with one or more samples.

Tuples that are stored in minimizer groups can be extended to include information about the bit vector of the pattern in C. In particular, the stored tuple can be the read score i, the position j of the minimizer in the mutated sequence read, and c and α, where c represents is a bit vector of samples, and α is a bit vector of mutations calculated by the function morphomutsMS(p ⁱ ,V _R ).

Subsequently, when processing groups of minimizers to calculate edge weights, each edge weight can be annotated with the corresponding set of samples. If the bitwise AND of the pattern bitvectors associated with a pair of mutated sequence reads is zero, then the corresponding edge may be discarded. If the edge score is less than a predefined threshold score value, the edge may be discarded. When there are multiple possible edges between a pair of mutated sequence reads, it becomes possible to store only the highest edge weight, and the associated set of bitvectors for that edge can be computed by bitwise ANDing the sample bitvectors. This approach has the advantage that natural sequence variation in different samples can be distinguished from mutations introduced during sample processing. The step S240 of counting the number of matched and/or mismatched mutations in aligning the minimizers of the two mutated sequence reads may be performed on any pair of one or more mutation positions identified for the two mutated sequence reads only if there is overlap in the patterns associated with the respective pair. from one or more mutation positions identified for two mutated sequence reads, i.e. only if a pair of one or more mutation positions identified for two mutated sequence reads is associated with at least one common pattern.

If only the ends of the mutated template target nucleic acid molecules contain the pattern tag, then some of the many mutated P sequence reads will carry the pattern tag. In particular, mutated sequence reads resulting from sequencing the ends of mutated target nucleic acid template molecules will carry the pattern label. After clustering the mutated sequence reads, it becomes possible to associate patterns with read clusters simply by evaluating the presence of pattern-tagged reads in each cluster. When only one pattern label occurs in a cluster, the association with the pattern is simple and unambiguous. Performing graph clustering of an undirected weighted graph may include associating with each cluster of mutated sequence reads a pattern label contained in at least one of the mutated sequence reads in the corresponding cluster.

Occasionally, multiple sample marks may occur in the same cluster, either due to noise or an error in the sequencing or data analysis procedures. In this case, it may still be possible to reliably bind to a sample if there is a large excess of one sample label over other labels. In cases where a one-to-one association is not possible, it may still be possible to remove the ambiguity for a sample by applying a semi-supervised graph decomposition procedure that decomposes a cluster of multiple samples into a number of smaller clusters, one cluster per sample label. Even if a cluster contains no reads bearing a pattern label, it may still be possible to associate a cluster with a pattern if most of the pattern masks associated with inter-read joins point to the same pattern. Performing graph clustering of an undirected weighted graph may include identifying, in each cluster of mutated sequence reads, one or more pattern labels contained in the mutated sequence reads in the corresponding cluster of mutated sequence reads. Each cluster of mutated sequence reads can be associated with a pattern label that occurs most frequently in the mutated sequence reads in the corresponding cluster. Optionally, if two or more different pattern marks are identified in a cluster of mutated sequence reads, the cluster of mutated sequence reads may be divided into two or more clusters, with each of the two or more clusters associated with a corresponding one of two or more different pattern marks and containing different sequences of mutated sequence reads.

Sample preparation and sequencing

Method 10 for determining at least a portion of the sequence of at least one template target nucleic acid molecule may include sequencing 100 regions of at least one template target nucleic acid molecule containing mutations to obtain multiple mutated sequence reads. The method 10 for determining at least a portion of the sequence of at least one template target nucleic acid molecule may further comprise performing a method 200 for determining a score correlated with the probability that two mutated sequence reads are from the same mutated sequence, based on the set mutated sequence reads obtained by sequencing 100.

The sequencing step may include:

a) providing a pair of samples, each sample containing at least one template target nucleic acid molecule;

(b) sequencing regions of at least one template target nucleic acid molecule in a first sample of a pair of samples to obtain multiple unmutated sequence reads;

(c) introducing mutations into at least one template target nucleic acid molecule in a second sample from a pair of samples to obtain at least one mutated template target nucleic acid molecule;

(d) sequencing regions of at least one mutated template target nucleic acid molecule to obtain multiple mutated sequence reads.

In a preferred embodiment, the step of introducing mutations comprises introducing transition-type mutations into at least one template target nucleic acid molecule in the second of a pair of samples.

The sequencing step may include:

(a) providing a plurality of sample pairs, each sample pair comprising a first sample and a second sample, each sample containing at least one template target nucleic acid molecule;

(b) inserting a sample barcode into at least one template target nucleic acid molecule from each sample pair such that each sample pair is associated with a corresponding barcode;

(c) sequencing regions of at least one template target nucleic acid molecule in each first sample to obtain a plurality of unmutated sequence reads, the sequencing being performed separately for each first sample, thereby providing a corresponding set of unmutated sequence reads for each first sample;

(d) pooling the second samples to form a pool of second samples;

(e) introducing mutations into template target nucleic acid molecules in the sample pool to produce mutated template target nucleic acid molecules;

(d) sequencing regions of the mutated template target nucleic acid molecules to obtain multiple mutated sequence reads.

The sequencing step may optionally further comprise fragmenting at least one template target nucleic acid molecule in each first sample after insertion of the sample barcode and before sequencing regions of the at least one template target nucleic acid molecule. The sequencing step may optionally further comprise fragmenting at least one template target nucleic acid molecule or mutated template target nucleic acid molecules in the sample pool prior to sequencing regions of the mutated template target nucleic acid molecules.

In the methods of the invention, any number of at least one template target nucleic acid molecule can be sequenced simultaneously. Thus, in one embodiment of the invention, at least one template target nucleic acid molecule comprises a plurality of template target nucleic acid molecules. At least one template target nucleic acid molecule optionally contains at least 10, at least 20, at least 50, at least 100, or at least 250 template target nucleic acid molecules. At least one template target nucleic acid molecule optionally contains 10 to 1000, 20 to 500, or 50 to 100 template target nucleic acid molecules.

Step S110: sample preparation

Since the first sample of the sample pair and the second sample of the sample pair both contain at least one template target nucleic acid molecule, the sample pair can be obtained from the same target organism or taken from the same source sample.

For example, if the user plans to sequence at least one template target nucleic acid molecule in a sample, the user may take a pair of samples from the same source sample. The user may optionally replicate the at least one template target nucleic acid molecule in the original sample before extracting a sample pair therefrom. The user may plan to sequence different nucleic acid molecules of a particular organism, such as E. coli. In this case, the first sample of the sample pair may be an E. coli sample from one source, and the second sample of the sample pair may be an E. coli sample from a second source.

The sample pair may be from any source that contains, or is suspected to contain, at least one template target nucleic acid molecule. The sample pair may comprise a sample of nucleic acid molecules obtained from a human, such as a sample swabbed from the skin of a human patient. In an alternative embodiment, a pair of samples may be obtained from other sources, such as a water source. Such samples can contain billions of template nucleic acid molecules. It may be possible to sequence each of these billions of nucleic acid molecules simultaneously using the methods of the invention, and therefore there is no upper limit on the number of template target nucleic acid molecules that can be used in the methods of the invention.

In one embodiment, multiple pairs of samples may be provided. For example, at least 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 15, 20, 25, 50, 75, 100, 500, 1000, or 5000 sample pairs may be provided. Optionally provide less than 10000, less than 5000, less than 1000, less than 100, less than 75, less than 50, less than 25, less than 20, less than 15, less than 11, less than 10, less than 9, less than 8, less than 7, less than 6, less than 5 or less than 4 samples. Optionally, 2 to 100, 2 to 75, 2 to 50, 2 to 25, 5 to 15, or 7 to 15 sample pairs are provided.

By providing multiple sample pairs, at least one template target nucleic acid molecule in different sample pairs can be labeled with different sample marks (also referred to herein as barcodes). For example, if the user plans to provide 2 sample pairs, then all or substantially all of at least one template target nucleic acid molecule in the first sample pair may be labeled with sample label A, and all or substantially all of at least one molecule template target nucleic acids in the second pair of samples can be labeled with the label B of the sample.

Suitable methods for amplifying at least one template target nucleic acid molecule are known in the art. For example, PCR is commonly used. PCR is described in more detail below under the heading "introduction of mutations in at least one molecule of the template target nucleic acid".

Introduction of mutations into at least one molecule of the template target nucleic acid

The method may include the step of introducing mutations into at least one template target nucleic acid molecule in a second sample from a pair of samples to obtain at least one mutated template target nucleic acid molecule.

Mutations can be substitution type mutations, insertion type mutations or deletion type mutations. For the purposes of the present invention, the term "mutation by type of substitution" should be interpreted as that the nucleotide is replaced by another nucleotide. For example, converting an ATCC sequence to an AGCC sequence introduces one substitution type mutation. For the purposes of the present invention, the term "insert mutation" should be interpreted as meaning that at least one nucleotide is added to the sequence. For example, the conversion of an ATCC sequence to an ATTCC sequence is an example of an insertion-type mutation (with the insertion of an additional T nucleotide). For the purposes of the present invention, the term "deletion type mutation" should be interpreted as that at least one nucleotide is removed from the sequence. For example, the conversion of the ATTCC sequence to ATCC is an example of a deletion-type mutation (with the removal of the T nucleotide). Preferably the mutations are substitution type mutations.

The phrase "introducing mutations into at least one template target nucleic acid molecule" means subjecting at least one template target nucleic acid molecule to the second of a pair of samples of conditions in which at least one template target nucleic acid molecule mutates. This can be achieved using any suitable method. For example, mutations can be introduced by chemical mutagenesis and/or enzymatic mutagenesis.

The step of introducing mutations into at least one template target nucleic acid molecule optionally results in a mutation of 1% to 50%, 3% to 25%, 5% to 20%, or about 8% of the nucleotides of at least one template nucleic acid molecule. target acids. At least one mutated template target nucleic acid molecule optionally contains 1% to 50%, 3% to 25%, 5% to 20%, or about 8% mutations.

The user can determine how many mutations are introduced into at least one mutated template target nucleic acid molecule and/or the level to which, at the stage of introducing mutations into at least one template target nucleic acid molecule, mutations are introduced into at least one template target nucleic acid molecule. target nucleic acid by performing the step of introducing mutations into a nucleic acid molecule with a known sequence, sequencing the resulting nucleic acid molecule and determining the percentage of the total number of nucleotides that have undergone a change compared to the original sequence.

The step of introducing mutations into at least one template target nucleic acid molecule optionally results in the mutation of at least one template target nucleic acid molecule in a substantially random manner. At least one mutated template target nucleic acid molecule optionally contains a substantially random pattern of mutations.

At least one mutated template target nucleic acid molecule contains a substantially random pattern of mutations if it contains mutations along its entire length in substantially similar amounts. For example, a user can determine whether at least one mutated template target nucleic acid molecule contains a substantially random pattern of mutations by introducing mutations into a test nucleic acid molecule with a known sequence to produce a mutated test nucleic acid molecule. The sequence of the mutated test nucleic acid molecule can be compared to the test nucleic acid molecule and the positions of each of the mutations determined. The user can then determine if mutations occur throughout the length of the mutated test nucleic acid in substantially similar numbers as follows:

(i) calculate the distances between each of the mutations;

(ii) calculate the average for the distances;

(iii) form a subsample of distances without substitution for a smaller number, for example from 500 or 1000 distances;

(iv) construct a model set of 500 or 1000 geometrically distributed distances, with a moment mean to compare with that previously obtained at the observed distances; And

(v) run the Kolmogorov-Smirnov test on these two distributions.

At least one mutated template target nucleic acid molecule can be considered to contain a substantially random pattern of mutations at D<0.15, D<0.2, D<0.25 or D<0.3 depending on the length of unmutated reads.

Similarly, in the step of introducing mutations into at least one template target nucleic acid molecule, a mutation is introduced into at least one template target nucleic acid molecule in a substantially random manner, if the resulting at least one mutated template target nucleic acid contains essentially random pattern of mutations. It can be determined whether, in the step of introducing mutations into at least one template target nucleic acid molecule, mutations are introduced into at least one template target nucleic acid molecule in a substantially random manner, it is possible by performing the step of introducing mutations into at least one template target nucleic acid molecule. target nucleic acid on a test nucleic acid with a known sequence to obtain a mutated test nucleic acid molecule. The user can then sequence the mutated test nucleic acid molecule to identify which mutations have been introduced and determine if the mutated nucleic acid molecule contains a substantially random pattern of mutations.

At least one mutated template target nucleic acid molecule optionally contains an unbiased mutation pattern. In the step of introducing mutations, at least one molecule of the template target nucleic acid is optionally mutated in an unbiased manner. At least one mutated template target nucleic acid molecule contains an unbiased mutation pattern if the types of mutations introduced are random. If the introduced mutations are substitution type mutations, then said mutations are random if similar proportions of A (adenosine), T (thymine), C (cytosine) and G (guanine) are introduced. The phrase "similar proportions of A (adenosine), T (thymine), C (cytosine), and G (guanine) are administered" means that the amounts of adenosine, thymine, cytosine, and guanine nucleotides administered are within 20% of each other (e.g., 20 A-nucleotides, 18 T-nucleotides, 24 C-nucleotides and 22 G-nucleotides).

It can be determined whether, in the step of introducing mutations into at least one template target nucleic acid molecule, mutations are introduced into at least one template target nucleic acid molecule in an unbiased manner, by performing the step of introducing mutations into at least one template target nucleic acid molecule targeting a test nucleic acid with a known sequence to produce a mutated test nucleic acid molecule. The user can then sequence the mutated test nucleic acid molecule, identify which mutations were introduced, and determine if the mutated nucleic acid molecule contains an unbiased mutation pattern.

In general, methods for obtaining the sequence of at least one template target nucleic acid molecule can be used even if the step of introducing mutations into at least one template target nucleic acid molecule introduces non-uniformly distributed mutations. Thus, in one embodiment, at least one mutated template target nucleic acid molecule contains non-uniformly distributed mutations. In the step of introducing mutations, at least one mutated template target nucleic acid molecule is optionally introduced with mutations that are unevenly distributed. Mutations are considered "unequally distributed" if the mutations are introduced in a biased manner, i.e. the amounts of adenosine, thymine, cytosine and guanine nucleotides introduced are not within 20% of each other. It can be determined whether at least one mutated template target nucleic acid molecule actually contains non-uniformly distributed mutations or whether mutations are introduced into at least one template target nucleic acid molecule at the stage of introducing mutations that are not uniformly distributed, in the same way as described above in with respect to determining whether mutations are introduced in the step of introducing mutations into at least one template target nucleic acid molecule in an unbiased manner.

Similarly, methods for obtaining the sequence of at least one template target nucleic acid molecule can be used even if mutated sequence reads and/or non-mutated sequence reads contain unevenly distributed sequencing errors. Thus, in one embodiment, mutated sequence reads and/or non-mutated sequence reads contain sequencing errors that are unevenly distributed. Similarly, in one embodiment, the step of sequencing regions of at least one template target nucleic acid molecule and/or sequencing regions of at least one mutated template target nucleic acid molecule introduces sequencing errors that are unevenly distributed.

Whether a particular step of sequencing regions of at least one template target nucleic acid molecule and/or sequencing regions of at least one mutated target nucleic acid molecule introduces sequencing errors that are unevenly distributed is likely to depend on the accuracy of the sequencing tool and, likely to be known to the user. However, the user can ascertain whether the step of sequencing regions of at least one template target nucleic acid molecule and/or sequencing regions of at least one mutated template target nucleic acid molecule introduces sequencing errors that are unevenly distributed by performing the sequencing method on the molecule nucleic acid with a known sequence and comparing the obtained sequence reads with sequence reads obtained using the original nucleic acid molecule with a known sequence. The user can then apply the probability function described in Example 6 and determine values for M and E. target nucleic acid introduces sequencing errors that are unevenly distributed.

The introduction of mutations into at least one template target nucleic acid molecule by chemical mutagenesis can be achieved by exposing at least one template target nucleic acid molecule to a chemical mutagen. Suitable chemical mutagens include mitomycin C (MMC), n-methyl-N-nitrosourea (MNU), nitrous acid (NA), diepoxybutane (DEB), 1,2,7,8,-diepoxyoctane (DEO), ethyl methanesulfonate ( EMS), methyl methanesulfonate (MMS), N-methyl-N'-nitro-N-nitrosoguanidine (MNNG), 4-nitroquinoline-1-oxide (4-NQO), 2-methyloxy-6-chloro-9(3-[ ethyl-2-chloroethyl]-aminopropylamino)-acridine dihydrochloride (ICR-170), 2-aminopurine (2A), bisulfite and hydroxylamine (HA). For example, when nucleic acid molecules are exposed to bisulfite, bisulfite deaminates cytosine to form uracil, effectively introducing a C-T substitution mutation.

As noted above, the step of introducing mutations into at least one molecule of the template target nucleic acid can be carried out by enzymatic mutagenesis. Enzymatic mutagenesis is optionally carried out using DNA polymerase. For example, some DNA polymerases are error prone (low fidelity polymerases) and replication of at least one template target nucleic acid molecule using the error prone DNA polymerase will introduce mutations. Taq polymerase is an example of a low fidelity polymerase, and the step of introducing mutations into at least one molecule of the template target nucleic acid can be carried out by replicating at least one molecule of the template target nucleic acid using Taq polymerase, for example, by PCR .

The DNA polymerase may be a low bias DNA polymerase.

If the step of introducing mutations into at least one template target nucleic acid molecule is carried out using a DNA polymerase, the at least one template target nucleic acid molecule can be incubated with DNA polymerase and suitable primers under conditions suitable to The DNA polymerase has catalyzed the formation of at least one mutated template target nucleic acid molecule.

Suitable primers comprise short nucleic acid molecules that are complementary to regions flanking at least one template target nucleic acid molecule or flanking regions of nucleic acid molecules that are complementary to at least one template target nucleic acid molecule. For example, if at least one template target nucleic acid molecule is part of a chromosome, the primers will be complementary to regions of the chromosome immediately 3' to 3' of at least one template target nucleic acid molecule and immediately 5' to 5 the '-terminus of at least one template target nucleic acid molecule, or the primers will be complementary to regions of the chromosome directly from the 3' to the 3'-terminus of the nucleic acid molecule complementary to at least one molecule of the template target nucleic acid, and immediately from 5 ' to the 5' end of the nucleic acid molecule complementary to at least one template target nucleic acid molecule.

Suitable conditions include a temperature at which the DNA polymerase can replicate at least one molecule of the template target nucleic acid. For example, a temperature of 40°C to 90°C, 50°C to 80°C, 60°C to 70°C, or about 68°C.

The step of introducing mutations into at least one nucleic acid molecule may include multiple cycles of replication. For example, the step of introducing mutations into at least one molecule of the template target nucleic acid preferably comprises:

i) replicating at least one template target nucleic acid molecule to produce at least one nucleic acid molecule that is complementary to at least one template target nucleic acid molecule; And

ii) a cycle of replication of at least one template target nucleic acid molecule to produce replicas of at least one template target nucleic acid molecule.

The step of introducing mutations into at least one molecule of the template target nucleic acid optionally comprises at least 2, at least 4, at least 6, at least 8, at least 10, less than 10, less than 8, about 6, 2 to 8 or 1 to 7 cycles of replication of at least one template target nucleic acid molecule. The user may choose to use a low number of replication cycles to reduce the possibility of introducing amplification bias.

The step of introducing mutations into at least one molecule of the template target nucleic acid optionally comprises at least 2, at least 4, at least 6, at least 8, at least 10, less than 10, less than 8, about 6, 2 to 8 or 1 to 7 replication cycles at 60°C to 80°C.

The step of introducing mutations into at least one molecule of the template target nucleic acid is optionally carried out using polymerase chain reaction (PCR). PCR is a process that includes many cycles of the following steps in the replication of a nucleic acid molecule:

a) melting;

b) annealing; And

c) completion and extension.

The nucleic acid molecule (eg, at least one template target nucleic acid molecule) is mixed with suitable primers and a polymerase. In the melting step, the nucleic acid molecule is heated to a temperature above 90° C. so that the double-stranded nucleic acid molecule is denatured (split into two strands). In the annealing step, the nucleic acid molecule is cooled to a temperature below 75°C, for example from 55°C to 70°C, about 55°C or about 68°C, so that the primers are attached to the nucleic acid molecule. At the stages of completion and elongation, the nucleic acid molecule is heated to a temperature of more than 60°C, so that the DNA polymerase catalyzes the completion of the primer - the addition of nucleotides complementary to the template chain.

The step of introducing mutations into at least one template target nucleic acid molecule optionally includes replication of at least one template target nucleic acid molecule using Taq polymerase under reaction conditions that make it error prone. For example, the step of introducing mutations into at least one molecule of the template target nucleic acid may include PCR using Taq polymerase in the presence of Mn ²⁺ , Mg ²⁺ , or unequal concentrations of dNTPs (for example, excess cytosine, guanine, adenine, or thymine ).

Step S120: sequencing

Obtaining Data Containing Unmutated Sequence Reads and Mutated Sequence Reads

The methods of the invention may include the step of receiving mutated sequence reads and optionally receiving unmutated sequence reads. Unmutated sequence reads and mutated sequence reads can be obtained from any source.

Unmutated sequence reads are optionally obtained by sequencing regions of at least one template target nucleic acid molecule in the first sample of a pair of samples. Mutated sequence reads are optionally obtained by introducing mutations into at least one template target nucleic acid molecule in the second sample of a pair of samples to obtain at least one mutated template target nucleic acid molecule and sequencing regions of at least one mutated template nucleic acid molecule -targets.

Unmutated sequence reads optionally comprise the region sequences of at least one template target nucleic acid molecule in the first sample of a sample pair, and the mutated sequence reads comprise the region sequences of at least one mutated template target nucleic acid of the second sample of the sample pair, and the sample pair was taken from the same original sample or was obtained from the same organism.

Sequencing regions of at least one template target nucleic acid molecule or at least one mutated template target nucleic acid molecule

The method for determining the sequence of at least one template target nucleic acid molecule may include the step of sequencing regions of at least one template target nucleic acid molecule in a first sample from a pair of samples to obtain unmutated sequence reads and/or the step of sequencing at least one mutated molecule template target nucleic acid to obtain mutated sequence reads.

The sequencing steps can be performed using any sequencing method. Examples of possible sequencing methods include Maxam-Gilbert sequencing, Sanger sequencing, sequencing involving bridged amplification (e.g. bridged PCR) or any high throughput sequencing (HTS) technique as described in Maxam AM, Gilbert W (February 1977) , "A new method for sequencing DNA", Proc. Natl. Acad. sci. U.S.A. 74 (2): 560-4, Sanger F, Coulson AR (May 1975), "A rapid method for determining sequences in DNA by primed synthesis with DNA polymerase)}, J. Mol. Biol. 94(3):441-8; and Bentley DR, Balasubramanian S, et al. (2008), ((Accurate whole human genome sequencing using reversible terminator chemistry", Nature, 456 (7218): 53-59.

In a typical embodiment, at least one or preferably both sequencing steps include bridged amplification. The bridge amplification step is optionally performed using a build-in time of greater than 5, greater than 10, greater than 15, or greater than 20 seconds. An example of the use of bridging amplification are sequencers for analyzing the Illumina genome. Preferably, paired-end sequencing is used.

The steps of (i) sequencing regions of at least one template target nucleic acid molecule in a first sample of a pair of samples to obtain unmutated sequence reads and (ii) sequencing at least one mutated template target nucleic acid molecule to obtain mutated sequence reads are optionally performed using the same sequencing method.

The steps of (i) sequencing regions of at least one template target nucleic acid molecule in a first sample of a pair of samples to obtain unmutated sequence reads and (ii) sequencing at least one mutated template target nucleic acid molecule to obtain mutated sequence reads are optionally performed using different sequencing methods.

The steps of (ii) sequencing regions of at least one template target nucleic acid molecule in a first sample of a pair of samples to obtain unmutated sequence reads and (ii) sequencing at least one mutated template target nucleic acid molecule to obtain mutated sequence reads may optionally be be performed using more than one sequencing method. For example, a fraction of at least one template target nucleic acid molecule in a first sample of a sample pair can be sequenced using a first sequencing method, and a fraction of at least one template target nucleic acid molecule in a first sample of a sample pair can be sequenced with using the second sequencing method. Similarly, a fraction of at least one mutated template target nucleic acid molecule can be sequenced using a first sequencing method, and a fraction of at least one mutated template target nucleic acid molecule can be sequenced using a second sequencing method.

The steps of (i) sequencing regions of at least one template target nucleic acid molecule in a first sample of a pair of samples to obtain unmutated sequence reads and (ii) sequencing at least one mutated template target nucleic acid molecule to obtain mutated sequence reads are optionally performed at different times. In an alternative embodiment, steps (i) and (ii) can be carried out relatively simultaneously, for example within 1 year of each other. The first sample of the sample pair and the second sample of the sample pair need not be taken simultaneously with each other. If two samples are obtained from the same organism, they can be provided at essentially different times, even with a multi-year interval, and therefore the two sequencing steps can also be separated by a period of several years. In addition, even if the first sample of a sample pair and the second sample of a sample pair were obtained from the same original sample, biological samples can be stored for some time, and therefore it is not necessary to perform the sequencing steps at the same time.

Mutated sequence reads and/or non-mutated sequence reads can be single-ended or double-ended sequence reads.

Mutated sequence reads and/or unmutated sequence reads are optionally greater than 50 nt, greater than 100 nt, greater than 500 nt, less than 200,000 nt, less than 15,000 nt, less than 1000 nt, 50 to 200,000 nt, 50 to 15,000 nt, or 50 up to 1000 nt.

The sequencing steps are optionally performed using a sequencing depth of 0.1 to 500 reads, 0.2 to 300 reads, or 0.5 to 150 reads per nucleotide per at least one template target nucleic acid molecule. The deeper the sequencing, the higher the accuracy of the determined/generated sequence will be, but the assembly may be more difficult.

Choice of options

Preferably, the parameters used in the method 200 are selected as follows.

In a preferred embodiment, the weight w(ψ) of each seed pattern is in the range of 5 to 50, preferably 10 to 30, further preferably 13 to 23. This ensures that each seed pattern is large enough to ensure that each k -measure masked by each seed pattern ψ is unique with a high probability. For example, for bacterial genomes with a typical length of 5 million nucleotides, the weight w(ψ) of each seed pattern ψ is preferably in the range of 13-19, given that 4 ¹³ >5 million. For genomes similar in size to humans, with a typical length of about 3 billion nucleotides, the weight w(ψ) of each seed pattern is preferably in the range of 19-23, given that 4 ¹⁹ >3×10 ⁹ .

In a preferred embodiment, the size k of each k-mer used in step S230 to determine the positions of one or more mutations in each mutated sequence read is greater than the weight w(ψ) of each seed pattern. The size k of each k-measure can be less than 5 times, less than 4 times, less than 3 times, or less than 2 times the weight w(ψ) of each seed pattern ψ. The size k of each k-mer used in step S230 to determine the positions of one or more mutations in each mutated sequence read may range from 10 to 250, preferably from 13 to 100, further preferably from 15 to 50, most preferably from 20 to 40. This ensures that the size k is small enough to provide a low probability that any k-mer will include an insertion or deletion sequencing error, which is a disadvantage in the context of method 200.

Below is an example of a family of seed patterns containing seed patterns with weight w(ψ)=16 and k=27:

ψ ₁ ={0, 1, 2, 3, 5, 6, 9, 12, 13, 14, 16, 18, 20, 21, 22, 23},

ψ ₂ \u003d {0, 1, 2, 4, 5, 9, 10, 11, 13, 18, 19, 21, 23, 24, 25, 26},

ψ ₃ ={0, 1, 2, 3, 4, 5, 7, 8, 9, 10, 13, 15, 16, 18, 19, 20},

ψ ₄ ={0, 1, 2, 4, 6, 7, 12, 14, 16, 17, 20, 21, 23, 24, 25, 26},

In one embodiment, the k-measures used in step S220 of applying the overall minimizing function, i.e. the one or more minimizers determined for each mutated sequence read have a size k different from the k-mers used in step S230 to determine the positions of one or more mutations in each mutated sequence read. The size k of each minimizer can range from 5 to 50, preferably from 10 to 30, further preferably from 13 to 23. The size k of each minimizer can be chosen based on the same considerations as choosing the weight w(ψ) of the seed patterns. The k-size of each minimizer can range from 13 to 19 for bacteria and from 19 to 23 for genomes similar in size to humans.

Implementation of method 200

Method 200 may have different implementations. The preferred approach is to first calculate the set U _M in an initial pass through some or all of the mutated P sequence reads and unmutated R sequence reads, then calculate W _M in a second pass through the mutated P sequence reads and unmutated R sequence reads. With W _M , in a third pass through the mutated P sequence reads, the positions of the minimizers can be computed along with the positions of one or more mutations, and these positions can be stored in minimizer groups either in RAM or in a persistent storage device (e.g., on a disk). Multiple minimizer groups may optionally be stored in the same file, either sorted or unsorted. Each minimizer group (or each file) can then be read sequentially or in parallel, processing the minimizer groups to calculate edge weights. Because each mutated sequence read may occur in multiple minimizer groups, it is possible that a pair of mutated sequence reads may have multiple computed weight estimates. In this case, it is necessary to use some indicator to select the preferred weight, usually the maximum. Finally, if the sequencing chemistry produced paired-end reads, and each read in a pair of paired-end reads has common minimizers, then the scores for the two ends can be summed to give one score for the pair of paired-end reads.

Experimental data

Method 200 was used to process multiple real SAM datasets, with each SAM dataset containing unmutated sequence reads and mutated sequence reads.

A SAM dataset of Arobacter butzleri JV22 was processed. This organism has a 2.3 million bp genome that exists as a single circular chromosome. The C++ implementation of method 200 was performed on an Amazon AWS service instance. The SAM dataset consists of 956,133 pairs of reference (unmutated) reads and 2,154,909 pairs of mutated reads, derived from approximately 8,000 mutated long templates. 2,087,506 mutated reads (96.9%) originate from the interior of the mutated long templates, while 67,403 (3.1%) originate from the ends of the long templates and contain sample barcodes. Each individual read is 150 nt or less in length. Pairs of reads were pre-adapter cropped and quality-enhanced cropped. Method 200 required 12 minutes of CPU time and 1.2 GB of RAM to process the dataset, yielding 30,033,939 potential read relationships. These links were then subjected to graph clustering using Markov clustering (mcl) and the resulting 6779 read clusters were assembled de novo using MEGAHIT (see Dinghua Li, Chi-Man Liu, Ruibang Luo, Kunihiko Sadakane, and Tak-Wah Lam. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph Bioinformatics, Oxford, England, 31(10): 1674{1676, May 2015) to obtain reconstructions of long mutated templates. Finally, long mutated templates were used along with unmutated reads in a hybrid genome assembly computed by the Unicycler software. The resulting assembly is shown in Fig. 4B compared to an assembly obtained from short reads only (shown in FIG. 4A). On FIG. 4A shows assembly from short reads of the 2.3 Mb Arcobacter butzlerii genome. using a Shovill assembly line prior to performing method 200. This gave an assembly of 78 scaffolds, with the largest scaffold covering 342 kb and the N50 scaffold at approximately 127 kb. On FIG. 4B shows the assembly of the 2.3 Mbp genome of Arcobacter butzlerii. using method 200. The ring chromosome is largely resolved into a single contig, leaving a copy number of a small 200 nt region unresolved.

The scalability and resolution power of the approach implemented in method 200 was measured using simulated data. Sequence 50 kb. from the CFTR gene was used to model increasing coverage values with mutated long templates and corresponding mutated short reads from these templates. Simulations were performed using newly developed scripts that first generate long mutated templates, then activate the well-known Illumina read simulator called artsim to simulate sequencing from mutated templates over short reads. In addition to the mutated data, a 30-fold coverage of the unmutated sequence was simulated in artsim. We simulated coverage with long mutated templates in the range from 10 ¹ to 10 ⁶ with increments of an order of magnitude. The proportion of mutations was fixed at 6%. For each long template, 10x coverage by short reads was simulated. The results for the simulated data were evaluated by measuring the proportion of false positive associations in method 200.

On FIG. Figure 5 shows the effect of coverage depth with short reads of a long template. The amount of data with short reads per long template is shown on the x-axis, and the y-axis shows various performance measures for the results of method 200. It can be seen that with low template coverage with short reads, e.g. <4x, poor and incomplete reconstructions of the original long templates are obtained. However, when the coverage of the mutated template is in the 5-10x range, good reconstructions can be obtained.

• links num: number of links between mutated reads according to method 200. links fp: number of flagged false positive links.

• links fp rate: proportion of false positive links from all received links.

• mcl num: number of clusters created by Markov clustering for the graph, as reported by mmdreaming.

• idba scaf num: number of scaffold sequences reconstructed by clustering mutated short reads.

• idba scaf bp: the sum of the lengths of all scafs assembled.

Claims (54)

1. A computer-implemented method for determining the sequence of at least a portion of at least one template target nucleic acid by determining whether two mutated sequence reads originate from the same mutated sequence, comprising: receiving a plurality of mutated sequence reads, each mutated sequence read corresponding to a subsequence of a mutated sequence, wherein the mutated sequence contains mutations compared to a non-mutated sequence; applying a common minimizer function for each mutated sequence read, thereby determining one or more corresponding minimizers for each mutated sequence read; determining the positions of one or more corresponding minimizers in each mutated sequence read; determining the positions of one or more mutations in each mutated sequence reading; And for at least two mutated sequence reads with a common minimizer, counting the number of match and/or mismatch mutations when the respective minimizers are aligned to determine a score correlated with the probability that said at least two mutated sequence reads occur from the same mutation-containing sequence, assembling said at least two mutated sequence reads based on said index; And determining the sequence of at least a portion of at least one template target nucleic acid based on said assembly. 2. The method of claim 1, further comprising receiving a plurality of unmutated sequence reads, each unmutated sequence read corresponding to a subsequence of the unmutated sequence. 3. The method of claim 1 or 2, wherein the step of applying a common minimizing function to each mutated sequence reading comprises identifying i) one or more k-mers in the corresponding mutated sequence reading that(s) are(s) listed first( -i) in an ordered list of possible k-mers, or ii) one or more k-mers that occur in a predefined set of possible k-mers, wherein the one or more minimizers defined for the corresponding mutated sequence reading are the identified one or more than k-mers. 4 . The method of claim 3, wherein i) in the ordered list of possible k-mers, the k-mers are ordered based on the probability that the k-mers occur in a mutation-containing sequence and do not occur in a mutation-free sequence, or ii) a predetermined the set of candidate k-mers contains k-mers that are relatively unlikely to occur in the mutated sequence but not in the non-mutated sequence, and optionally the predefined set of possible k-mers does not contain k-mers that are relatively unlikely to occur in the mutated sequence . 5. The method according to paragraphs. 2 and 3 or 4, wherein the ordered list of possible k-mers or the predefined set of possible k-mers consists of k-mers that occur more frequently in the set of mutated sequence reads than in the set of unmutated sequence reads, optionally k-mers, which occur more frequently in the set of mutated sequence reads than in the set of unmutated sequence reads are relatively likely to occur in the sequence containing the mutation. 6. The method according to p. 2 and any of paragraphs. 3-5, in which the predetermined set of possible k-mers consists of k-mers that occur n or more times in the set of mutated reads of the sequence, and occur less than n times in the set of unmutated reads of the sequence, where n is an integer, which is greater than or equal to 1, optionally k-mers that occur n or more times in the set of mutated sequence reads and occur less than n times in the set of unmutated sequence reads are relatively likely to occur in the sequence containing the mutation. 7. The method of claim 6, wherein the predetermined set of possible k-mers consists of k-mers that do not occur in the set of unmutated sequence reads. 8. The method according to claim 6 or 7, in which n is 2. 9. The method according to p. 2 and any of paragraphs. 3-8, further comprising generating an ordered list of possible k-mers or a predefined set of possible k-mers based on a comparison of the k-mers in the set of mutated sequence reads and the k-mers in the set of unmutated sequence reads. 10. A method according to any one of the preceding claims, wherein each minimizer is a k-mer longer than 5, preferably greater than 10. 11. The method of any one of the preceding claims, further comprising distributing the mutated sequence reads into one or more minimizer groups such that each minimizer group contains mutated sequence reads having a common minimizer and does not contain mutated sequence reads that do not have a common minimizer, and however, the step of counting the number of mutations with a matching position and/or with a mismatching position can be performed only on mutated sequence reads that are in the same minimizer group. 12. The method according to any one of paragraphs. 1-11, wherein the step of determining the positions of one or more mutations in each mutated sequence read includes: obtaining a set of unmutated primed masked k-mers by applying each of the one or more seed patterns to the k-mers in a plurality of unmutated sequence reads; for each mutated sequence read, apply one or more seed patterns to the k-mers in the corresponding mutated sequence read to obtain a plurality of mutated primed-masked k-mers, and determine the positions of one or more mutations by identifying one or more positions in the mutated sequence read that are masked by all seed patterns corresponding to the mutated primed masked k-mers from the set of mutated seeded masked k-mers that occur in the set of unmutated seeded masked k-mers. 13. The method of claim 12, wherein the one or more seed patterns are selected such that the probability of obtaining identical seed masked k-mers when at least one of the one or more seed patterns is applied to any k-mer of a plurality of mutated sequence reads and the corresponding k-mer of the set of unmutated sequence reads is greater than 90%, preferably greater than 99%. 14. The method of claim 12 or 13, wherein the mutated sequence contains transition-type mutations compared to a sequence that does not contain mutations; And wherein the one or more seed patterns may be one or more transition seed patterns. 15. The method according to any one of paragraphs. 12-14, wherein each of the plurality of mutated sequence reads corresponds to a subsequence of a mutated sequence associated with one of the plurality of samples, and each of the plurality of unmutated sequence reads corresponds to a subsequence of a mutant-free sequence associated with one of the plurality of samples, each containing mutations the sequence contains mutations compared to the corresponding non-mutated sequence; wherein obtaining a set of unmutated primed masked k-mers comprises obtaining a corresponding set of unmutated primed masked k-mers for each sample; the method further comprising generating a set of unmutated sample bit vectors, wherein each unmutated sample bit vector determines, for a corresponding k-mer in the set of unmutated primed masked k-mers, in which of the set of samples the corresponding k-mer occurs; And wherein, for each mutated sequence read and for each set and/or each combination of sets of unmutated primed masked k-mers, determining the positions of one or more mutations includes identifying one or more positions in the mutated sequence reading that are masked by all seed patterns corresponding to the mutated primed masked k-mers from a plurality of mutated primed-masked k-mers that occur in the corresponding set or combination of unmutated primed-masked k-mer sets, and linking the identified one or more positions to one or more patterns associated with the corresponding set or combination of unmutated primed masked sets k-mers. 16. The method according to any one of paragraphs. 2-15, wherein the step of determining the positions of one or more mutations in each mutated sequence read includes: for one or more of the mutated sequence reads, aligning the corresponding mutated sequence read with the reference assembly; And determining the positions of one or more mutations in the corresponding mutated sequence reading by identifying, in the corresponding mutated sequence reading, position differences between the corresponding mutated sequence reading and the reference assembly. 17. The method according to p. 16 and depending on any of paragraphs. 12-15, wherein for each mutated sequence read, the step of determining the positions of one or more mutations in each mutated sequence read includes: if the position in the corresponding mutated sequence reading is aligned with the reference assembly, determining the position in the corresponding mutated sequence reading as the position of the mutation in the corresponding mutated sequence reading if the position in the corresponding mutated sequence reading is the position at which the corresponding mutated sequence reading differs from the reference assembly; And if the position in the corresponding mutated sequence reading is not aligned with the reference assembly, determining the position in the corresponding mutated sequence reading as the position of the mutation in the corresponding mutated sequence reading if the position in the corresponding mutated sequence reading is a position that is masked by all seed patterns that match the masked primers mutated k-mers from the set of seeded mutated k-mers that occur in the set of seeded unmutated k-mers. 18. The method of any one of the preceding claims, comprising determining a score correlated with the probability that at least two mutated sequence reads are from the same mutated sequence based on the number of matched and/or mismatched mutations. 19. The method of claim 18, wherein the metric correlated with the probability that at least two mutated sequence reads are from the same mutated sequence is one of: i) a probability density that at least two mutated sequence reads are from the same mutated sequence, and ii) a scoring function that is correlated with a probability density that at least two mutated sequence reads are from the same mutated sequence. 20. The method according to any of the preceding claims, further comprising creating an undirected weighted graph from a plurality of mutated sequence reads, moreover, the undirected weighted graph contains nodes corresponding to a plurality of mutated reads of the sequence, and while the edges between the nodes are associated with the corresponding weight values of the edges, and the weight of each edge is determined based on the number of mutations with the same position and/or with the mismatch position determined for the two mutated reads sequences corresponding to two nodes associated with the corresponding edge. 21. The method of claim 20, wherein the edge weights correspond to a score correlated with the probability that at least two mutated sequence reads corresponding to two nodes associated with the respective edge come from the same mutated sequence. 22. The method of claim 20 or 21, further comprising performing a graph clustering operation on the undirected weighted graph, thereby obtaining clusters of mutated sequence reads that are expected to come from the same mutated sequence. 23. The method of claim 22, wherein the graph clustering includes Markov or Infomap clustering. 24. The method of claim 22 or 23, further comprising reconstructing at least a portion of the mutated sequence by assembling the mutated sequence reads into clusters. 25. The method according to paragraphs. 2 and 24 further comprising reconstructing at least a portion of the mutant-free sequence by deriving at least a portion of the probable mutant-free sequence from the reconstructed portion of the mutated sequence, optionally using multiple unmutated sequence reads. 26. The method of obtaining at least part of the sequence of the target nucleic acid template molecule, including the method according to any one of paragraphs. 20–25. 27. A method for determining at least a portion of the sequence of at least one template target nucleic acid molecule, comprising sequencing regions of at least one template target nucleic acid molecule containing mutations to obtain a plurality of mutated sequence reads, performing the method of any preceding claim on the resulting set of mutated sequence reads. 28. The method of claim 27, wherein the sequencing step comprises (a) providing a pair of samples, each sample containing at least one template target nucleic acid molecule; (b) sequencing regions of at least one template target nucleic acid molecule in a first sample of a pair of samples to obtain multiple unmutated sequence reads; (c) introducing mutations into at least one template target nucleic acid molecule in a second sample from a pair of samples to obtain at least one mutated template target nucleic acid molecule; (d) sequencing regions of at least one mutated template target nucleic acid molecule to obtain a plurality of mutated sequence reads. 29. The method of claim 28, wherein the step of introducing mutations comprises introducing transition-type mutations into at least one template target nucleic acid molecule in a second of a pair of samples.