Patent Description:
Commercial sequencing systems generally produce either short reads having low error rates (i.e., Illumina sequencers) or long reads having high error rates (i.e., Pacific Biosciences sequencers). Consequently, most sequence alignment tools have been developed and optimized for both these use cases: (<NUM>) identifying short motifs from short reads having low error rates or (<NUM>) identifying long motifs from long reads having high error rates. However, in certain assays it is desirable to be able to identify short motifs from long sequences having high error rates. <CIT> and <CIT>, both from the same applicant, discloses a sequence analysis apparatus for executing a sequence analysis pipeline on genetic sequence data and provide systems, apparatuses, and methods for implementing bioinformatic protocols. However, neither <CIT> and <CIT> explicitly disclose that the selection of the first, initial, seed length in such a method is, or should be, based on an error rate of the sequence reads, nor do these documents disclose of such a method with the iterative repetition of masking and aligning an additional set of sequence reads with each subsequent set having a smaller seed length.

The present invention relates generally to aligning sequences, and more specifically to efficiently identifying short motifs in long sequences with high error rates.

In some embodiments, a method for aligning sequence reads to a reference sequence is provided. The method can include aligning a first set of sequence reads from an entire population of sequence reads to a reference sequence with a Burrows-Wheeler transform using a first seed length, wherein the first seed length is selected based on an error rate of the sequence reads; masking the first set of sequence reads such that the entire population of sequence reads comprises a subset of masked sequence reads and unmasked sequence reads; aligning a second set of sequence reads from the unmasked sequence reads to the reference sequence with the Burrow-Wheeler transform using a second seed length, wherein the second seed length is smaller or shorter than the first seed length, in order to achieve greater sensitivity of the second alignment step; and determining an alignment of the sequence reads to the reference sequence based on the first set of sequence reads and the second set of sequence reads, wherein the method further includes iteratively masking and aligning additional sets of sequence reads with each subsequent set of reads having a greater seed length, and determining the alignment of the sequence reads with the additional sets of sequence reads.

In some embodiments, the first seed length is less than <NUM> bases. In some embodiments, the first seed length is less than <NUM> bases. In some embodiments, the first seed length is <NUM> bases.

In some embodiments, the error rate of the sequence reads is at least <NUM>%. In some embodiments, the error rate of the sequence reads is at least <NUM>%. In some embodiments, the error rate of the sequence reads is at least <NUM>%.

In some embodiments, the sequence reads are sequenced from a plurality of concatamers, wherein each concatamer is formed of oligonucleotide sequences that have been joined together, wherein the oligonucleotide sequences correspond to a plurality of loci from a set of chromosomes. In some embodiments, the set of chromosomes comprises chromosome <NUM>, <NUM>, <NUM>, X, and Y. In some embodiments, the set of chromosomes is selected from the group consisting of chromosome <NUM>, <NUM>, <NUM>, X, and Y. In some embodiments, the method further includes calculating a frequency each loci is found in the sequence reads.

Also disclosed but not part of the invention, a method for aligning sequence reads to a reference sequence is provided. The method can include aligning a first set of sequence reads from an entire population of sequence reads to a reference sequence with a Burrows-Wheeler transform using a first set of sensitivity parameters, wherein the first set of sensitivity parameters is selected based on an error rate of the sequence reads; masking the first set of sequence reads such that the entire population of sequence reads comprises a subset of masked sequence reads and unmasked sequence reads; aligning a second set of sequence reads from the unmasked sequence reads to the reference sequence with the Burrow-Wheeler transform using a second set of sensitivity parameters, wherein the second set of sensitivity parameters results in a higher sensitivity than the first set of sensitivity parameters; determining an alignment of the sequence reads to the reference sequence based on the first set of sequence reads and the second set of sequence reads.

In some embodiments, the method further includes iteratively masking and aligning additional sets of sequence reads with each subsequent set of reads having a set of sensitivity parameters that results in higher sensitivity, and determining the alignment of the sequence reads with the additional sets of sequence reads.

In some embodiments, the sensitivity parameters are selected from the group consisting of seed generation, chaining and filtering, and thresholding.

In some embodiments, a computer product includes a computer readable medium storing a plurality of instructions for controlling a computer system to perform any of the methods above.

In some embodiments, a system includes the computer product above and one or more processors for executing instructions stored on the computer readable medium.

In some embodiments, a system includes means for performing any of the above methods.

In some embodiments, a system includes one or more processors configured to perform any of the above methods.

In some embodiments, a system includes modules that respectively perform the steps of any of the above methods.

A nanopore based sequencing system can generate vast amounts of long read sequencing data, depending on the number of nanopore based sensors that are fabricated in the sensing chip. In some embodiments, each sensing chip can have millions of cells, each with a nanopore sensor. One advantage of using a nanopore sequencer is the ability to generate long reads. However, the current raw read accuracy of nanopore sequencers is generally lower (i.e., between <NUM> to <NUM>%) than established short read technologies (i.e., accuracy greater than <NUM>%).

To take advantage of the long read capability of a nanopore sequencer, relatively short target sequences can be concatemerized into a longer linear sequence that can be efficiently sequenced in the nanopore sequencer. The concatemers can be sequenced in parallel by the nanopore sequencer to generate a consensus sequence with the desired degree of accuracy (i.e., greater than <NUM>%, <NUM>%, or <NUM>%, for example). In order to generate the consensus sequence, the sequence fragments need to be aligned.

Although the alignment methods and systems are described in the context of a nanopore sequencer, other sequences generated by other types of sequencing devices can also be aligned using the alignment methods and systems described herein. For example, Non-limiting examples of sequence assays that are suitable for use with the methods disclosed herein include nanopore sequencing (<CIT>, <CIT>, <CIT>, <CIT> and <CIT>), Sanger sequencing, capillary array sequencing, thermal cycle sequencing (<NPL>)), solid-phase sequencing (<NPL>)), sequencing with mass spectrometry such as matrix-assisted laser desorption/ionization time-of-flight mass spectrometry (MALDI-TOF/MS;<NPL>)), sequencing by hybridization (<NPL>), and NGS methods, including but not limited to sequencing by synthesis (e.g., HiSeq™, MiSeq™, or Genome Analyzer, each available from Illumina), sequencing by ligation (e.g., SOLiD™, Life Technologies), ion semiconductor sequencing (e.g., Ion Torrent™, Life Technologies), and SMRT® sequencing (e.g., Pacific Biosciences).

Commercially available sequencing technologies include: sequencing-by-hybridization platforms from Affymetrix Inc. (Sunnyvale, Calif. ), sequencing-by-synthesis platforms from Illumina/Solexa (San Diego, Calif. ) and Helicos Biosciences (Cambridge, Mass. ), sequencing-by-ligation platform from Applied Biosystems (Foster City, Calif). Other sequencing technologies include, but are not limited to, the Ion Torrent technology (ThermoFisher Scientific), and nanopore sequencing (Genia Technology from Roche Sequencing Solutions, Santa Clara, Cal. ); and Oxford Nanopore Technologies (Oxford, United Kingdom).

While alignment analysis is not a new topic, certain sequencing applications present a challenge. For example, the noninvasive prenatal testing (NIPT) assay presents a case that is not well developed in the world of bioinformatics. Identification of short motifs, for example in one embodiment NIPT indices of <NUM> base-pair, from long and noisy read sequences (in one embodiment concatemer nanopore reads with <NUM>-<NUM> base-pairs and <NUM>-<NUM>% of error) is a new consideration to algorithm design because of the needs to satisfy both compute time / resource constraints as well as to meet strict sensitivity and specificity requirements of the assay. Most existing alignment tools and methods are designed for short motifs sequenced at low error rates (e.g., Illumina sequencers) and for long motifs with high error rate (e.g., Pacific Bioscience sequencers): in fact, the NIPT assay's use-case of short motifs and high error rate is exactly the weakness of these approaches, which suffer in terms of sensitivity, compute resource or both.

As a concrete consideration, genomic aligners are designed for cases wherein only a few very long reference sequences (e.g., chromosomes of millions base-pairs) are used as the alignment target and each read sequence is expected to be aligned uniquely to one single target locus. In contrast, in one embodiment the NIPT assay probes biological samples with artificial short sequence motifs (NIPT Indices), cleaves them with a restriction enzyme (i.e. Pstl) and concatenates them into a long template (i.e., greater than <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, or <NUM> bases), containing nongenomic short oligos, and sequences the template as long, high error reads (i.e. at least <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, or <NUM> percent). The NIPT assay can have over <NUM> very short reference sequences (NIPT indices of <NUM> base-pairs) and each concatemer read sequence aligns as a set of non-overlapping segments to multiple targets at up to <NUM>% error rate, as shown in <FIG>. The goal of secondary analysis is to identify those indices in concatemer reads, i.e. to find frequency of each database index in given patient sample. Index frequencies are then sent to a tertiary analysis module (i.e., "Forte") which builds statistical model from chromosomal frequencies and establishes screening test outcome (probability), recommendation for clinician, etc..

Similar reasons, as well as low noise resilience, are applicable against analytical tools covering other classes of alignment problems, making them unfit or non-optimal for the type of the alignment problems found in certain types of sequencing based assays. For example, short motifs, and well-defined ground truth (e.g., set of <NUM>+ target NIPT indices without possible allelic polymorphisms and/or other inconsistencies typical for genomic analysis) would make this type of problem seemingly fit for an alignment-free approach, e.g. using sequence k-mer seeding/hashing. However, the presence of high noise in reads would require very short k-mer length to maintain high sensitivity. For example, internal data showed in practice that changing k-mer length from <NUM> base-pairs to <NUM> base-pairs could lead to loss of sensitivity of <NUM>% or more, depending on the read error-rate. Short(er) k-mers also create additional computational burden (typically complexity is O(n<NUM>), where n is number of k-mers, inversely proportional to k-mer length).

We conceived a novel alignment method which leverages multi-stage secondary analysis, with each stage progressively reducing the amount of data to be analyzed in the next stage(s), but increasing exhaustiveness of the search on the remaining data received from previous stage(s). This way, less noisy alignments can be quickly identified from the initially large data-pools in early stage(s), while very noisy alignments can be identified equally fast from smaller data-pools in latter stage(s) of computation, thus maintaining target sensitivity while reducing overall compute times.

We found that the alignment method requires a non-trivial computational implementation, because we showed that a naive, non-optimized implementation of the alignment method would lead to additional auxiliary data/compute cost (which negates the alignment method's benefit) and also potential incompatibility with downstream processing.

To address all these, we implemented our alignment method as a highly optimized and concise extension to the industry-standard BWA-MEM, a well-established open-source genomic alignment tool under the Apache license. We had to change and tune its core algorithms, data structures and context-dependent optimizations to fit the unique nature of alignment problem we had at hand for the NIPT assay. We also had to change its overall infrastructure to accommodate our multi-stage analysis functionality. More specifically, the following parts of BWA-MEM were changed or tuned: core algorithms (seed generation, chaining and filtering, thresholding, etc.), data structures (reads, references, alignment info, SAM records, reduced I/O), memory usage and access patterns (I/O buffering, arrays creation and propagation), and the overall infrastructure and logic of the flow of execution (to accommodate multi-stage analysis functionality we added a new iterative alignment procedure to BWA, added iteration-specific parameter creation, checkpointing, etc.). A flow chart of the modified and optimized workflow of BWA-MEM is presented in <FIG>.

For example, in some embodiments the seed parameter is reduced from a <NUM> base or base pair kmer (i.e., BWA-MEM default) to a <NUM> base or base pair kmer. The size of the reduced kmer is based on the error rate. If the error rate of the sequence reads is decreased, then the seed length can be increased to reduce the computation time. An increased first pass seed size, as compared to the second pass seed size, allows the first pass to quickly identify oligo sequences that match an NIPT index (but with lower sensitivity), which allows these sequences to be removed from the unaligned sequence pool, which can be processed using higher sensitivity parameters in the next pass, such as a smaller seed size (seed generation, chaining and filtering, thresholding, etc.). The process of seeding involves converting the kmer to an integer value (i.e., each base type has a particular value and the kmer is the sum of the base values). These integer values can be much more easily and quickly compared with the corresponding kmer integer values from the reference NIPT index than trying to match the kmer strings themselves.

Along with seed generation, chaining can be tuned to improve performance. Chaining is the process of finding groups of seeds, called chains, that are overlapping seeds or seeds that are co-linear and close to each other. The chains can be identified by various grouping algorithms, such as the pigeon hole method or greedy clustering algorithm.

This project was exciting because we have developed a novel use for arguably the most popular bioinformatics tool of today. Importance of this project lies in the fact that it has enabled nanopore sequencing based NIPT assay secondary analysis to be rapidly completed (i.e., less than <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, or <NUM> hours). It also enabled analysis on virtually any kind of commodity hardware, a stark contrast to initial project estimates for very expensive high-end computing hardware.

Work on this project also lays ground for similar algorithmic approaches in other areas (i.e. sequencing based assays) where long reads carrying short and noisy motifs could be leveraged, but lack highly optimized tools and algorithms, for example, inference of splicevariant events in RNA-seq isoforms, meta-transcriptomics, etc..

One of the impacts of the optimized alignment method is in direct cost reduction for the final product. Initial estimates before optimization for hardware configuration required for secondary analysis for the NIPT assay were predicting compute workstation with <NUM>-<NUM> GB of RAM and <NUM>-<NUM> CPU cores. The optimized alignment method reduced compute time by 2X (reduction in CPU cores) and expected RAM memory usage by <NUM>%. It is hard to estimate exact cost-reduction impact, but it is fair to say it will be in tens of thousands of dollars per each installed configuration, potentially saving up to <NUM>% of computing hardware cost.

A second impact is reduction in compute time. By reducing compute time, it was also possible to fit per-run cycle analysis into a <NUM>-hour time allotment, so that secondary analysis on previous run cycle can be performed while instrument is sequencing next run cycle, effectively making secondary run analysis results ready within one hour after sequencing experiment completion. This significant product requirement was not feasible and at risk with the previous, unoptimized approach.

The alignment method described herein is particularly suited for identifying short motifs from long, high error reads. In some embodiments, a short motif is less than about <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, or <NUM> bases or base pairs. In some embodiments, a high error read is a raw read with an error rate of greater than about <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, or <NUM>%. In some embodiments, long reads have at least <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, or <NUM> bases or base pairs. In some embodiments, the algorithm can be applied to shorter length high error reads, depending on the output of the sequencer. For example, shorter reads can be less than <NUM>, <NUM>, <NUM>, <NUM>, <NUM>, or <NUM> bases.

The alignment methods and algorithms described herein can be implemented on a computer system. For example, <FIG> is a block diagram illustrating one embodiment of a computer system <NUM> configured to implement one or more aspects of the present invention. As shown, computer system <NUM> includes, without limitation, a central processing unit (CPU) <NUM> and a system memory <NUM> coupled to a parallel processing subsystem <NUM> via a memory bridge <NUM> and a communication path <NUM>. Memory bridge <NUM> is further coupled to an I/O (input/output) bridge <NUM>via a communication path <NUM>, and I/O bridge <NUM> is, in turn, coupled to a switch <NUM>.

In operation, I/O bridge <NUM> is configured to receive user input information from input devices <NUM> (e.g., a keyboard, a mouse, a video/image capture device, etc.) and forward the input information to CPU <NUM> for processing via communication path <NUM> and memory bridge <NUM>. In some embodiments, the input information is a live feed from a camera/image capture device or video data stored on a digital storage media on which object detection operations execute. Switch <NUM> is configured to provide connections between I/O bridge <NUM> and other components of the computer system <NUM>, such as a network adapter <NUM> and various add-in cards <NUM> and <NUM>.

As also shown, I/O bridge <NUM> is coupled to a system disk <NUM> that may be configured to store content and applications and data for use by CPU <NUM>and parallel processing subsystem <NUM>. As a general matter, system disk <NUM> provides non-volatile storage for applications and data and may include fixed or removable hard disk drives, flash memory devices, and CD-ROM (compact disc read-only-memory), DVD-ROM (digital versatile disc-ROM), Blu-ray, HD-DVD (high definition DVD), or other magnetic, optical, or solid state storage devices. Finally, although not explicitly shown, other components, such as universal serial bus or other port connections, compact disc drives, digital versatile disc drives, film recording devices, and the like, may be connected to I/O bridge <NUM> as well.

In various embodiments, memory bridge <NUM> may be a Northbridge chip, and I/O bridge <NUM> may be a Southbrige chip. In addition, communication paths <NUM> and <NUM>, as well as other communication paths within computer system <NUM>, may be implemented using any technically suitable protocols, including, without limitation, AGP (Accelerated Graphics Port), HyperTransport, or any other bus or point-to-point communication protocol known in the art.

In some embodiments, parallel processing subsystem <NUM> comprises a graphics subsystem that delivers pixels to a display device <NUM> that may be any conventional cathode ray tube, liquid crystal display, light-emitting diode display, or the like. In such embodiments, the parallel processing subsystem <NUM> incorporates circuitry optimized for graphics and video processing, including, for example, video output circuitry. Such circuitry may be incorporated across one or more parallel processing units (PPUs) included within parallel processing subsystem <NUM>. In other embodiments, the parallel processing subsystem <NUM> incorporates circuitry optimized for general purpose and/or compute processing. Again, such circuitry may be incorporated across one or more PPUs included within parallel processing subsystem <NUM> that are configured to perform such general purpose and/or compute operations. In yet other embodiments, the one or more PPUs included within parallel processing subsystem <NUM> may be configured to perform graphics processing, general purpose processing, and compute processing operations. System memory <NUM> includes at least one device driver <NUM> configured to manage the processing operations of the one or more PPUs within parallel processing subsystem <NUM>. The system memory <NUM> also includes a software application <NUM> that executes on the CPU <NUM> and may issue commands that control the operation of the PPUs.

In various embodiments, parallel processing subsystem <NUM> may be integrated with one or more other the other elements of <FIG> to form a single system. For example, parallel processing subsystem <NUM> may be integrated with CPU <NUM> and other connection circuitry on a single chip to form a system on chip (SoC).

It will be appreciated that the system shown herein is illustrative and that variations and modifications are possible. The connection topology, including the number and arrangement of bridges, the number of CPUs <NUM>, and the number of parallel processing subsystems <NUM>, may be modified as desired. For example, in some embodiments, system memory <NUM> could be connected to CPU <NUM> directly rather than through memory bridge <NUM>, and other devices would communicate with system memory <NUM> via memory bridge <NUM> and CPU <NUM>. In other alternative topologies, parallel processing subsystem <NUM> may be connected to I/O bridge <NUM> or directly to CPU <NUM>, rather than to memory bridge <NUM>. In still other embodiments, I/O bridge <NUM> and memory bridge <NUM> may be integrated into a single chip instead of existing as one or more discrete devices. Lastly, in certain embodiments, one or more components shown in <FIG> may not be present. For example, switch <NUM> could be eliminated, and network adapter <NUM> and add-in cards <NUM>, <NUM> would connect directly to I/O bridge <NUM>.

<FIG> compares key cost/performance metrics between traditional BWA (blue) and v0. <NUM> Multistage Aligner (orange) executing on vhmem hardware of RSS-SC1 cluster to process San Jose's test case 20180503_uzuki_000001_WSK18R04C8_L03. Cost (first five bars) are drastically reduced, while integrity (last four bars) are equivalent.

Claim 1:
A method for aligning sequence reads to a reference sequence, the method comprising:
aligning a first set of sequence reads from an entire population of sequence reads to a reference sequence with a Burrows-Wheeler transform using a first seed length, wherein the first seed length is selected based on an error rate of the sequence reads;
masking the first set of sequence reads such that the entire population of sequence reads comprises a subset of masked sequence reads and unmasked sequence reads;
aligning a second set of sequence reads from the unmasked sequence reads to the reference sequence with the Burrow-Wheeler transform using a second seed length, wherein the second seed length is smaller than the first seed length;
determining an alignment of the sequence reads to the reference sequence based on the first set of sequence reads and the second set of sequence reads; and wherein the method further comprises:
iteratively masking and aligning additional sets of sequence reads with each subsequent set of reads having a smaller seed length, and determining the alignment of the sequence reads with the additional sets of sequence reads.