Methods and systems for sequence alignment with reduced sequence representation
The methods and systems for aligning sequence reads using reduced sequence representations address the challenges of high error rates in long-read sequencing by encoding with a reduced alphabet and generating syncstrobes, enhancing alignment accuracy and efficiency, especially in error-prone sequencing technologies like SAM.
Patent Information
- Application Number
- PCT/US2025/030029
- Authority / Receiving Office
- WO · WO
- Patent Type
- Applications
- Current Assignee / Owner
- Priority Date
- 2024-06-28
- Filing Date
- 2025-05-19
- Publication Date
- 2026-01-02
AI Technical Summary
Long-read sequencing technologies face high error rates and computational challenges in sequence alignment due to introduced mutations and sequencing errors, complicating the alignment and analysis of nucleic acid sequences.
Methods and systems for aligning sequence reads using reduced sequence representations, including encoding with a reduced alphabet, homopolymer quantization, and generating syncstrobes for efficient alignment, which can tolerate high mutation rates and indel differences without requiring a reference sequence.
Improves accuracy and recall in identifying mutated sites, enhances alignment efficiency, and supports high sensitivity and specificity in sequence matching, particularly in error-prone sequencing methods like Sequencing Aided by Mutagenesis (SAM), leading to better quality of sequence reads and downstream analysis.
Smart Images

Figure US2025030029_02012026_PF_FP_ABST
Abstract
Description
PATENTMETHODS AND SYSTEMS FOR SEQUENCE ALIGNMENT WITH REDUCED SEQUENCE REPRESENTATIONCROSS-REFERENCE TO RELATED APPLICATIONS
[0001] This application claims priority to U.S. Prov. No. 63 / 665707 filed June 28, 2025 entitled “METHODS AND SYSTEMS FOR SEQUENCE ALIGNMENT WITH REDUCED SEQUENCE REPRESENTATION” which is incorporated by reference herein in its entirety.REFERENCE TO SEQUENCE LISTING
[0002] The present application is being filed along with a Sequence Listing in electronic format. The Sequence Listing is provided as a file entitled ILLINC841WOSEQ.XML, created on April 14, 2025, which is 11,086 bytes in size. The information in the electronic format of the Sequence Listing is incorporated herein by reference in its entirety.BACKGROUNDField
[0003] The present disclosure generally relates to the field of sequence alignment. In particular, the present disclosure relates to the field of aligning nucleic acid sequence reads by encoding sequences with a reduced sequence representation.Description
[0004] Sequencing Aided by Mutagenesis (SAM) is a technique used to aid sequence alignment by intentionally introducing mutations into nucleic acid molecules during library preparation, and using these mutations to align sequence reads. Long read sequencing technologies may have high error rates which complicate sequence alignment and analysis. Successful implementation of these technologies depends heavily on computational methods for sequence data analysis.SUMMARY
[0005] Disclosed herein are methods for aligning sequence reads based on reduced sequence representations. In some embodiments, the method includes receiving sequence reads generated from nucleic acids; encoding sequence reads with a reduced sequence representation, thereby generating transformed sequence reads; generating seeds from the transformed sequence reads; and matching the seeds, thereby aligning the transformed sequence reads.
[0006] In some embodiments, the method comprises encoding the sequence reads using a reduced alphabet. In some embodiments, the method comprises encoding the sequence reads in RY format.
[0007] In some embodiments, the method comprises reducing a length of homopolymers within the transformed sequence reads. In some embodiments, the method comprises quantizing the homopolymers by rounding the length of homopolymers within the sequence reads to a predetermined value of a set of two or more predetermined values.
[0008] In some embodiments, generating seeds from the transformed sequence reads comprises selecting strand-symmetric pairs of minimizers from the transformed sequence reads. In some embodiments, the method comprises generating a sequence index based on the transformed sequence reads. In some embodiments, generating a sequence index comprises indexing strobemers by selecting an outer pair of k-mers and selecting one or more additional k-mers between the outer pair of k-mers. In some embodiments, the outer pair of k- mers, or the one or more additional k-mers, are selected strand-symmetrically. In some embodiments, selecting k-mers comprises using a minimizer function to select k-mers. In some embodiments, the minimizer function selects k-mers based on a parameterized syncmer scheme. In some embodiments, the minimizer function evaluates k-mers based on whether s- mer minimizers occur in one or more pre-determined positions specified in the parameterized syncmer scheme.
[0009] In some embodiments, generating a sequence index comprises generating a set of syncstrobes by successively removing one or more individual k-mer components from one or more strobemers.
[0010] In some embodiments, the method comprises matching seeds from the sequence index, and generating one or more alignments based on seed matches.
[0011] In some embodiments, the method comprises using k-mer components of a strobemer seed match as anchors in a gapped sequence alignment. In some embodiments, the method comprises scoring an alignment based on whether a distance between two or more k- mer components of a strobemer seed match is consistent in a matched strobemer.
[0012] In some embodiments, the method comprises identifying repetitive strobemers and storing a single representative strobemer to be used as a seed in seed matching. In some embodiments, the repetitive strobemers are marked as repetitive.
[0013] In some embodiments, the method comprises aligning short sequence reads with a length of 50-500 bp to each other. In some embodiments, the method comprises aligning short sequence reads with a length of 50-500 bp to long sequence reads with a length of 50- 250,000 bp. In some embodiments, aligning the transformed sequence reads comprises aligning without a reference sequence.
[0014] In some embodiments, the sequence reads comprise mutations or sequence errors. In some embodiments, 4% to 12% of nucleotides comprise a mutation or a sequence error. In some embodiments, the mutation comprises a transition mutation.
[0015] In some embodiments, the nucleic acids comprise mutated nucleic acids and unmutated nucleic acids, and wherein the sequence reads comprise mutated sequence reads and unmutated sequence reads. In some embodiments, the method comprises identifying mutated positions in mutated sequence reads by comparing alignments between a mutated sequence read set and an unmutated read sequence set. In some embodiments, the method comprises selecting a set of alignments of unmutated reads for each mutated sequence read based on an alignment score. In some embodiments, the method comprises selecting a set of unmutated sequence reads which are consistent with each other, which have a sequence identity above a threshold to a mutated sequence read covering a corresponding sequence, or which together cover at least a predetermined percentage of the mutated sequence read.
[0016] In some embodiments, the method comprises marking mutated positions at sites which differ between the mutated sequence read set and the unmutated sequence read set.
[0017] In some embodiments, the method comprises: generating a sequence index from transformed sequence reads; identifying seed matches from index entries; and generating alignments based on the seed matches. In some embodiments, generating a sequence index from transformed sequence reads comprises: selecting an outer pair of k-mers and selectingone or more additional k-mers between the outer pair of k-mers using a syncmer-preferred minimizer function; and removing one or more k-mer components from one or more strobemers, thereby generating sets of syncstrobes.
[0018] In some embodiments, the method comprises: generating a sequence index comprising seeds generated from sequence reads encoded with reduced sequence representation; matching seeds from the sequence index and generating alignments; scoring alignments; and combining separate alignments of short sequence read pairs into a paired-end alignment.
[0019] In some embodiments, the method comprises: generating a sequence index of syncstrobes from transformed sequence reads, wherein generating a sequence index comprises: selecting an outer pair of k-mers and selecting one or more additional k-mers between the outer pair of k-mers using a syncmer-preferred minimizer function; removing one or more k-mer components from one or more strobemers, thereby generating sets of syncstrobes; ordering sets of syncstrobes; sorting and scanning indexes of the sequence index, wherein the indexes correspond to mutated sequence reads and unmutated sequence reads, thereby identifying one or more matches; generating one or more ungapped, anchored, or gapped alignments based on the one or more matches; scoring one or more ungapped, anchored, or gapped alignments; and combining separate alignments of short sequence read pairs into a paired-end alignment.
[0020] Further disclosed herein are electronic systems for aligning sequence reads based on reduced sequence representations. In some embodiments, the electronic system comprises a processor configured to perform a method comprising: receiving sequence reads generated from nucleic acids; encoding sequence reads with a reduced sequence representation, thereby generating transformed sequence reads; generating seeds from the transformed sequence reads; and matching the seeds, thereby aligning the transformed sequence reads.
[0021] Further disclosed herein are non-transitory computer-readable media. In some embodiments, the non-transitory computer-readable medium comprises a plurality of instructions, which when executed by at least one processor, cause the at least one processor to: receive sequence reads generated from nucleic acids; encode sequence reads with a reduced sequence representation, thereby generating transformed sequence reads; generate seeds fromthe transformed sequence reads; and match the seeds, thereby aligning the transformed sequence reads.BRIEF DESCRIPTION OF THE DRAWINGS
[0022] Features of examples of the present disclosure will become apparent by reference to the following detailed description and drawings, in which like reference numerals correspond to similar, though perhaps not identical, components. For the sake of brevity, reference numerals or features having a previously described function may or may not be described in connection with other drawings in which they appear. In addition to the features described herein, additional features and variations will be readily apparent from the following descriptions of the drawings and exemplary embodiments. It is to be understood that these drawings depict typical embodiments, and are not intended to be limiting in scope.
[0023] FIG. 1 is a flow diagram that schematically illustrates an exemplary method for aligning sequence reads based on reduced sequence representations.
[0024] FIG. 2A is a flow diagram that schematically illustrates a process aligning transformed sequence reads that may take place within the method of FIG. 1.
[0025] FIG. 2B is a flow diagram that schematically illustrates a process of generating a sequence index from transformed sequence reads that may take place within the process of FIG. 2A.
[0026] FIG. 3 A is a block diagram of an exemplary sequencing system that may be used to perform the disclosed methods.
[0027] FIG. 3B is a block diagram of an exemplary computing device that may be used in connection with the exemplary sequencing system of FIG. 3 A.
[0028] FIG. 4A illustrates syncmer pairs and component selection for a syncstrobe.FIG. 4A includes a 30-nucleotide sequence (SEQ ID NO: 01).
[0029] FIG. 4B illustrates strobe sets created by packing the syncstrobes in a set by dropping a component. The upper portion of FIG. 4B includes the 30-nucleotide sequence (SEQ ID NO: 01) and nucleotide sequences of strobe set 1 including (from left to right): SEQ ID NO: 02, SEQ ID NO: 03, and SEQ ID NO: 04. The lower portion of FIG. 4B includes the 30-nucleotide sequence (SEQ ID NO:01) and nucleotide sequences of strobe set 2 including (from left to right): SEQ ID NO:05, SEQ ID NO:06, and SEQ ID NO:07.
[0030] FIG. 4C illustrates syncmer selection with transition mutation and syncstrobe selection with RY-encoding. FIG. 4C includes an unmutated nucleotide sequence (SEQ ID NO:01), a mutated nucleotide sequence (SEQ ID NO:08), and an embodiment of a RY encoded nucleotide sequence (SEQ ID NO: 09).
[0031] FIG. 5 illustrates an exemplary logical layout of entries in the RYE aligner sequence index.
[0032] FIG. 6 illustrates an exemplary logical layout of syncstrobes in the syncstrobe aligner sequence index.
[0033] FIG. 7 is a conceptual illustration of alignment between two reads with a seed-match.
[0034] FIG. 8 illustrates plots of the number and the cumulative number of unique k-mers against their number of occurrences in a genome, and the distribution of unique k-mers in the HG002 T2T assembled genome.
[0035] FIG. 9 is an illustration showing contained-in and paired-with (best paired) reads.
[0036] FIG. 10 is a pair of bar graphs illustrating a comparison of recall of paired- with (best paired) and contained-in reads using matches of various types of seeds. Panel A shows short mutated reads matched with best reference reads and Panel B shows long mutated reads matched with contained reference reads.
[0037] FIG. 11 is a set of bar graphs illustrating a comparison of the accuracy and recall of morphomutation markings in short mutated sequence reads using different methods. Panel A shows data for FNTR, a dataset based on large tandem repeats (>5kbp in length). Panel B shows data for MUTMARK, a dataset based on any length tandem repeats and including 5kbp flanks either side. Panel C shows data for WG T2T, a simulated data set based on the whole genome of the T2T HG002 reference genome.
[0038] FIG. 12 is a bar chart illustrating a comparison of the accuracy and recall based on false negatives (FN) and false positives (FP), of small variants called from sequence reads using different methods for mutation markings.
[0039] FIG. 13 is a pair of bar graphs illustrating a comparison of the accuracy and recall of morphomutation markings in short mutated sequence reads with different mutation rate. Panel A shows data for CHR012122, a simulated dataset based on chromosomes 1, 21,and 22 with a 5% mutation rate. Panel B shows data for CHR012122 MRZ, a simulated dataset based on chromosomes 1, 21, and 22 with a 10% mutation rate.DETAILED DESCRIPTION
[0040] The foregoing and other aspects of the present disclosure will now be described in more detail with respect to the description and methodologies provided herein. This description is not intended to be a detailed catalogue of all the ways in which the embodiments of the present disclosure may be implemented, or of all the features that may be added to the present disclosure. For example, features illustrated with respect to one embodiment may be incorporated into other embodiments, and features illustrated with respect to a particular embodiment may be deleted from that embodiment. In addition, numerous variations and additions to the various embodiments suggested herein, which do not depart from the instant disclosure, will be apparent to those skilled in the art in light of the instant detailed description, figures and claims. Hence, the following specification is intended to illustrate some particular embodiments, and not to exhaustively specify all permutations, combinations and variations thereof.
[0041] All patents, patent applications, and other publications, including all sequences disclosed within these references, referred to herein are expressly incorporated herein by reference, to the same extent as if each individual publication, patent or patent application was specifically and individually indicated to be incorporated by reference. All documents cited are, in relevant part, incorporated herein by reference in their entireties for the purposes indicated by the context of their citation herein. However, the citation of any document is not to be construed as an admission that it is prior art with respect to the present disclosure.
[0042] In the following detailed description, reference is made to the accompanying drawings, which form a part hereof. In the drawings, similar symbols typically identify similar components, unless context dictates otherwise. The illustrative embodiments described in the detailed description, drawings, and claims are not meant to be limiting. Other embodiments may be utilized, and other changes may be made, without departing from the spirit or scope of the subject matter presented herein. It will be readily understood that the aspects of the present disclosure, as generally described herein, and illustrated in the figures,can be arranged, substituted, combined, separated, and designed in a wide variety of different configurations, all of which are explicitly contemplated herein.
[0043] Sequencing Aided by Mutagenesis (SAM) is a technique used to aid sequence alignment by intentionally introducing mutations into nucleic acid molecules (for example, into about 5% of nucleotides during library preparation), and using these introduced mutations to help align sequence reads. For example, the Illumina Complete Long Read (ICLR) sequencing technology employs SAM to determine the sequences of synthetic long reads from short read sequence data. These long reads generated from short reads are referred to as “ICLRs.” The current biochemistry used in ICLR technology introduces transition mutations such as A to G and C to T. Some long-read sequencing technologies produce sequence reads which have a relatively high sequencing error rate, such as in about 5% or more of base calls or even more than about 10% of base calls. Sequencing errors may be a base swap where one base is converted to another, or instead can be an insertion or deletion of nucleotides. However, intentionally introduced mutations and sequencing errors may complicate sequence alignment and analysis because such mutations and errors may make it computationally difficult to align two different versions of a nucleic acid sequence, such as where some sequence reads include mutations and / or sequencing errors, and other sequence reads do not contain these mutations and / or sequencing errors.
[0044] The present disclosure relates to methods and systems for aligning nucleotide sequences that may be used, for example, on data from SAM (such as ICLR) and / or error-prone sequencing technologies. For example, the methods and systems described herein may receive sequence reads from nucleic acids and then employ a variety of techniques to improve the alignment of those sequences. For example, the sequences may be encoded with a reduced sequence representation, such as a reduced alphabet. The reduced alphabet may be selected based on the mutation (or sequencing error) biochemistry. For example, an A to G or Gto A conversion may be encoded with a single character such as “R”. Similarly, a C to T or T to C conversion may be encoded with a single character, such as “Y”. In some embodiments, use of the reduced sequence representation may allow for a more efficient alignment of a sequence read which includes a mutation and / or sequence error to a sequence read which does not include a mutation and / or sequence error, because both sequence reads may have the same encoding under the reduced sequence representation, thereby facilitating matching and creatingmore complete alignments. The methods and systems disclosed herein may then align any transformed sequence reads (sequence reads which have been encoded with the reduced sequence representation).
[0045] The present disclosure also relates to other techniques for efficient sequence alignment, including homopolymer quantization, where the length of homopolymers within sequence reads may be rounded to a predetermined value. Sequencing methods commonly err in determining a length of a homopolymer run. In some embodiments, a set of predetermined values may be created for homopolymer lengths, and when a homopolymer run is detected in a sequence read, its length may be approximated to a value from the set of predetermined values. This technique may facilitate an alignment that tolerates sequencing errors, while retaining information about homopolymer length.
[0046] Syncmers are a pre-defined subset of k-mers, defined by the presence of an .s-minimizer in the k-mer. The concept of a syncmer was first introduced in Edgar (2021). The concept of a syncmer was later generalized by Dutta et al. (2022), who refer to the generalization as Parameterized Syncmer Schemes. In addition to the lengths k and .s of the k- mer and s-minimizer, a Parameterized Syncmer Scheme defines a variable to specify the required position of the s-minimizer within the k-mer and allows specification of more than one s minimizer position within a syncmer. In Dutta et al. (2022) it was shown that the overlap of neighboring syncmers can be minimized by careful selection of the positions of the s-minimizers within the k-mer. When a single s-minimizer is used, the optimal position for the s-minimizer is in the middle of the k-mer.
[0047] In some embodiments, a sequence index is generated from transformed sequence reads, and seeds (strings of sequence information constructed based on sequence reads) from index entries are matched. One additional embodiment relates to the construction of a new type of seed referred to as a “syncstrobe” which is a collection of nearby k-mers that is constructed by selecting a pair of k-mers that are approximately a pre-determined window apart from each other in a nucleic acid sequence (“outer pair of k-mers”), and selecting one or more additional k-mers between the outer pair of k-mers.
[0048] Additional embodiments may relate to specific techniques for selecting k- mers when constructing a syncstrobe. For example, an additional embodiment includes selecting an outer pair of k-mers regardless of which strand of a double-stranded nucleic acidthe k-mers correspond to (“strand-symmetric minimizer pairs”). Another embodiment may include using a function to select k-mers that determines whether syncmers (another type of seed) are available in a desired window, and if a syncmer is available, selecting a syncmer for use as a k-mer in the construction of the syncstrobe (“syncmer-preferred minimizer”).
[0049] Another embodiment may include local repeat marking, where seeds such as syncstrobes which are locally repeated within a nucleic acid sequence are marked as locally repetitive, and only a single representative seed is used in the sequence index, thereby improving alignment efficiency.
[0050] A final embodiment may include generating a set of syncstrobes for use as index entries in a sequence index, which is a data structure used for storing, sorting, and scanning strings of sequence information derived from sequence reads. For example, to generate sequence index entries, in some embodiments a syncstrobe is constructed as described above, one or more k-mer components is removed to generate a new syncstrobe, and the removing step is repeated by removing different k-mer components, to generate a set of syncstrobes which include different combinations of k-mers (“strobe drop indexing”). These and other innovations disclosed herein may further contribute to improved alignment efficiency and mutation recall.
[0051] In some embodiments, encoding sequences in a reduced sequence representation such as RY-encoding may increase sensitivity for sequence matching because it enables matching sequences with or without a transversion mutation. Other techniques discussed herein relate to preserving specificity when matching sequences that are encoded in a reduced sequence representation. For example, a reduced sequence representation alone may, in some embodiments, decrease specificity in sequence matching because sequences with a different nucleotide sequence (not due to transversion mutations) may have the same encoding with a reduced alphabet, producing false matches. Thus, the present disclosure also includes several techniques which enable efficient filtration of sequence matches, including the syncstrobe seeds discussed above and other techniques discussed further herein.
[0052] Sequence alignments produced by the methods and systems disclosed herein may be useful when comparing mutated sequences to unmutated sequences, for example in processes that aim to identify sites that contain introduced mutations via comparison amongst mutated and unmutated sequences. The methods and systems described herein mayfind particular utility in the field of Sequencing Aided by Mutagenesis (SAM). They may also be useful in error-prone sequencing methods, such as long-read sequencing.
[0053] In some embodiments, the disclosed methods and systems may be used with sequence reads that have a transition mutation at up to 100% of sites. In some embodiments, the disclosed methods and systems provide for aligning any matching sequence reads that have 1 or fewer transversion mutations. In some embodiments, the methods and systems can advantageously identify seed matches in the presence of multiple transversions and / or indel differences between sequence reads.
[0054] The methods and systems described herein have several advantages. First, they improve the accuracy and recall in identifying mutated sites or sequencing errors. In the context of ICLR technology, for example, these improvements may lead to better quality of ICLRs and eventually in better downstream secondary analysis like variant calling.
[0055] Second, the methods and systems described herein do not require a reference sequence. Some previous methods rely on a high-quality reference genome for the organism being sequenced, which is disadvantageous in the case of non-modal organisms or in applications like metagenomics. Some previous methods attempt to solve this problem by using de-novo assemblies from short unmutated reads as a reference in such cases, but this is computationally intensive. Further, the low-quality reference decreases the quality of the ICLRs produced by the alignments and thus the downstream secondary analysis. In contrast, the proposed methods and systems may be “reference-free” in that no reference may be required for identifying mutations in reads (though a reference may be used if available).
[0056] Third, the methods and systems disclosed herein may tolerate increased mutation rates, such as transition mutations at up to 100% of sites. This feature may be particularly advantageous as increasing the mutation rate in sequencing aided by mutagenesis (SAM) techniques is expected to help the assembly stage, and may yield more accurate and complete ICLRs, thus further improving variant calling. In order to take advantage of SAM techniques with high mutation rates, methods and systems that are more tolerant to the increased mutations are needed. The methods and systems described herein are tolerant of increased mutation rates, and hence can support the potential gain in sequence assembly and variant calling with an increased mutation rate.
[0057] The methods and systems described herein are advantageous because in some embodiments they enable long chains of minimizers to efficiently filter sequence matches with high sensitivity for short sequences even in the presence of high mutation rates. The concept of an index on values constructed from collections of nearby k-mers (known as “strobemers”) was previously introduced by Sahlin (2021). However, strobemers were introduced without the advantages of strand-symmetry and syncmer-preferred minimizer strobes that are described herein. In the work of Sahlin (2021), only simple minimizers (minstrobes) and a random selection function (randstrobes) were described for use with strobemer construction. As such, the method of Sahlin (2021) would typically require every position of a sequence to be indexed, with each strand being separately indexed. In some embodiments, the methods described herein can produce a much smaller index because the combination of strand-symmetry and syncmer-preferred minimizer strobes means only a small fraction of all positions become indexed, yet sensitivity and specificity are preserved. For example, while Sahlin (2021) introduced the concept of strobemers, the present disclosure introduces several advantages, including strand-symmetric minimizer pairs, syncmer-preferred minimizers, homopolymer quantization, reduced-alphabet encoded strobemers, local repeat marking, and strobe drop indexing.
[0058] Some prior work, such as Blassel et al. (2022), generally discussed a homopolymer compression technique in sequence analysis which may compress a homopolymer run to a single base representation of that nucleotide. However, in contrast, in some embodiments, the methods and systems described herein instead quantize (round) the length of the homopolymer to a nearby value. When homopolymer lengths are approximately correct, the quantization approach has the significant advantage of preserving information about the length of a homopolymer run, while still tolerating error in the run length.
[0059] Furthermore, some embodiments have additional advantages because they may readily generalize to support matches containing indel differences.Definitions
[0060] Although the following terms are believed to be well understood by one of skill in the art, the following definitions are set forth to facilitate understanding of the presently disclosed subject matter.
[0061] All technical and scientific terms used herein, unless otherwise defined below, are intended to have the same meaning as commonly understood by one of ordinary skill in the art. References to techniques employed herein are intended to refer to the techniques as commonly understood in the art, including variations on those techniques or substitutions of equivalent techniques that would be apparent to one of skill in the art.
[0062] As used herein, the terms “a” or “an” or “the” may refer to one or more than one. For example, “a” marker can mean one marker or a plurality of markers.
[0063] As used herein, the term “about,” when used in reference to a measurable value such as an amount of mass, dose, time, temperature, and the like, is meant to encompass variations of 20%, 10%, 5%, 1%, 0.5%, or even 0.1% of the specified amount.
[0064] As used herein, the term “and / or” refers to and encompasses any and all possible combinations of one or more of the associated listed items, as well as the lack of combinations when interpreted in the alternative (“or”).
[0065] Throughout this specification, unless the context requires otherwise, the words “comprise,” “comprises,” and “comprising” will be understood to imply the inclusion of a stated step or element or group of steps or elements but not the exclusion of any other step or element or group of steps or elements.
[0066] As used herein, the term “consists essentially of’ (and grammatical variants thereof), as applied to the compositions and methods of the present disclosure, means that the compositions / methods may contain additional components so long as the additional components do not materially alter the composition / method.
[0067] The term “nucleic acid” or “polynucleotide” refers to a deoxyribonucleotide or ribonucleotide polymer in either single- or double-stranded form, and unless otherwise limited, encompasses known analogs of natural nucleotides that hybridize to nucleic acids in manner similar to naturally occurring nucleotides, such as peptide nucleic acids (PNAs) and phosphorothioate DNA. Unless otherwise indicated, a particular nucleic acid sequence includes the complementary sequence thereof. Nucleotides include, but are not limited to, ATP, dATP, CTP, dCTP, GTP, dGTP, UTP, TTP, dUTP, 5-methyl-CTP, 5-methyl-dCTP, ITP, diTP, 2-amino-adenosine-TP, 2-amino-deoxyadenosine-TP, 2-thiothymidine triphosphate, pyrrolo-pyrimidine triphosphate, and 2-thiocytidine, as well as the alphathiotriphosphates for all of the above, and 2'-O-methyl-ribonucleotide triphosphates for all the above bases.Modified bases include, but are not limited to, 5-Br-UTP, 5-Br-dUTP, 5-F-UTP, 5-F-dUTP, 5-propynyl dCTP, and 5-propynyl-dUTP.
[0068] As used herein, the term “reference genome” or “reference sequence” refers to any particular known genome sequence, whether partial or complete, of any organism or virus which may be used to reference identified sequences from a subject. For example, a reference genome used for human subjects as well as many other organisms is found at the National Center for Biotechnology Information at ncbi.nlm.nih.gov. In various embodiments, the reference sequence is significantly larger than the reads that are aligned to it. For example, it may be at least about 100 times larger, or at least about 1000 times larger, or at least about 10,000 times larger, or at least about 105times larger, or at least about 106times larger, or at least about 107times larger. In one example, the reference sequence is that of a full-length genome. Such sequences may be referred to as genomic reference sequences. Other examples of reference sequences include genomes of other species, such as of control organisms as disclosed herein, as well as chromosomes, sub-chromosomal regions (such as strands), etc., of any species. In various embodiments, the reference sequence is a consensus sequence or other combination derived from multiple individuals. However, in certain applications, the reference sequence may be taken from a particular individual.
[0069] The term “nucleic acid sample” herein may refer to a sample, typically derived from one or more biological fluids, cells, tissues, organs, or organisms, comprising a nucleic acid or a mixture of nucleic acids comprising at least one nucleic acid sequence that is to be screened for copy number variation. In certain embodiments the nucleic acid sample comprises at least one nucleic acid sequence whose copy number is suspected of having undergone variation. Such samples may include, but are not limited to sputum / oral fluid, amniotic fluid, blood, a blood fraction, or fine needle biopsy samples (such as surgical biopsy, fine needle biopsy, etc.), urine, peritoneal fluid, pleural fluid, and the like. Although the sample is often taken from a human subject (such as a patient), the sample may be from any mammal, including, but not limited to dogs, cats, horses, goats, sheep, cattle, pigs, etc. The sample may be used directly as obtained from the biological source or following a pretreatment to modify the character of the sample. For example, such pretreatment may include preparing plasma from blood, diluting viscous fluids and so forth. Methods of pretreatment may also involve, but are not limited to, filtration, precipitation, dilution, distillation, mixing, centrifugation,freezing, lyophilization, concentration, amplification, nucleic acid fragmentation, inactivation of interfering components, the addition of reagents, lysing, etc. If such methods of pretreatment are employed with respect to the sample, such pretreatment methods are typically such that the nucleic acid(s) of interest remain in the test sample, sometimes at a concentration proportional to that in an untreated test sample (such as namely, a sample that is not subjected to any such pretreatment method(s)). Such “treated” or “processed” samples are still considered to be biological “test” samples with respect to the methods described herein. A “nucleic acid sample” may also include nucleic acid sequence information stored in a memory, and which was originally obtained from a source such as one or more biological fluids, cells, tissues, organs, or organisms.
[0070] The term “read” or “sequence read” (or sequencing reads) refers to a sequence obtained from a portion of a nucleic acid sample. A read may be represented by a string of nucleotides sequenced from any part or all of a nucleic acid molecule. Typically, though not necessarily, a read represents a short sequence of contiguous base pairs in the sample. The read may be represented symbolically by the base pair sequence (in A, T, C, or G) of the sample portion. It may be stored in a memory device and processed as appropriate to determine whether it matches a reference sequence or meets other criteria. A read may be obtained directly from a sequencing apparatus or indirectly from stored sequence information concerning the sample. In some cases, a read is a DNA sequence of sufficient length (such as at least about 25 bp) that can be used to identify a larger sequence or region, for example, that can be aligned and specifically assigned to a chromosome or genomic region or gene. For example, a sequence read may be a short string of nucleotides (such as 20-150 bases) sequenced from a nucleic acid fragment, a short string of nucleotides at one or both ends of a nucleic acid fragment, or the sequencing of the entire nucleic acid fragment that exists in the biological sample. Sequence reads may be obtained by any method known in the art. For example, a sequence read may be obtained in a variety of ways, such as using sequencing techniques or using probes, such as in hybridization arrays or capture probes, or amplification techniques, such as the polymerase chain reaction (PCR) or linear amplification using a single primer or isothermal amplification. Sequence reads can be generated by techniques such as sequencing by synthesis, sequencing by binding, or sequencing by ligation. Sequence readscan be generated using instruments such as MINISEQ, MISEQ, NEXTSEQ, HISEQ, and NOVASEQ sequencing instruments from Illumina, Inc. (San Diego, CA).
[0071] As used herein, a “short sequence read” refers to a sequence read of between 50-500 bp, for example, about 50 - 300 bp, and includes paired end sequence reads.
[0072] As used herein, a “long sequence read” refers to a sequence read of more than about 500 bp, for example 500 - 250,000 bp or more. A long sequence read may be obtained from a long-read sequencing technology, or may be synthetically constructed from multiple short sequence reads (for example, ICLRs).
[0073] As used herein, a “transformed sequence read” refers to a sequence read which has been encoded with a reduced sequence representation, for example RY encoding or WS encoding.
[0074] ‘RY encoding” or “RY format” refers to encoding a nucleic acid sequence with only two letters, where “R” stands for purines (adenine [A] and guanine [G]), and “Y” stands for pyrimidines (cytosine [C] and thymine [T] in DNA, or uracil [U] in RNA).
[0075] As used herein, a “file” includes digital files. In some embodiments, a file is on a computer storage medium (such as a computer hard drive, for example a spinning magnetic disk drive or a solid state drive). In some embodiments, the digital file is stored in the format of a BAM, FASTQ, SAM, CRAM, JSON, CIGAR, or VCF file.
[0076] As used herein, a k-mer refers to a sequence of length k.
[0077] As used herein, a “minimizer” refers to a type of k-mer that is selected by a minimizer function from a set of all possible k-mers within a given window. A minimizer may be selected by dividing a nucleic acid sequence into overlapping windows of a fixed length, and selecting a “smallest” k-mer within the window from the set of all possible k-mers. The term “smallest” can refer to the lexicographic order, a numeric hash value, or some other predetermined criterion. Minimizers reduce computational burden, as compared to using all possible k-mers of a sequence, because only one k-mer is selected per window. Minimizers are described in [Roberts et al., Reducing storage requirements for biological sequence comparison, Bioinformatics, 20: 18, 3363-3369 (2004)].
[0078] As used herein, a “syncmer” refers to a k-mer that includes a shorter s-mer minimizer at a position from a set of one or more pre-specified positions. A [k, s, {«}]-syncmeris a A mer that contains a shorter 5-mer minimizer at a position of the ir-mer that is in the set {«}•
[0079] As used herein, a “strobemer” refers to a collection of nearby k-mers.
[0080] As used herein, the term “syncstrobe” refers to a strobemer that includes syncmer-preferred minimizers. For example, a syncstrobe may be a strobemer generated by selecting an outer pair of k-mers and one or more k-mers between the outer pair, where syncmers are preferred for selection as the k-mers.Methods for Sequence Alignment
[0081] In an aspect, disclosed herein are methods for aligning sequence reads based on reduced sequence representations. FIG. 1 is a flow diagram that schematically illustrates an exemplary method 100 for aligning sequence reads based on reduced sequence representations. In some embodiments, the method 100 is implemented on a computer. The method 100 may be embodied in a set of executable program instructions stored on a computer-readable medium, such as one or more disk drives, of a computing system. When the method 100 is initiated, the executable program instructions can be loaded into a memory and executed by one or more processors of a server device.Receiving sequence reads
[0082] As shown in FIG. 1, the method 100 for aligning sequence reads based on reduced sequence representations may start from start block 110. The method 100 may proceed to block 120, wherein sequence reads from nucleic acids are received. The sequence reads may be generated using any method known to those of skill in the art. The sequence reads may be in any digital file format.
[0083] In some embodiments, the sequence reads comprise mutations or sequence errors. In some embodiments, at least 4% of the nucleotides comprise a mutation or a sequence error. In some embodiments, the mutation comprises a transition mutation.
[0084] For example, in some embodiments, the sequence reads were generated through an intentional mutation process that introduces random mutations (such as transition mutations) in nucleic acid molecules, such as a Sequencing Aided by Mutagenesis (SAM) methods. In some embodiments, transition mutations are introduced in about 1 %- 10% of bases,for example about 5% of bases. In some embodiments, the nucleic acids comprise mutated nucleic acids and unmutated nucleic acids, and the sequence reads comprise mutated sequence reads and unmutated sequence reads.
[0085] In addition, the method 100 may also be used in other applications that do not necessarily involve SAM, such as in long-read sequencing, which may be prone to sequencing errors.Encoding sequence reads with a reduced sequence representation
[0086] The method 100 may proceed to block 130, wherein sequence reads are encoded with a reduced sequence representation, thereby generating transformed sequence reads.
[0087] In some embodiments, sequence reads are encoded using a reduced alphabet. In some embodiments, the method 100 includes encoding the sequence reads in RY format (“RY-encoding”). For example, RY-encoding may advantageously allow for recording and processing of sequence reads in a format that is tailored to the relevant biochemistry used in a SAM process involving transition mutations, where an unmutated base and a mutated base would be recorded with the same character. In such an embodiment, a transition mutation would effectively produce no change in the recorded sequence using the reduced alphabet. Similarly, for sequencing techniques which introduce accidental sequence errors, the sequencing errors may be recorded via a reduced alphabet in a way that does not change the sequence. This may allow for better alignment and / or comparison between unmutated sequence reads and mutated sequence reads (and similarly between sequence reads which include a sequencing error and sequence reads which do not).
[0088] In some embodiments, the method 100 includes WS-encoding the sequence reads. In some embodiments, reduced-alphabet representations of sequence reads are stored in a digital file.
[0089] In some embodiments, the length of homopolymers is reduced within the transformed sequence reads. In some embodiments, homopolymers are quantized. For example, in some embodiments, the length of homopolymers within the sequence reads is rounded to a predetermined value of a set of two or more predetermined values. In someembodiments, this step of reducing length of homopolymers takes place before encoding sequence reads with a reduced alphabet.Aligning transformed sequence reads
[0090] The method 100 may proceed to process block 140, wherein transformed sequence reads are aligned. The process block 140 is explained in more detail with reference to FIG. 2A, which is a block diagram that illustrates additional details on the process taking place within the process block 140.
[0091] The method 100 may be used to align sequence reads of various lengths to each other, optionally without a reference sequence. For example, in some embodiments, short sequence reads are aligned to each other. In some embodiments, short sequence reads to long sequence reads, such as ICLRs. In some embodiments, the short sequence reads are aligned to a reference sequence, such as the human genomic sequence GRCh38 from the Genome Reference Consortium.
[0092] In some embodiments, aligning the transformed sequence reads comprises aligning without a reference sequence. For example, the method 100 may advantageously still perform an alignment when a reference sequence is not available (de novo sequencing) or is of low quality.
[0093] The method 100 may end at block 150.
[0094] Alignments of transformed sequence reads may be used for further processing and sequence analysis (not shown). For example, the method 100 may further include identifying mutated positions in mutated sequence reads by comparing a mutated sequence read set and an unmutated read sequence set.Generating a sequence index from transformed sequence reads
[0095] FIG. 2A is a block diagram that illustrates additional details on the specific methods taking place within the process block 140, wherein transformed sequence reads are aligned. As shown in FIG. 2 A, the process 140 may include process block 210, wherein a sequence index is generated from transformed sequence reads. At process block 210, the method may include processing transformed sequence reads to generate index entries, and sorting index entries. The process block 210 is explained in more detail with reference to FIG.2B, which is a block diagram that illustrates additional details on the process taking place within the process block 210.Identifying one or more seed matches
[0096] The process 140 may proceed to block 220, wherein one or more seed matches are identified from index entries. For example, seed matches may be identified by sorting and scanning indexes of the sequence index.
[0097] In some embodiments, each seed match triggers a processing step. In some embodiments where sequence reads include unmutated sequence reads and mutated sequence reads, the processing step could collect all matching unmutated sequence reads for each mutated read and store only top N matches (where N is a predetermined number) for use in generating alignments as described below with reference to block 230. In some embodiments, a set of N highest-scoring unmutated reads (where N is a predetermined number) may be selected for each mutated sequence read. In some embodiments, the method comprises selecting a set of unmutated sequence reads which are consistent with each other, which have a sequence identity above a threshold to a mutated sequence read covering a corresponding sequence, and / or which together cover at least a predetermined percentage of the mutated sequence read. In some embodiments, the method comprises, based on a comparison of one or more alignments of unmutated sequence reads and one or more alignments of mutated sequence reads, marking mutated positions at sites which differ between the mutated sequence read set and the unmutated sequence read set.Generating initial alignments
[0098] The process 140 may proceed to block 230, wherein the initial alignments are generated based on the one or more seed matches. In some embodiments, the initial alignment is an ungapped alignment.Scoring alignments
[0099] The process 140 may proceed to block 240, wherein alignments (such as the initial alignments generated at block 230) are scored. In some embodiments, scoring comprises using a scoring scheme which accounts for transition mutations. In some embodiments, the scoring comprises bit-parallel alignment scoring. For example, at block 240, genomic positionsof k-mer components of a strobemer seed match may be used to guide ungapped bit-parallel alignment scoring.
[0100] The process 140 may proceed to decision state 250, wherein the methods and systems determine whether an alignment score is above a threshold. If the score is above a threshold, the process 140 may proceed to block 260, wherein the seed match (and / or an associate alignment) is stored.
[0101] If the score is not above a threshold, the process 140 may proceed to block 265, wherein an alignment may be modified, an alternative alignment may be computed, or a seed match may be discarded. For example, if an initial alignment is an ungapped alignment, an alternative alignment such as a soft-clipped or (partially) gapped alignment may be computed. For example, in some embodiments, at block 265, k-mer components of a strobemer seed match may be used in a gapped sequence alignment. In the case of where an alignment is modified or an alternative alignment is computed, the process 140 may return to block 240 and the process 140 may proceed as previously described. Seed matches which do not produce alignments over a scoring threshold may ultimately be discarded.Computing further sequence alignments
[0102] The process 140 may proceed to block 270, wherein further sequence alignments are computed based on stored seed matches. For example, alignments may be combined to form a paired-end alignment. Chaining and extending from seeds may also proceed.Processing Transformed Sequence Reads to Generate Index Entries
[0103] FIG. 2B is a block diagram that illustrates the details on the methods taking place within the process block 210, wherein a sequence index is generated from transformed sequence reads. In some embodiments, the index entries comprise syncstrobes. However, in alternative embodiments, other types of seeds for use as index entries may be generated, such as other types of strobemers.
[0104] As shown in FIG. 2B, the process 210 may include block 215, wherein an outer pair of k-mers and one or more additional k-mers between the outer pair of k-mers are selected. In some embodiments, a minimizer function is used to select k-mers. In someembodiments, the minimizer function evaluates k-mers based on a parameterized syncmer scheme. For example, in some embodiments, the minimizer function evaluates k-mers based on whether s-mer minimizers occur in one or more positions specified in the parameterized syncmer scheme.
[0105] For example, in some embodiments, k-mers and s-mers are computed from transformed sequence reads, and syncmers are identified. A minimizer function may be used to select a pair of k-mers that are approximately a pre-determined distance apart from each other. In some embodiments, the minimizer function will determine whether there are any syncmers selects the minimum among the syncmers present in the window, but if no syncmers are present then the minimizer function selects a minimum among the k-mers (a minimizer). In some embodiments, the minimizer function is strand- symmetric and may use k-mers that are from either strand of a double-stranded nucleic acid sequence.
[0106] Furthermore, at block 215, locally repetitive strobemers (such as syncstrobes) may be identified, and a single representative strobemer (such as a syncstrobe) may be stored to be used for seeding. This step may further improve alignment efficiency and speed. In some embodiments, an index entry (such as a strobemer / syncstrobe) may be marked as locally repetitive.
[0107] The process 210 may proceed to block 217, wherein one or more k-mer components are removed (dropped) from syncstrobes, thereby generating sets of syncstrobes. In some embodiments, the one or more individual k-mer components are removed to improve one or more of sensitivity, specificity, or index compression. For example, generating sets of syncstrobes may reduce the size of the index compared to using all possible k-mers. In some embodiments, this step block 217 may improve alignment processes because it allows multiple version of syncstrobes to be created in a set, and therefore if one k-mer component includes an error that would cause problems with alignments, there will be other syncstrobes in the set which do not include that component, thereby increasing sensitivity.
[0108] The process 140 may include block 219, wherein index entries are ordered. For example, sets of index entries such as syncstrobes may be ordered, for example, using a hash function to select a representative for each set. In some embodiments, index entries are ordered by representative and a top N number of sets, wherein N is a predetermined value, are chosen to be stored to a sequence index file.Systems for Sequence Alignment
[0109] Further disclosed herein are electronic systems for aligning sequence reads based on reduced sequence representations. In some embodiments, the system includes a processor configured to perform a method comprising: receiving sequence reads from nucleic acids, encoding sequence reads from nucleic acids with a reduced sequence representation, thereby generating transformed sequence reads, and aligning the transformed sequence reads.
[0110] Further disclosed herein are non-transitory computer-readable media. In some embodiments, the non-transitory computer-readable medium includes a plurality of instructions, which when executed by at least one processor, cause the at least one processor to: receive sequence reads from nucleic acids, encode sequence reads from nucleic acids with a reduced sequence representation, thereby generating transformed sequence reads, and align the transformed sequence reads.
[0111] FIG. 3 A illustrates a diagram of an environment in which a sequence alignment system can operate in accordance with one or more implementations. The following paragraphs describe the sequence alignment system with respect to illustrative figures that portray example implementations and embodiments. For example, FIG. 3A illustrates a schematic diagram of a computing system 3000 in which a sequence alignment application 3106 operates in accordance with one or more implementations. As illustrated, the computing system 3000 includes one or more server device(s) 3102 connected to a user client device 3108, a local device 3118, and a sequencing device 3114 via a network 3112. The network 3112 can comprise any suitable network over which computing devices can communicate.
[0112] As shown in FIG. 3 A, the computing system 3000 includes the server device(s) 3102. In various implementations, the server device(s) 3102 may generate, receive, analyze, store, and transmit digital data, such as data for nucleobase calls or sequenced nucleic- acid polymers. In some implementations, the server device(s) 3102 receive various data from the sequencing device 3114, such as data from a sample genome and / or sequence reads. The server device(s) 3102 may also communicate with the user client device 3108. In particular, the server device(s) 3102 can send data for sequence reads, direct nucleobase calls, nucleobase calls, and / or sequencing metrics to the user client device 3108.
[0113] As shown, the server device(s) 3102 includes a sequencing application 3110. In general, the sequencing application 3110 analyzes the data (such as call data) receivedfrom the sequencing device 3114 or elsewhere to determine nucleobase sequences for nucleic- acid polymers. For example, the sequencing application 3110 can receive raw data from the sequencing device 3114 and determine a nucleobase sequence for a sample genome or a nucleic-acid segment. In some implementations, the sequencing application 3110 determines the sequences of nucleobases in DNA and / or RNA segments or oligonucleotides.
[0114] As also shown, the sequencing application 3110 includes the sequence alignment application 3106. As described below, in some embodiments, the sequence alignment application 3106 can align sequence reads based on reduced sequence representations. For example, in some embodiments, the sequence alignment application 3106 receives sequence reads from nucleic acids, encodes sequence reads from nucleic acids with a reduced sequence representation, thereby generating transformed sequence reads, and aligns the transformed sequence reads.
[0115] While the sequencing application 3110 has been described as including the sequence alignment application 3106, other systems or methods may be included within the sequencing application 3110, such as an application to determine taxonomic groups associated with sample nucleic acids (not illustrated).
[0116] Moreover, while the sequence alignment application 3106 is described being implemented on the server device(s) 3102, as part of the sequencing application 3110, in some implementations, the sequence alignment application 3106 is implemented by (such as located entirely or in part) on the user client device 3108, the sequencing device 3114, and / or the local device 3118. As mentioned, in some implementations, the sequence alignment application 3106 is implemented by one or more other components of the computing system 3000, such as the sequencing device 3114. In particular, the sequence alignment application 3106 can be implemented in a variety of different ways across the server device(s) 3102, the network 3112, the user client device 3108, the local device 3118, and the sequencing device 3114.
[0117] As further shown in FIG. 3 A, the computing system 3000 includes the user client device 3108. In various implementations, the user client device 3108 can generate, store, receive, and send digital data. In particular, the user client device 3108 can receive the data from the sequencing device 3114. As further illustrated, the user client device 3108 includes a sequencing application 3110. The sequencing application 3110 may be a web application or anative application stored and executed on the user client device 3108 (for example, a mobile application, desktop application, or web application). The sequencing application 3110 can receive data from the sequencing application 3110 and / or sequence alignment application 3106. For example, the user client device 3108 can receive variant call files and / or alignment files from the sequencing application 3110.
[0118] The sequencing application 3110 can also include instructions that (when executed) cause the user client device 3108 to receive data from the sequence alignment application 3106 and present data from the sequencing device 3114 and / or the server device(s) 3102. Furthermore, the sequencing application 3110 can instruct the user client device 3108 to display data for variant calls, such as nucleobase calls or an indication of a copy number variant. Indeed, the user client device 3108 can display nucleobase call results for a genome sample and / or an indication of a predicted copy number variant.
[0119] As further shown in FIG. 3A, the computing system 3000 includes the sequencing device 3114. In various implementations, the sequencing device 3114 can sequence a genomic sample or other nucleic-acid polymer. For example, the sequencing device 3114 analyzes nucleic-acid segments or oligonucleotides extracted from genomic samples to generate data either directly or indirectly on the sequencing device 3114. More particularly, the sequencing device 3114 receives and analyzes, within nucleotide-sample slides (such as flow cells), nucleic-acid sequences extracted from genomic samples. In one or more implementations, the sequencing device 3114 utilizes sequencing by synthesis (SBS) to sequence a genomic sample or other nucleic-acid polymers. In addition to, or in the alternative, to communicating across the network 3112, in some implementations, the sequencing device 3114 bypasses the network 3112 and communicates directly with the user client device 3108.
[0120] As further depicted in FIG. 3A, in some implementations, the server device(s) 3102 includes a distributed collection of servers, where the server device(s) 3102 include several server devices distributed across the network 3112 and located in the same or different physical locations. For instance, the server device(s) 3102 can be implemented, in whole or in part, on the local device 3118. To illustrate, the local device 3118 may implement the sequencing application 3110 and / or the sequence alignment application 3106. Further, the server device(s) 3102 and / or the local device 3118 can include a content server, an application server, a communication server, a web-hosting server, or another type of server.
[0121] The user client device 3108 illustrated in FIG. 3 A can include various types of client devices. For example, in some implementations, the user client device 3108 includes non-mobile devices, such as desktop computers or servers, or other types of client devices. In various implementations, the user client device 3108 includes mobile devices, such as laptops, tablets, mobile telephones, or smartphones.
[0122] Though FIG. 3A illustrates the components of the computing system 3000 communicating via the network 3112, in certain implementations, the components of computing system 3000 can also communicate directly with each other, bypassing the network 3112. For instance, in some implementations, the user client device 3108 communicates directly with the sequencing device 3114. Additionally, in some implementations, the user client device 3108 communicates directly with the sequence alignment application 3106 and / or the server device(s) 3102. In some implementations, the user client device 3108 communicates directly with the local device 3118. Moreover, the sequence alignment application 3106 can access one or more databases housed on or accessed by the server device(s) 3102 or elsewhere in the computing system 3000.
[0123] FIG. 3B is a block diagram of an exemplary server device 3102 that may be used in connection with the computing system 3000 of FIG. 3 A. The server device 3102 may be configured to align sequence reads based on reduced sequence representations. The general architecture of the server device 3102 depicted in FIG. 3B includes an arrangement of computer hardware and software components. The server device 3102 may include many more (or fewer) elements than those shown in FIG. 3B. It is not necessary, however, that all of these generally conventional elements be shown in order to provide an enabling disclosure. As illustrated, the server device 3102 includes a processing unit 310, a network interface 320, a computer readable medium drive 330, an input / output device interface 340, a display 350, and an input device 360, all of which may communicate with one another by way of a communication bus. The network interface 320 may provide connectivity to one or more networks or computing systems. The processing unit 310 may thus receive information and instructions from other computing systems or services via a network. The processing unit 310 may also communicate to and from memory 370 and further provide output information for an optional display 350 via the input / output device interface 340. The input / output device interface 340 may also accept input from the optional input device 360, such as a keyboard,mouse, digital pen, microphone, touch screen, gesture recognition system, voice recognition system, gamepad, accelerometer, gyroscope, or other input device.
[0124] The memory 370 may contain computer program instructions (grouped as modules or components in some embodiments) that the processing unit 310 executes in order to implement one or more embodiments. The memory 370 generally includes RAM, ROM and / or other persistent, auxiliary or non-transitory computer readable media. The memory 370 may store an operating system 372 that provides computer program instructions for use by the processing unit 310 in the general administration and operation of the server device 3102. The memory 370 may store a reference genome 373, such as for use by the sequencing application 3110. The memory 370 may further include computer program instructions and other information for implementing aspects of the present disclosure.
[0125] For example, in one embodiment, the memory 370 includes a sequencing application 3110, which may include a sequence alignment application 3106. The sequence alignment application 3106 can perform the methods disclosed herein. In addition, memory 370 may include or communicate with the data store 390 and / or one or more other data stores that store one or more inputs, one or more outputs, and / or one or more results (including intermediate results) of aligning sequence reads based on reduced sequence representations,, such the sequence reads, transformed sequence reads, seeds, alignments, and / or one or more reference genomes.
[0126] In some embodiments, the disclosed systems and methods may involve approaches for shifting or distributing certain sequence data analysis features and sequence data storage to a cloud computing environment or cloud-based network. User interaction with sequencing data, genome data, or other types of biological data may be mediated via a central hub that stores and controls access to various interactions with the data. In some embodiments, the cloud computing environment may also provide sharing of protocols, analysis methods, libraries, sequence data as well as distributed processing for sequencing, analysis, and reporting. In some embodiments, the cloud computing environment facilitates modification or annotation of sequence data by users. In some embodiments, the systems and methods may be implemented in a computer browser, on-demand or on-line.
[0127] In some embodiments, software written to perform the methods as described herein is stored in some form of computer readable medium, such as memory, CD-ROM, DVD-ROM, memory stick, flash drive, hard drive, SSD hard drive, server, mainframe storage system and the like.
[0128] In some embodiments, the methods may be written in any of various suitable programming languages, for example compiled languages such as C, C#, C++, Fortran, and Java. Other programming languages could be script languages, such as Perl, MatLab, SAS, SPSS, Python, Ruby, Pascal, Delphi, R and PHP. In some embodiments, the methods are written in C, C#, C++, Fortran, Java, Perl, R, Java or Python. In some embodiments, the method may be an independent application with data input and data display modules. Alternatively, the method may be a computer software product and may include classes wherein distributed objects comprise applications including computational methods as described herein.
[0129] In some embodiments, the methods may be incorporated into pre-existing data analysis software, such as that found on sequencing instruments. Software comprising computer implemented methods as described herein are installed either onto a computer system directly, or are indirectly held on a computer readable medium and loaded as needed onto a computer system. Further, the methods may be located on computers that are remote to where the data is being produced, such as software found on servers and the like that are maintained in another location relative to where the data is being produced, such as that provided by a third party service provider.
[0130] An assay instrument, desktop computer, laptop computer, or server which may contain a processor in operational communication with accessible memory comprising instructions for implementation of systems and methods. In some embodiments, a desktop computer or a laptop computer is in operational communication with one or more computer readable storage media or devices and / or outputting devices. An assay instrument, desktop computer and a laptop computer may operate under a number of different computer based operational languages, such as those utilized by Apple-based computer systems or PC based computer systems. An assay instrument, desktop and / or laptop computers and / or server system may further provide a computer interface for creating or modifying experimental definitions and / or conditions, viewing data results and monitoring experimental progress. In some embodiments, an outputting device may be a graphic user interface such as a computer monitor or a computer screen, a printer, a hand-held device such as a personal digital assistant (such asPDA, Blackberry, iPhone), a tablet computer (such as iP D), a hard drive, a server, a memory stick, a flash drive and the like.
[0131] A computer readable storage device or medium may be any device such as a server, a mainframe, a supercomputer, a magnetic tape system and the like. In some embodiments, a storage device may be located onsite in a location proximate to the assay instrument, for example adjacent to or in close proximity to, an assay instrument. For example, a storage device may be located in the same room, in the same building, in an adjacent building, on the same floor in a building, on different floors in a building, etc. in relation to the assay instrument. In some embodiments, a storage device may be located off-site, or distal, to the assay instrument. For example, a storage device may be located in a different part of a city, in a different city, in a different state, in a different country, etc. relative to the assay instrument. In embodiments where a storage device is located distal to the assay instrument, communication between the assay instrument and one or more of a desktop, laptop, or server is typically via Internet connection, either wireless or by a network cable through an access point. In some embodiments, a storage device may be maintained and managed by the individual or entity directly associated with an assay instrument, whereas in other embodiments a storage device may be maintained and managed by a third party, typically at a distal location to the individual or entity associated with an assay instrument. In embodiments as described herein, an outputting device may be any device for visualizing data.
[0132] An assay instrument, desktop, laptop and / or server system may be used itself to store and / or retrieve computer implemented software programs incorporating computer code for performing and implementing computational methods as described herein, data for use in the implementation of the computational methods, and the like. One or more of an assay instrument, desktop, laptop and / or server may comprise one or more computer readable storage media for storing and / or retrieving software programs incorporating computer code for performing and implementing computational methods as described herein, data for use in the implementation of the computational methods, and the like. Computer readable storage media may include, but is not limited to, one or more of a hard drive, a SSD hard drive, a CD-ROM drive, a DVD-ROM drive, a floppy disk, a tape, a flash memory stick or card, and the like. Further, a network including the Internet may be the computer readable storage media. In some embodiments, computer readable storage media refers to computational resourcestorage accessible by a computer network via the Internet or a company network offered by a service provider rather than, for example, from a local desktop or laptop computer at a distal location to the assay instrument.
[0133] In some embodiments, computer readable storage media for storing and / or retrieving computer implemented software programs incorporating computer code for performing and implementing computational methods as described herein, data for use in the implementation of the computational methods, and the like, is operated and maintained by a service provider in operational communication with an assay instrument, desktop, laptop and / or server system via an Internet connection or network connection.
[0134] In some embodiments, a hardware platform for providing a computational environment comprises a processor (such as CPU) wherein processor time and memory layout such as random access memory (such as RAM) are systems considerations. For example, smaller computer systems offer inexpensive, fast processors and large memory and storage capabilities. In some embodiments, graphics processing units (GPUs) can be used. In some embodiments, hardware platforms for performing computational methods as described herein comprise one or more computer systems with one or more processors. In some embodiments, smaller computer are clustered together to yield a supercomputer network.
[0135] In some embodiments, computational methods as described herein are carried out on a collection of inter- or intra-connected computer systems (such as grid technology) which may run a variety of operating systems in a coordinated manner. For example, the CONDOR framework (University of Wisconsin-Madison) and systems available through United Devices are exemplary of the coordination of multiple stand-alone computer systems for the purpose dealing with large amounts of data. These systems may offer Perl interfaces to submit, monitor and manage large sequence analysis jobs on a cluster in serial or parallel configurations.References
[0136] Each of the following references is incorporated by reference in its entirety.
[0137] Blassel L, Medvedev P, Chikhi R. Mapping-friendly sequence reductions: Going beyond homopolymer compression. iScience. 2022 Oct; 25(11): 105305. doi: 10.1016 / j.isci.2022.105305.
[0138] Chin CS, Behera S, Khalak A, Sedlazeck FJ, Sudmant PH, Wagner J, ZookJM. Multiscale analysis of pangenomes enables improved representation of genomic diversity for repetitive and clinically relevant genes. Nature Methods. 2023 Aug; 20(8): 1213- 1221. doi: 10.1038 / s41592-023-01914-y, number: 8 Publisher: Nature Publishing Group.
[0139] Chin CS, Khalak A, Human Genome Assembly in 100 Minutes. bioRxiv; 2019. doi: 10.1101 / 705616, pages: 705616 Section: New Results.
[0140] Dutta A, Pellow D, Shamir R. Parameterized syncmer schemes improve long-read mapping. PLoS computational biology. 2022 Oct; 18(10):el010638. doi: 10.1371 / j ournal . pcbi .1010638.
[0141] Edgar R. Syncmers are more sensitive than minimizers for selecting conserved k-mers in biological sequences. PeerJ. 2021; 9:el0805. doi: 10.7717 / peerj.10805.
[0142] Ekim B, Sahlin K, Medvedev P, Berger B, Chikhi R. Efficient mapping of accurate long reads in minimizer space with mapquik. Genome Research. 2023 Jul; 33(7): 1188- 1197. doi: 10.1101 / gr.277679.123.
[0143] Lederman R. A random-permutations-based approach to fast read alignment. BMC Bioinformatics. 2013 Apr; 14(Suppl 5):S8. doi: 10.1186 / 1471-2105-14-S5- S8.
[0144] Pandey P, Bender MA, Johnson R, Patro R. A General-Purpose Counting Filter: Making Every Bit Count. In: Proceedings of the 2017 ACM International Conference on Management of Data SIGMOD ‘ 17, New York, NY, USA: Association for Computing Machinery; 2017. p. 775-787. doi.org / 10.1145 / 3035918.3035963, doi: 10.1145 / 3035918.3035963.
[0145] Rautiainen M, Nurk S, Walenz BP, Porubsky D, Rhie A, Phillippy AM. Telomere-to-telomere assembly of diploid chromosomes with Verkko. Nat Biotechnol. 2023; 9:el0805. doi: 10.1038 / s41587-023-01662-6.
[0146] Sahlin K. Effective sequence similarity detection with strobemers. Genome Research. 2021 Nov; 31 (l l):2080-2094. doi: 10.1101 / gr.275648.121.
[0147] Sahlin K, Baudeau T, Cazaux B, Marchet C. A survey of mapping algorithms in the long-reads era. Genome Biology. 2023 Jun; 24(1):133. doi.org / 10.1186 / S13059-023-02972-3, doi: 10.1186 / s 13059-023 -02972-3.EXAMPLES
[0148] Some aspects of the embodiments discussed above are disclosed in further detail in the following examples, which are not in any way intended to limit the scope of the present disclosure. Those in the art will appreciate that many other embodiments also fall within the scope of the disclosure, as it is described herein above and in the claims.Example 1
[0149] The following example describes embodiments of methods for sequence alignment using a reduced alphabet. As mentioned above, the reduced alphabet may be selected based on particular mutations within a nucleotide sequence. For example, an A to G or G to A conversion may be encoded with a single character such as “R”. Similarly, a C to T or T to C conversion may be encoded with a single character, such as “Y” as discussed in more detail below.Mutation biochemistry tailored sequence encodings
[0150] Given an arbitrary starting alphabet S of arbitrary size |S| (typically 4, but possibly 5 or 6 in the case of methylation readouts), a mutation process can be defined as a stochastic matrix of dimension |S |x |S | that gives the probability of any character representing a nucleotide being changed into any other nucleotide character as a result of the mutation process. An alphabet recoding maps the characters of S to an alternative alphabet S', which may be larger or smaller in size. When |S'|x|S| the recoding is referred to as a reduced alphabet encoding. A reduced alphabet encoding that is tailored to a mutation biochemistry, such as converting one nucleotide to another, is an encoding where the values in the diagonal of the stochastic mutation probability matrix of S' are closer to 1. Take, for example, a mutation process that introduces transition mutations for DNA sequences (S = {A, C, G, T}), interconverting A <=> G and C <=> T at a rate of 5%. The stochastic mutation matrix is:0.95 0 0.05 0 i 0 0.95 0 0.05=0.05 0 0.95 0 (1)L 0 0.05 0 0.95-
[0151] An exemplary reduced alphabet encoding would map {A, G] -» R and {C, T] -» Y . The resulting reduced alphabet S' = {R, Y] and has the stochastic mutation matrix:
[0152] In this embodiment, it can be seen that when the sequence is recoded such that A <=> G and C <=> T conversions are each represented by a single character in a reduced alphabet representation, the mutation process does not change the sequence in that there is still a one-to-one mapping of characters to nucleotides. Other mutation processes and reduced alphabet combinations may give values < 1.0 in the diagonal of the stochastic mutation matrix.Windowing guarantee
[0153] Dutta et al. (2022) introduced the concept of a windowing guarantee to improve sensitivity in sequence regions that lack syncmers. A window guarantee is an extension of a parameterized syncmer scheme that ensures every window of length w contains at least one selected k-mer. A typical approach would be to use a minimizer function to select a k-mer if no syncmer exists in a window.Syncstrobes and strobe sets
[0154] At a given k-mer size, the information content of an RY-encoded sequence is half that of a standard DNA sequence. That is, in a simple coding scheme without compression, an RY-encoded k-mer requires k bits of information to represent, while a standard DNA sequence would require 2k bits in a typical encoding. The lower information content means that for any given k, the number of possible k-mers in an RY-encoding is much lower than a standard DNA encoding, and therefore the probability that two random sequences of length k are identical is higher. For example, assuming uniformly random sequences, the probability that two sequences are identical in an RY-encoding is 2'k, while for a standard DNA sequence it is 'k. In the context of match filtration for sequence alignment, the effect when processing RY-encoded sequences is that the size of k must be approximately doubled to obtain equivalent match specificity compared to standard DNA sequences. Typical values of k for processing standard human genome sequences range from 19 to 31 , depending on application context (such as short reads or long reads, mapping or assembly). The problem with simplydoubling the size of k to process RY-encoded sequences is that the probability that a k-mer contains some kind of sequencing error grows. The problem is particularly acute when working with short reads (e.g. 150nt reads), because there are only a small number of k-mers that can serve as match seeds. This problem motivates several innovations to enable high sensitivity and specificity sequence matching in the presence of substitution and indel errors in sequences.
[0155] The concept of syncstrobes is described below. Syncstrobes are a new type of seed for finding matches in sequences in the presence of mutations. FIGS. 4A-4C illustrate some of the ideas and concepts described in the rest of this section.
[0156] FIGS. 4A-4C illustrate the ideas of syncstrobes and strobe sets. Here, syncmers with k = 5, s = 2, are used, with the s-mer located in the center of the k-mer and lexicographic order as the minimizer function; number of components n = 3, target distance between pairs 5 = 18, and window size w = 10. The example is constructed such that it has syncmers in every window and minimizer augmentation was not needed.
[0157] FIG. 4A illustrates syncmer pairs and component selection for a syncstrobe. Each colored bar represents a syncmer. Each syncmer tries to find a pair after a certain distance. Here, syncmers 401 and 405 each make a pair with syncmers 403 and 406 (respectively). Syncmer 402 could not find the syncmer to pair with (because its forward mate syncmer 404 has syncmer 405 as its backward mate). Therefore, there are two syncstrobes here corresponding to two pairs. The pair of syncmers 401 and 403 is selected with syncmer 402 as their internal component since it is the only one appearing in the window after syncmer 401. Similarly, the pair of syncmers 405 and 406 is selected with syncmer 402 as their internal component since it is the minimum of the two syncmers 402 and 403 appearing in the window between them. The construction of this example is such that the left syncmer of the pair is the smaller one; thus the construction of windows is left to right.
[0158] FIG. 4B illustrates strobe sets created by packing the syncstrobes in a set by dropping a component. The two syncstrobes components are then packed to give two strobe sets. Each set has strobemers created by dropping one of the components. For example, Syncstrobe 1 (Strobe Set 1) from above has three strobemers: one created by dropping syncmer 403 and using only syncmers 401 and 402; second created by dropping syncmer 402 and using only syncmers 401 and 403; third created by dropping syncmer 401 and using only syncmers 402 and 403.
[0159] FIG. 4C illustrates how different sequences with transition mutations will produce different syncmers and how using RY-encoded sequence instead would result in the same syncstrobes from a mutated and an unmutated sequence of a locus. Transitions are highlighted here.Generalization of syncmers to transformed sequence encodings
[0160] The concept of a syncmer is extended to apply to transformed sequences. Instead of identifying syncmers on the original sequence, a transformation is first applied to the sequence, for example a reduced alphabet encoding that is tailored to a mutation process (see above). The syncmers in the transformed sequence can then be used in the construction of a sequence index, possibly in conjunction with additional processing (see the notions described below). The sequence index can contain reference to the original untransformed sequence, so that any matches found via the sequence index can be evaluated on the original untransformed sequence.Syncmer-preferred minimizers
[0161] The concept of a minimizer function that prefers syncmers over other k- mers is introduced. In any given window of length w, the minimizer function selects the minimum among the syncmers present in the window, but if no syncmers are present then the minimizer function selects a minimum among the k-mers. This concept is distinct from the windowing guarantee of Dutta et al. (2022) because only a single k-mer can be selected in any given window, whereas application of their windowing guarantee does not select among syncmers in a window when multiple syncmers are present. Additionally, when the minimizer function is strand-symmetric it can be seen that the same set of syncmer-preferred minimizers are produced regardless of whether the forward or reverse complement strand of DNA is processed.Strand-symmetric syncmer-preferred minimizer pairs
[0162] The first concept introduced herein to increase the sensitivity and specificity of RY-encoded sequence match filtration is that of a strand-symmetric syncmer-preferred minimizer pair. The general idea is to find a pair of k-mers that are approximately located some specified distance apart from each other in the sequence. One way this can be done is bydefining a pair of syncmer-preferred minimizers such that the distance between their starting positions is closest to the specified distance.
[0163] Define the parameter 8 > 0 to be the distance between the syncmer- preferred minimizers pairs. Then for position i of each syncmer-preferred minimizer within sequence s, its pair in forward direction will be (if any) the syncmer-preferred minimizer at position j such that j > i + k and 8 — (j — i) is minimum for all such j. Similarly, we can find such pairings in the backward direction for each syncmer-preferred minimizer. The pairings can be made strand- symmetric if only those pairings (i, j) are selected where forward pairing for i is j and backward pairing for j is i (i and j are positions of syncmer-preferred minimizer such that i < j). When implemented in code, it may be necessary to select one of the two kmers to be the “first” in the data structure. To preserve strand symmetry in that case, the minimizer function can be applied again to select the lesser of the two k-mers to be used as the “first” in the data structure.Syncstrobes: syncmer-preferred minimizer strobemers
[0164] When 8 has been chosen to be large enough, it becomes possible to select one or more additional kmers from the sequence between the initial kmer pair. The proposed method selects n — 2 additional kmers in the region between the initial kmer pair, to give n total kmers for indexing. By increasing the value of nk, the specificity of the seed filtration can be increased.
[0165] To define the concept more formally, first denote the distance between the initial k-mer pair as A(i,j) = j — i + k, where i and j are the left- end coordinates in the sequence s of the initial kmer pair and i < j. The proposed method then generates n — 2 additional (mostly) non-overlapping windows of size from which syncmer-preferredminimizers are selected. In order to preserve strand symmetry in the case where division of d(i, j~) by n — 2 does not produce an integer, the windows are allowed to overlap by up to 1 bp on each side. A syncstrobe is then defined by an ordered set of syncmer-preferred minimizers m ... mnwhere m and mnare derived from the initial kmer pair, and when n > 2, m2are derived from the windows between the initial kmer pair. In order to ensure strand symmetry,... mnand the reverse sequence of kmers mn... m are compared to each other and the lesser (under some function, e.g. lexicographic) is chosen for indexing. The keyused for the indexing may be a simple concatenation of the minimizer sequences or may be transformed by a (optionally invertible) hash function.Strobe drop indexing
[0166] In order to further increase sensitivity for sequence match filtration, we introduce the concept of strobe drop indexing. Given a syncstrobe consisting of n kmer components, the idea is to index strobemers that contain subsets of the n kmers. The subsets can be of arbitrary size and one or all or an arbitrary number of subsets can be selected for indexing. A typical size would be n — 1, in which case indexing all n possible subsets of n — 1 kmers yields a guarantee that a match will be generated in the region spanned by the initial kmer pair even when a substitution error exists somewhere in the region. The guarantee is a lower bound guarantee and is not tight. Typically, in the case of RY-encoded syncstrobes, a seed match can be identified in the presence of an arbitrarily large number of transition differences and multiple transversion or indel differences, depending on the exact positions of the non-transition differences.
[0167] Indexing multiple subset strobemers has two potential downsides. One is that the index becomes larger. This downside is to some extent mitigated by the innovations described above which greatly reduce the size of the index.
[0168] A second potential downside is that the approach generates significant extra computation in terms of seed processing. In particular, if all of the subset strobemers produce sequence matches, then in effect the same match could be processed n times. In the case where subset strobemers are generated for all n — 1 subsets, this downside can be mitigated by storing the identity and index of the dropped kmer. In cases where the dropped kmer matches, all n subset strobemers can be expected to trigger matches. In those cases, n — 1 of the matches can be ignored by only processing the match for one particular drop index.Accelerating alignments with syncstrobes
[0169] Concepts that can speed-up the alignments between sequences with matching synctrobes are introduced below.Local repeat marking
[0170] In some low complexity and locally repetitive sequences the identified syncstrobe sequences may not be unique. The local region may produce a large number of syncstrobes. In these cases, the syncstrobe is still useful for seed match filtration, but the coordinates of the seed match may not produce the highest scoring alignment. Rather than indexing a large number of identical syncstrobes that differ only in (nearby) sequence coordinates, only a single representative is indexed and marked as locally repetitive. A downstream alignment process can then consider alternative nearby sequence alignments, rather than taking the seed match coordinates as fixed alignment anchors.Strobe k-mer delta position filtration
[0171] In the process of generating a syncstrobe, the positions of the individual k- mer components of the strobemer can be recorded. The positions can be recorded either as absolute positions in a sequence, or in terms of their distance to the first k-mer in the syncstrobe or distance to the previous k-mer in the syncstrobe. The concept of “delta position filtration” is introduced, where the positions of the k-mer components of the strobemer are considered as part of the match criteria during match filtration. The idea is that when a syncmer sequence match exists, differences in the distances between the k-mer components can imply the existence of alignment gaps. If the k-mer component matches are taken as alignment anchors, then the minimum number and length of those alignment gaps can be inferred from the differences in k-mer position. For some alignment scoring schemes, it becomes possible to compute an upper bound on the alignment score, which could be used as a filtration criterion. Alternatively, filtration could be applied based on the number of alignment gaps, or the minimum implied length of one or more alignment gaps. In practice, filtration on the number of alignment gaps appears to increase specificity and reduce compute time with negligible cost in terms of sensitivity.Drop index match filtration
[0172] In the subsection Strobe drop indexing we described a method that indexes subsets of the n components in a syncstrobe. Drop index filtration is a filtration criterion that checks whether the corresponding components were dropped from two matching syncstrobes.That is, a match would be allowed when both syncstrobes had component i dropped, but not when one syncstrobe had the i-th component dropped and the other had the i-th component dropped, where i #= j.Strobe k-mer alignment anchoring
[0173] In the subsection above entitled “Strobe k-mer delta position filtration” the concept of recording the (relative) positions of the individual k-mer components of a strobemer is described. When available, the positions can be used to constrain further alignment between two sequences that have a syncstrobe match. In cases where the difference in positions is identical among the matched sequences, it is possible to simply score an ungapped alignment, which can be done in a bit-parallel manner. In cases where the relative positions differ, or where an ungapped alignment scores poorly, it is possible to trigger gapped alignment in the region between the matching k-mer components. Any gapped alignment method known in the art can be used, such as Needleman-Wunsh or wavefront alignment, under arbitrary scoring schemes.
[0174] This approach of constraining an alignment using component k-mer matches can greatly reduce the alignment search space. Take as an example two 150 nucleotide (nt) sequences. In the case of a single 21 bp seed match (typical of short read alignment methods such as DRAGEN), in the best case where the match is in the center of the 15 Ont sequences it leaves an area in the dynamic programming matrix of 642+652or about 8300 cells to compute. In the case of a syncstrobe match with four evenly spaced matching k-mer components, where k = 15 and a match span of 100 nt, there would be 252+l 32+l 32+l 42+252or about 1784 cells to compute. In practice, the number of cells actually computed is likely to be much smaller in both cases due to the use of additional heuristics and approaches such as wavefront alignment that further constrain the alignment search space.Homopolymer quantization
[0175] Homopolymer error is a type of sequencing error triggered by various aspects of sequencing biochemistry, such as polymerase strand slippage, 3’ deblocking failure, or uncertainty in signal processing of various sequencing systems. When runs of a single base exist in the template DNA (e.g. AAAAAAAAAA - a 10-base long run of A), homopolymererror can cause the sequence to be misread to have a longer or shorter run of the nucleotide. Previous work introduced homopolymer compression, which is a sequence transformation that reduces homopolymer runs to a single instance of the base, e.g. AAAAAAAAAA becomes A.
[0176] Homopolymer quantization may be introduced to improve the ability of the system to correctly read homopolymers in long homopolymer regions. Homopolymer quantization maps a homopolymer run of length I to one or more nearby values. An example of a homopolymer quantization function is shown in Table 1.Table 1Table 1. An example of homopolymer quantization. Homopolymers in the input sequence are transformed based on their length to have the quantized lengths given in the second column. Any given homopolymer can be quantized to multiple values.
[0177] Unlike homopolymer compression, homopolymer quantization preserves some information about the length of the homopolymer run. Ideally, the homopolymerquantization function will be tailored to the particular error profile of a sequence biochemistry. For example, in some sequencing technologies, homopolymer errors are extremely rare in homopolymer runs of length four or less. In such cases, the optimal quantization function would preserve the lengths of the short homopolymer runs, while quantizing the lengths of longer runs. As noted in Blassel et al. (2022), homopolymer run length information is important for accurate read mapping in some parts of the human genome, such as in centromeric repeat sequences.
[0178] In the context of match seed filtration with syncstrobes, homo polymer quantization would typically be applied prior to application of recoding scheme such as RY- encoding. In the case where the quantization function specifies multiple output values, the number of possible output sequences grows exponentially. However, when used with syncstrobes, it may be sufficient to consider only the combinatorial possibilities that occur within the span of a syncstrobe, i.e. 2k + 8 base pairs.Example 2: RYE aligner algorithm
[0179] In the following example, a specific software implementation of some of the described methods, referred to as the RY-encoded aligner (RYE aligner) is described.Sequence indexing
[0180] The RYE aligner is designed to align a set of short read sequence pairs against a set of 1 or more longer sequences. A typical dataset might have on the order of 300 million short read sequence pairs, with short read sequences ranging from 100-500 nucleotides (nt), and a typical size of 150nt. The long read sequence dataset might typically consist of 30 - 50 million sequences with lengths ranging from 1000 - 30000nt, with a modal length around 5000nt.
[0181] The initial step of the RYE aligner generates indexes of the short and the long sequences. To do so, up to 5 suffixes of the 150nt read are encoded as an RY bitvector, and the corresponding WS bitvector is also generated. The 5 suffixes enable a naive parallel streaming implementation of seed matching / filtration that will be described below. A typical value of 5 is 8. The logical layout of the data structure is shown in FIG. 5. The bitvectors are stored in-memory or in one-or-more files on disk. The assignment of bitvector to a file can becarried out using an arbitrary hash function H (•) on the bitvector. In one embodiment, the hash function produces a uniform distribution of values so that the data can be spread evenly among the files, which can later be processed separately in parallel. In another implementation the hash function groups all values with a common prefix of bits in the same file.
[0182] The long reads are then processed to generate indexes. Long read processing involves indexing subsequences of length equal to the read sequence, e.g. 150nt, spaced every 5 nucleotides apart. Both strands (forward and reverse complement) of the long read sequence are indexed every 5 base pairs. Each entry in the index has a similar layout to that which was generated for the short reads, see FIG. 5. As with the short reads, the entries are hashed into an identical number of files or in-memory buffers using the same hash function H(-).
[0183] Finally, the entries of each index file (or in-memory buffer) are sorted on the length b prefix of the RY bitvector. FIG. 5 illustrates a logical layout of entries in the RYE aligner sequence index. The RY bitvector stores the sequence encoded in a reduced binary RY alphabet, while the WS bitvector stores the same sequence encoded in a reduced binary WS alphabet. When combined, the RY and WS bitvectors encode the original base 4 nucleotide sequence. In the case of indexing short reads L would typically be 151. R would depend on the dataset size but in the case of human genomes a typical value would be in the range of 32 to 40. R must be large enough to support assignment of a unique numeric identifier to each read in the dataset, for example 2Rmust be larger than the number of reads in the dataset. M depends on the read length - in the case of short 151nt reads it may be set to 8 but when indexing long reads it should be set such that 2Mis larger than the longest read. In total for short reads, a typical entry would consume 342 bits or 48 bytes in a word-padded representation on a 64-bit machine. The collection of these records is sorted on the first b bits of the RY bitvector, from which match seeds are generated.Seed matching
[0184] In some embodiments, the seed matching or seed-hits step consists of a linear scan through a pair of sorted files or buffers of RYE sequence index. When identical length b prefixes are found in the files, a match processing step is triggered. Alternative approaches to seed matching are possible, such as processing batches of query sequences against an index of target sequences, or hashing the keys to find matches rather than scanningsorted lists. The match processing consists of computing an ungapped alignment score using a transition-aware scoring scheme. Further processing depends on the alignment score. In the case of a high scoring alignment, the seed hit may be logged for further processing. In the case of a low-scoring ungapped alignment, the seed match could be either (1) subjected to alignment soft-clipping at the 3’ (right) end, (2) subjected to a gapped alignment process anchored from the 3’ end of the length b seed match, and / or (3) discarded. In general, any seed matches that fail to produce a high scoring alignment (whether ungapped, soft-clipped, or partially gapped) will not be passed on to downstream processing steps.Bit-parallel ungapped alignment scoring
[0185] A bit-parallel alignment scoring method is applied to sequences that have produced a seed match. Bit-parallel alignment scoring is possible with any f = |S|x|S| scoring matrix, where each entry provides a score for alignment of two characters in a (possibly recoded) alphabet X. In a preferred implementation we employ a “transition” scoring scheme with identical characters assigned a value of 2 (the diagonal in ), transition differences (A: G, G: A, C: T, T: C) assigned a value of -1, and all other values set to -5. Then the score can be defined as: z = 2(popcount(BRyA ®lvs)j (popcount(-i'BM / s-)}where RYX, RYy, WSX, WSydenote the RY and WS-encoded bitvectors of sequences x and y, where © is the bitwise logical XOR operator, popcount( -) is the population count instruction which counts the number of 1 entries in a bitvector in a single CPU instruction, and z is the resulting alignment score.
[0186] In other words, the approach uses the bitwise logical operators XOR, AND, and NOT to count the number of identical sequence positions, the number of transition differences, and the number of transversion differences. The counting can be done in a small number of CPU instructions without the use of any conditional instructions. The exact degree of parallelism depends on the details of the CPU architecture, but at a minimum the approach could be expected to process 64 sites in parallel on a 64-bit machine, but possibly more whenSIMD operations are available. One implementation focuses on a scoring scheme with three parameters, where it can be seen that the counts for any of the entries of can be produced by combinations of bitwise logical operations and the popcount instruction, supporting bit parallel scoring for any JVC .Paired-end hit selection
[0187] When the query sequences include short read pairs it can be advantageous to combine the separate alignments of each read into a paired-end alignment. One way to do so is to first collect some or all of the seed hits for the query sequence, and then sort the records on target sequence coordinate. Pairs of nearby hits from the same paired-end read can then be found as nearby entries in the sorted list, and when they are in the expected distance and orientation they are combined into paired-end hits. Optionally, the reported alignments can be limited to some number of top scoring alignments.Example 3: Syncstrobe alignment algorithm
[0188] In the following example, an implementation of syncstrobe seeded alignment incorporating the ideas described above is described.Syncstrobe indexing
[0189] The initial step of the method computes an index of syncstrobes in the sequences to be aligned. To do so, the sequence is first recoded in the (optionally reduced) alphabet. Next, the canonical k-mers and 5-mers at each position are computed, where a canonical k-mer (and 5-mer) has the usual definition of being the lesser of the forward and reverse complement k-mer at a position according to some ordering function. In an exemplary implementation, each canonical k-mer is then successively evaluated on the criteria of the chosen parameterized syncmer scheme, namely whether 5-mer minimizers occur in the one- or-more positions specified in the scheme. In practice, the above steps of computing the canonical k-mers, 5-mers, and identifying syncmers can be combined into a single pass through the input sequence.
[0190] The syncmer set is then augmented with minimizers where the gap between the right-end coordinate of a syncmer and the left-end coordinate of the next syncmer is largerthan the selected window size w. This gap is divided into windows of size w (last window having length < w) and minimizer is selected from each window using a hash function. In one incarnation, the hash function used by minimap2 can be used.
[0191] Having identified syncmer-preferred minimizer positions, the next step computes their strand-symmetric pairs and proceeds to identify the internal n — 2 components (syncmer-preferred minimizer) between the pair, n — 1 syncstrobes from each pair are then created by packing (concatenating) the components while dropping / leaving out one from each syncstrobe. The lesser of the sequences (based on any ordering, such as lexicographic) created by concatenating the components in the forward and reverse order of the components is selected as a syncstrobe. Each of the n — 1 syncstrobes such created is then added to the set corresponding to the pair if it does not already exist in the set, or is marked repeated.
[0192] To limit the number of syncstrobes saved per sequence, the sets can be ordered using a function. One solution could be using the lesser of all the syncstrobes in a set, based on some hash function, as its representative. The sets can then be ordered in an increasing / decreasing order of their representatives and the top N could be chosen to be saved. All syncstrobes from such chosen sets are added to the index for the given sequence.
[0193] FIG. 6 illustrates a logical layout of syncstrobes in the syncstrobe aligner sequence index. It could be combined with other read related information like read id and optionally the read sequence (like that in RYE index). The strobemer would typically be under 64 bits. The position bits (P) would depend on the read length - in the case of short 151 nt reads it may be set to 8 but when indexing long reads it should be set such that 2Pis larger than the longest read (usually 16 should suffice). The length bits (L) would depend on different strobemer parameters - a value of 15 should be sufficient for a typical use. The dropped index bits (I) would depend on the number of components n but 8 bits should be more than sufficient. The delta positions bits (D) would be dependent on n but 8 bits should be sufficient for one delta position and delta for the dropped index is not needed to be saved. Assuming n = 4 and k = 15 for strobemers, one syncstrobe might be packed in 64 + 16 + 15 + 1 + 8 + 3 * 8 = 128 bits. If it is saved with only the read-id (which could be assigned 64 bits), one entry in the index can be packed in 192 bits or 24 bytes (no padding needed).Alignment of short reads
[0194] The indexes of mutated and unmutated short reads are then processed in similar fashion of seed-matching step of RYE index (explained in the previous section). The indexes could be sorted and scanned where each syncstrobe match triggers the processing step.
[0195] In some embodiments, the processing step could collect all matching unmutated reads for each mutated read and store only top N matches. One criterion for such a selection could try to maximize the overlap of two reads when the seed serves as an anchor. The rationale behind such criterion is to avoid the processing required in aligning all mutated reads of which only a subset will be used downstream. The downstream application of marking morphomutations (mutations intentionally introduced during library preparation) in the mutated read requires only the “best” aligned unmutated reads covering a mutated read (explained in detail in the subsection “de novo mutation marking”).Ungapped, anchored, and gapped alignments
[0196] For the top N matches chosen for each mutated read, the following method is employed for producing an alignment:• Consider only the highest strobemer length hit at the same position that drops a specific component (e.g. 2nd), skip others (that might be coming from other dropped indices of the same strobe set).• If the syncstrobe is repeated, identify its best anchoring position.• When the syncstrobe / seed is repeated in the mutated read, the position in the original seed might not be the best. The method attempts to find the best alignment position by trying several positions. To avoid trying positions which will definitely not result in better alignments, the first k-mer in the syncstrobe will also be repeated and positions having that k-mer are found.• If the lengths of the syncstrobes are equal and their k-mer delta positions match, compute the ungapped alignment score (in the similar bit-parallel fashion as described in the subsection Bit-parallel ungapped alignment scoring).• If the ungapped alignment is successful (explained in the next subsection), the method may terminate here with the ungapped alignment cigar / score being used for this seed-match; otherwise, move onto next step.• Otherwise, try ungapped alignment within the anchored components.• The components with matched delta positions act as anchors.• The sequences between two anchors, on the left of the first anchor, and on the right of the last anchor are attempted for alignment.* If the sequences being compared are of equal length, attempt an ungapped alignment on them.* For unequal length sequences or unsuccessful ungapped alignment attempts, trigger a gapped alignment (using wavefront method).• Each alignment of a part returns cigar / score for that part which can be concatenated / summed in the final cigar / score for this seed-match.
[0197] FIG. 7 is a conceptual illustration of alignment between two reads with a seed-match. The thick solid black lines 703 represent the reads with components of the seedmatch represented by the rectangles labeled Ml, M3, M4. The seeds here are equal in length, with the same component dropped (M2), and with the same delta positions. The lightly dotted line 701 shows the regions / sequences where the first ungapped alignment attempt will take place. If it fails, then the second ungapped alignment attempt will be in the parts shown with the tightly dotted lines 702. Failure of any part will then be aligned in the gapped way using waveform method.Upper bound alignment score approximation
[0198] An ungapped alignment between two sequences is done and the alignment score is computed in the similar bit-parallel fashion as described in the subsection Bit-parallel ungapped alignment scoring. To decide whether an ungapped alignment is successful, the expected score may be estimated (without computing gapped alignments). The expected minimum score allows for a maximum rate of transition mutations and transversion mutations. In some embodiments, for a given length of sequences being aligned, it can be computed as follows:limit150* lengthExpected score150
[0199] If the ungapped alignment score > expected score, the ungapped alignment can be deemed successful and the score and cigar representing matches at all length positions can be thought of as the output of this attempt. Otherwise, the attempt can be deemed failed.Long read alignment
[0200] The alignments of short unmutated reads onto long mutated reads follows the same process as that of the short read alignments (explained in the previous subsection). Once the alignments of individual reads are available, the usual flow of RYE aligner (explained in the previous section) with paired-end alignments search can be followed.
[0201] The search space of top-scoring alignments can be greatly reduced (length of the long reads essentially means many more seed-matches in comparison to short reads) by any of the usual seed chaining and selection methods. The nature of the syncstrobes allows their seed-hits to be chained in the usual manner.Parameters selection
[0202] To make informed choices of the parameter selection, the relationship of k- mer length may be assessed with uniqueness of genome using the T2T assembly (v0.7) of well characterized human sample E GQQ2 Rautiainen etal. (2023). FIG. 8 shows the plots of number and cumulative number of unique k-mers against their number of occurrences in the genome. As expected, k-mer size 32 in 4-base space has nearly 97% k-mers occurring once or twice. The curve is closely followed by k-mer size of 64 in the 2-base space; providing evidence that the seed size may be based around 64. FIG. 8 illustrates the distribution of unique k-mers in the HG002 T2T assembled genome, k-mer of lengths 32, 64, 96, 120 have been shown in 4- base (prefix orig) and 2-base (prefix mut) space.
[0203] In this example, 60 bits of RY-encoded sequence were used for syncstrobe indexing, k = 15, s = 7, and a syncmer scheme was used that requires a single s-minimizer that is located in the center of the k-mer, 8 = 90 and a> = 20. n = 5 k-mers were used persyncstrobe, and drop indexing of 1 strobe k-mer was applied to produce 60-bit strobemers for indexing.Application to data de novo mutation marking
[0204] Marking morphomutations in short mutated reads is one of the early steps in the previous ICLR pipeline. The previous method uses alignments of the mutated short reads on the reference genome (if available, de novo assembled genome otherwise). For the unmapped reads or parts of the reads, it uses the spaced seed family method (stored in a counting quotient filter such as CQF Pandey et al. (2017)) which does not require a reference. The syncstrobe method enables an all vs all direct mapping of mutated and unmutated short reads which, consequently, allows for marking morphomutations without a reference. The advantage of a reference-free method is quite evident in case of a non-modal organism (with no or low-quality reference genome available) or in the repeat-regions where the mappers usually produce alignments with low quality scores.Process Overview
[0205] The alignments of short mutated reads produced by the Syncstrobe aligner (explained in the subsection Alignment of short reads) are used by this step. The alignments are ordered in decreasing order of their alignment scores. The best scoring alignment serves as the initial template for marking morphomutations and for finding consistent alignments. Intuitively, an alignment is consistent if it covers previously uncovered positions and the morphomutations deduced from it matches the morphomutations deduced so far in the overlapping positions. Alignments are processed in the decreasing order of their scores and the extra positions they cover that are uncovered by the best alignment. If an alignment covers extra positions and produces the same morphomutations in the already covered positions, it is used to mark morphomutations in those extra positions and include the extra positions in covered positions. The alignments are processed until all the alignments have been exhausted or all the positions have been covered. It can be seen that the steps involved in finding consistent reads can be done in a bit-parallel manner.
[0206] Similarly, the marking can be done in a bit-parallel manner as follows:where WSX, WSydenote the WS-encoded bitvectors of sequences x and y, © is the bitwise logical XOR operator, and BMM is the bit vector containing morphomutations (1 for positions with morphomutation, 0 otherwise).
[0207] In one implementation, the gapped alignment step can be moved to this step (from computing alignments step) because the unmutated reads ideally should be ungapped with respect to each other. Thus, only the ungapped alignments are used for identifying consistent reads and compute the (computationally expensive) gapped alignment of only those reads that will be used for morphomutation marking.Use in reference-free rendering
[0208] Similar to morphomutation markings, syncstrobes allow direct (reference- free) alignment of unmarked reads to mutated long sequences which are used in the rendering stage (the final stage of an ICLR pipeline). Prior ICLR methods would use the mappings of unmutated reads on the reference to extract the unmutated short reads to be used for rendering away the morphomutations. With syncstrobes, the output of the syncstrobe aligner for the long reads could be passed to the rendering module without any major changes required.Example 4: Performance evaluations
[0209] In the following example, the performance of the syncstrobes was compared against different kinds of seeds: k-mers, minimizers, ryemers (k-mers in the reduced space), k- mers and minimizers are the usual widely used seeds (see Sahlin et al. (2023) for different seeds and their comparison). The utility of the syncstrobes to the real- world applications was demonstrated by comparing their performance in the mutation marking step of ICLR pipeline against the previous ICLR approach that uses spaced seed method and mappings against a reference genome. These evaluations are done on simulated datasets. The advantage of evaluation on simulated datasets is that the correct matching reads or the morphomutation locations in a mutated read are known, and so the correctness metrics can be evaluated unambiguously. Additionally, to measure and compare the syncstrobe method’s performance on real data (where ground truth of mutations is not known), an enrichment panel dataset wasused (one of the targeted sequencing applications of the ICLR pipeline) and measured the errors (FP+FNs) in small variants called from the final output of the pipeline.Evaluation MethodSimulated data
[0210] Two simulated datasets derived from HG002 T2T assembly were created mimicking the ICLR generation pipeline. The simulator, given a reference genome, creates long reads from it at random positions targeting the given depth and the mean read length (unmutated long reads); introduces mutations in the long reads at the given mutation rate (mutated long reads); creates short reads from each mutated long reads (mutated short reads) following a given sequencing profile. Additionally, the simulator also creates short reads from the given reference genome (unmutated / reference reads). The simulator also outputs the locus from which a read has been generated / extracted along with the information on the mutations introduced.
[0211] Regions with structural variants (SVs) (>50bp insertion or deletion) were selected with respect to the human reference GRCh38 and created the following datasets:• FNTR Based on large tandem repeats (>5kbp in length)• MUTMARK Based on any length tandem repeats and include 5kbp flanks on either side• Tandem repeat regions with SVs were selected because these are relatively difficult regions for aligners to get the mappings correct and have direct impact on variant calling, a crucial genomic application.Evaluation metrics
[0212] Various seeding strategies were evaluated by comparing them on the following criteria:
[0213] Recall of seed-hits. Seed-hits (matches) of short unmutated reads (aka ref reads) against the mutated long reads and the mutated short reads were computed. The hits of the unmutated short reads against long or short mutated reads were converted into the sets contained-in or paired with, respectively. The contained-in for each mutated long read is a setcomprising of the read-ids of the ref reads sharing one or more seeds with it. The paired-with for each mutated short read is the read-id of the ref read sharing one or more seeds with it along with having the maximum overlap of loci. The locus information in the read-names was used to create the true-contained-in and true-paired-with data which was used to compute the recall. Precision or specificity was not computed because the nature of the data (tandem repeats) will essentially result in matches between the reads from different loci. Besides, sensitivity / recall is crucial for this stage of the alignment since filters may be used to avoid processing false matches but cannot retrieve a lost alignment if it is missing from the seed-hits. FIG. 9 conceptually illustrates the notions of contained-in and paired-with (best paired) reads.
[0214] FIG. 9 is an illustration of contained-in and paired-with (best paired) reads. The black line 901 is a reference sequence given to the simulator. The thick line 902 is a random long read generated from the reference sequence. Circles 903, 904, 905 are three mutations introduced by the simulator in the long read. Thin red lines are short mutated reads generated from the long mutated read. Grey lines 906 (above black line 901) are the short (ref) reads generated from the unmutated reference sequence. For a mutated short read, the best paired ref read is the one with the maximum overlap of their originating loci (shown with dashed line 907 for one of the reads). For a mutated long read, the contained reads are all those ref reads that have originated from the same locus.
[0215] De novo morphomutation marking performance Using the alignments from the seed-hits of the unmutated against the mutated short reads, the mutations were marked in the mutated reads. The marked mutations were compared against the true mutations (output by the simulator) to compute precision and recall. Not every base goes under the mutation marking process (different methods including truth generation by a simulator ignore bases based on different rules), recall and Fl scores were computed only for those bases that are covered by truth morphomutations and those generated by a particular method. For the sake of completeness, true recall was included in the evaluations, which considers every base (irrespective of whether it was covered / uncovered by the truth or morphomutation marking process).ResultsPerformance of various syncstrobes against different seeds
[0216] Syncstrobes (with parameters as described in the subsection entitled “Parameters selection”) were compared against three types of seeds, using 60 bits for each seed:• k-mer k-mers, the contiguous bases starting at each position of a given sequence, are the most widely used seeds but are susceptible to mutations. Canonical k-mers (lesser of the k-mer’ s and its reverse complement’s integer representations) were used with k = 30 (each base taking two bits makes it 60 bits).• minimizers Minimizers are the subsampled k-mers with a distance guarantee. A minimum k-mer of length is selected from a window of length w based on a hash function h. In this example, k = 30, w = 19 (minimap2’s default window size), and minimap2’s hash function was used.• ryemer A ryemer is a substring of k contigous bases in the RY-encoded sequence. In this example, k = 60 and the first 8 ryemers from each of the RY-encoded sequence and its reverse complement were used.
[0217] FIG. 10 shows that the syncstrobe seeds outperformed other types of seeds in both kinds of evaluation set-up (especially giving large gain in recall in short mutated vs unmutated reads matches). FIG. 10 is a bar graph illustrating a comparison of recall of paired- with and contained-in reads using matches of various types of seeds. In matching short mutated reads against best paired ref reads, syncstrobes were significantly better than others. In matching long mutated reads against contained ref reads, syncstrobes and ryemers were both significantly better than others with syncstrobes being marginally better amongst the two.Example 5: Performance of various morphomutation marking methods
[0218] In the following example, syncstrobes were compared against the previous method of marking morphomutations in the ICLR pipeline. To make the evaluations comprehensive, only the spaced seeds (CQF) and the mapping only method were included in the comparison. Another simulated dataset (WG T2T) was also included in these evaluations which is based on the whole genome T2T HG002 reference (and not just selected tandemrepeats) to show that the gains that syncstrobe method produces are not limited to small datasets or just the repeat regions.
[0219] FIG. 11 depicts the superior performance of the syncstrobe method over other methods. FIG. 11 is a bar graph illustrating a comparison of the accuracy and recall of morphomutation markings in short mutated reads using different methods. The syncstrobe method was the most accurate and the most sensitive method of all. The syncstrobe method resulted in a roughly 30 to 55 percentage point increase in recall with respect to another reference free method (CQFOnly). It also produced a significant improvement (more than 6 pp increase in recall) in accuracy and recall over the previous method (cqfANDmapping). It was also found to be on par with other mapping methods that relied on higher quality reference sequences. The synsctrobe method was also robust to the data (its size as well as nature) in comparison to other methods.Example 6: Performance of syncstrobe based morphomutation marking on real data
[0220] In the following example, downstream small variants calling was used on the output reads of ICLR pipeline and benchmarked the calls against GIAB’s gold standard truthset. The standard variant calling (included in the ICLR Pipeline itself) produces two sets of calls: one derived from ICLR long reads only and another that has been merged (using machine learning methods) with calls from the short unmutated reads. Both sets were benchmarked for ICLR reads produced by the pipeline using the previous method as well as the syncstrobe method for mutation marking in the first stage. FIG. 12 shows that the (reference-free) syncstrobes used only in the first stage were performing on par with the reference- based mapping method. Further, the syncstrobe method marginally outperformed the previous method, producing fewer FPs and FNs, in case of final merged calls.
[0221] FIG. 12 illustrates a comparison of accuracy and recall of small variants called from the reads using different methods for mutation markings. The reference-free syncstrobe method was found to be on par with the reference-dependent previous method.Example 7: Robustness of syncstrobe based morphomutation marking with respect to mutation rate
[0222] In the following example, the system’s tolerance for increasing mutation rates was tested. Increasing the mutation rate in the ICLR pipeline is expected to help the assembly stage to yield better ICLRs and thus further improve the variant calling. A middlesized simulated dataset was designed using chromosomes 1, 21, and 22 of HG002 (using the same T2T assembly as has been used for other datasets). The rationale of producing this dataset and not using TR based datasets was to compare how syncstrobe based methods would behave against the whole genome when mutation rates were doubled (without having to produce computationally expensive whole genome dataset with double mutation rate). Chr 1 being the largest chromosome captures the complications of the whole genome. To add the complications arising from different chromosomes having similar loci, chromosome 21 and 22 were added whose p-arms are quite short and bear very high similarity.
[0223] The previous ICLR pipeline (baseline) and the syncstrobe based methods used in morphomutation marking (syncstrobe) were compared on the two simulated datasets based on chromosome 1, 21, and 22: CHR012122 (usual 5% mutation rate) and CHR012122_MR2 (double mutation rate such as 10%). FIG. 13 presents the results of the comparison for the mutation-markings where it can be seen that for the usual mutation rate, syncstrobe performed better, following the pattern shown on the other datasets. This improvement in performance over the baseline increased further on the data with double mutation rate. While baseline degraded in recall on doubling the mutation rate, the syncstrobe based method maintains it (precision for both stay roughly the same). This shows the method’s robustness to the mutation rate in the sequence data.
[0224] FIG. 13 is a bar graph illustrating a comparison of accuracy and recall of morphomutation markings in short mutated reads with different mutation rate. Baseline sees a drop in recall of about 9 percentage points (on doubling the mutation rate) while syncstrobe recall is the same as has been on every dataset evaluated.Other Considerations
[0225] Conditional language used herein, such as, among others, “can,” “might,” “may,” “e.g.,” and the like, unless specifically stated otherwise, or otherwise understood withinthe context as used, is generally intended to convey that certain embodiments include, while other embodiments do not include, certain features, elements and / or states. Thus, such conditional language is not generally intended to imply that features, elements and / or states are in any way required for one or more embodiments or that one or more embodiments necessarily include logic for deciding, with or without author input or prompting, whether these features, elements and / or states are included or are to be performed in any particular embodiment. The terms “comprising,” “including,” “having,” “involving,” and the like are synonymous and are used inclusively, in an open-ended fashion, and do not exclude additional elements, features, acts, operations, and so forth. Also, the term “or” is used in its inclusive sense (and not in its exclusive sense) so that when used, for example, to connect a list of elements, the term “or” means one, some, or all of the elements in the list.
[0226] Disjunctive language such as the phrase “at least one of X, Y or Z,” unless specifically stated otherwise, is otherwise understood with the context as used in general to present that an item, term, etc., may be either X, Y or Z, or any combination thereof (such as X, Y and / or Z). Thus, such disjunctive language is not generally intended to, and should not, imply that certain embodiments require at least one of X, at least one of Y or at least one of Z to each be present.
[0227] The terms “about” or “approximate” and the like are synonymous and are used to indicate that the value modified by the term has an understood range associated with it, where the range can be ±20%, ±15%, ±10%, ±5%, or ±1%. The term “substantially” is used to indicate that a result (such as a measurement value) is close to a targeted value, where close can mean, for example, the result is within 80% of the value, within 90% of the value, within 95% of the value, or within 99% of the value.
[0228] Unless otherwise explicitly stated, articles such as “a” or “an” should generally be interpreted to include one or more described items.
[0229] While the above detailed description has shown, described, and pointed out novel features as applied to illustrative embodiments, it will be understood that various omissions, substitutions, and changes in the form and details of the devices or algorithms illustrated can be made without departing from the spirit of the disclosure. As will be recognized, certain embodiments described herein can be embodied within a form that does not provide all of the features and benefits set forth herein, as some features can be used orpracticed separately from others. All changes which come within the meaning and range of equivalency of the claims are to be embraced within their scope.
[0230] It should be appreciated that all combinations of the foregoing concepts (provided such concepts are not mutually inconsistent) are contemplated as being part of the inventive subject matter disclosed herein. In particular, all combinations of claimed subject matter appearing at the end of this disclosure are contemplated as being part of the inventive subject matter disclosed herein.
[0231] The scope of the present disclosure is not intended to be limited by the specific disclosures of examples in this section or elsewhere in this specification, and may be defined by claims as presented in this section or elsewhere in this specification or as presented in the future. The language of the claims is to be interpreted broadly based on the language employed in the claims and not limited to the examples described in the present specification or during the prosecution of the application, which examples are to be construed as nonexclusive.
Claims
WHAT IS CLAIMED IS:
1. A method for aligning sequence reads based on reduced sequence representations, the method comprising: receiving sequence reads generated from nucleic acids; encoding sequence reads with a reduced sequence representation, thereby generating transformed sequence reads; generating seeds from the transformed sequence reads; and matching the seeds, thereby aligning the transformed sequence reads.
2. The method of claim 1, wherein the method comprises encoding the sequence reads using a reduced alphabet.
3. The method of claim 2, wherein the method comprises encoding the sequence reads in RY format.
4. The method of any of claims 1-3, wherein the method comprises reducing a length of homopolymers within the transformed sequence reads.
5. The method of claim 4, wherein the method comprises quantizing the homopolymers by rounding the length of homopolymers within the sequence reads to a predetermined value of a set of two or more predetermined values.
6. The method of any of claims 1-5, generating seeds from the transformed sequence reads comprises selecting strand-symmetric pairs of minimizers from the transformed sequence reads.
7. The method of any of claims 1-6, wherein the method comprises generating a sequence index based on the transformed sequence reads.
8. The method of claim 7, wherein generating a sequence index comprises indexing strobemers by selecting an outer pair of k-mers and selecting one or more additional k-mers between the outer pair of k-mers.
9. The method of claim 8, wherein the outer pair of k-mers, or the one or more additional k-mers, are selected strand-symmetrically.
10. The method of claim 8 or claim 9, wherein selecting k-mers comprises using a minimizer function to select k-mers.
11. The method of claim 10, wherein the minimizer function selects k-mers based on a parameterized syncmer scheme.
12. The method of claim 11, wherein the minimizer function evaluates k-mers based on whether s-mer minimizers occur in one or more pre- determined positions specified in the parameterized syncmer scheme.
13. The method of any of claims 7-12, wherein generating a sequence index comprises generating a set of syncstrobes by successively removing one or more individual k- mer components from one or more strobemers.
14. The method of any of claims 7-13, wherein the method comprises matching seeds from the sequence index, and generating one or more alignments based on seed matches.
15. The method of claim 14, wherein the method comprises using k-mer components of a strobemer seed match as anchors in a gapped sequence alignment.
16. The method of any of claims 14-15, wherein the method comprises scoring an alignment based on whether a distance between two or more k-mer components of a strobemer seed match is consistent in a matched strobemer.
17. The method of any of claims 14-16, wherein the method comprises identifying repetitive strobemers and storing a single representative strobemer to be used as a seed in seed matching.
18. The method of claim 17, wherein the repetitive strobemers are marked as repetitive.
19. The method of any of claims 1-18, wherein the method comprises aligning short sequence reads with a length of 50-500 bp to each other.
20. The method of any of claims 1-19, wherein the method comprises aligning short sequence reads with a length of 50-500 bp to long sequence reads with a length of 50-250,000 bp.
21. The method of any of claims 1-20, wherein aligning the transformed sequence reads comprises aligning without a reference sequence.
22. The method of any of claims 1-21, wherein the sequence reads comprise mutations or sequence errors.
23. The method of claim 22, wherein 4% to 12% of nucleotides comprise a mutation or a sequence error.
24. The method of claim 22 or claim 23, wherein the mutation comprises a transition mutation.
25. The method of any of claims 1-24, wherein the nucleic acids comprise mutated nucleic acids and unmutated nucleic acids, and wherein the sequence reads comprise mutated sequence reads and unmutated sequence reads.
26. The method of any of claims 1-25, wherein the method comprises identifying mutated positions in mutated sequence reads by comparing alignments between a mutated sequence read set and an unmutated read sequence set.
27. The method of claim 26, wherein the method comprises selecting a set of alignments of unmutated reads for each mutated sequence read based on an alignment score.
28. The method of claim 27, wherein the method comprises selecting a set of unmutated sequence reads which are consistent with each other, which have a sequence identity above a threshold to a mutated sequence read covering a corresponding sequence, or which together cover at least a predetermined percentage of the mutated sequence read.
29. The method of any of claims 26 to 28, wherein the method comprises marking mutated positions at sites which differ between the mutated sequence read set and the unmutated sequence read set.
30. The method of any of claims 1-28, wherein the method comprises: generating a sequence index from transformed sequence reads; identifying seed matches from index entries; and generating alignments based on the seed matches.
31. The method of claim 30, wherein generating a sequence index from transformed sequence reads comprises: selecting an outer pair of k-mers and selecting one or more additional k-mers between the outer pair of k-mers using a syncmer- preferred minimizer function; and removing one or more k-mer components from one or more strobemers, thereby generating sets of syncstrobes.
32. The method of any of claims 1-29, wherein the method comprises: generating a sequence index comprising seeds generated from sequence reads encoded with reduced sequence representation; matching seeds from the sequence index and generating alignments;scoring alignments; and combining separate alignments of short sequence read pairs into a paired-end alignment.
33. The method of any of claims 1-29, wherein the method comprises: generating a sequence index of syncstrobes from transformed sequence reads, wherein generating a sequence index comprises: selecting an outer pair of k-mers and selecting one or more additional k- mers between the outer pair of k-mers using a syncmer-preferred minimizer function; removing one or more k-mer components from one or more strobemers, thereby generating sets of syncstrobes; ordering sets of syncstrobes; sorting and scanning indexes of the sequence index, wherein the indexes correspond to mutated sequence reads and unmutated sequence reads, thereby identifying one or more matches; generating one or more ungapped, anchored, or gapped alignments based on the one or more matches; scoring one or more ungapped, anchored, or gapped alignments; and combining separate alignments of short sequence read pairs into a paired-end alignment.
34. An electronic system for aligning sequence reads based on reduced sequence representations, comprising a processor configured to perform a method comprising: receiving sequence reads generated from nucleic acids; encoding sequence reads with a reduced sequence representation, thereby generating transformed sequence reads; generating seeds from the transformed sequence reads; and matching the seeds, thereby aligning the transformed sequence reads.
35. A non-transitory computer-readable medium comprising a plurality of instructions, which when executed by at least one processor, cause the at least one processor to: receive sequence reads generated from nucleic acids; encode sequence reads with a reduced sequence representation, thereby generating transformed sequence reads; generate seeds from the transformed sequence reads; and match the seeds, thereby aligning the transformed sequence reads.
Citation Information
Cited By
Method, system, equipment and medium for assembling third-generation sequencing data based on clustering and graph construction
CN122024841A