Bioinformatics systems, apparatuses, and methods for performing secondary and / or tertiary processing

Abstract

A system, method and apparatus for executing a bioinformatics analysis on genetic sequence data is provided. Particularly, a genomics analysis platform for executing a sequence analysis pipeline is provided. The genomics analysis platform includes one or more of a first integrated circuit, where each first integrated circuit forms a central processing unit (CPU) that is responsive to one or more software algorithms that are configured to instruct the CPU to perform a first set of genomic processing steps of the sequence analysis pipeline. Additionally, a second integrated circuit is also provided, where each second integrated circuit forming a field programmable gate array (FPGA), the FPGA being configured by firmware to arrange a set of hardwired digital logic circuits that are interconnected by a plurality of physical interconnects to perform a second set of genomic processing steps of the sequence analysis pipeline, the set of hardwired digital logic circuits of each FPGA being arranged as a set of processing engines to perform the second set of genomic processing steps. A shared memory is also provided.

Description

(19) *EP004682891A2* (11) EP 4 682 891 A2 (12) EUROPEAN PATENT APPLICATION (43) Date of publication: 21.01.2026 Bulletin 2026 / 04 (21) Application number: 25210051.6 (22) Date of filing: 27.10.2017 (51) International Patent Classification (IPC): G16B 50 / 50 (2019.01) (52) Cooperative Patent Classification (CPC): G16B 30 / 10; G16B 20 / 20; G16B 20 / 40; G16B 30 / 20; G16B 50 / 30; G16B 50 / 40; G16B 50 / 50; G06F 21 / 76 (84) Designated Contracting States: AL AT BE BG CH CY CZ DE DK EE ES FI FR GB GR HR HU IE IS IT LI LT LU LV MC MK MT NL NO PL PT RO RS SE SI SK SM TR (30) Priority: 28.10.2016 US 201662414637 P (62) Document number(s) of the earlier application(s) in accordance with Art. 76 EPC: 17798358.2 / 3 532 967 (71) Applicant: ILLUMINA, INC. San Diego, CA 92122 (US) (72) Inventors: • HAHM, Mark, David San Diego, CA 92122 (US) • DE BEER, Jacobus San Diego, CA 92122 (US) • JAIN, Varun San Diego, CA 92122 (US) • MEHIO, Rami San Diego, CA 92122 (US) • OJARD, Eric San Diego, CA 92122 (US) • RUEHLE, Michael San Diego, CA 92122 (US) • PTASHEK, Amnon San Diego, CA 92122 (US) • CATREUX, Severine San Diego, CA 92122 (US) • VISVANATH, Arun San Diego, CA 92122 (US) (74) Representative: Lowden, Samuel Robert James Marks & Clerk LLP 2nd Floor Wytham Court 11 West Way Oxford OX2 0JB (GB) Remarks: •The complete document including Reference Table(s) and the Sequence Listing(s) can be downloaded from the EPO website •This application was filed on 21‑10‑2025 as a divisional application to the application mentioned under INID code 62. •Claims filed after the date of filing of the application / after the date of receipt of the divisional application (Rule 68(4) EPC). (54) BIOINFORMATICS SYSTEMS, APPARATUSES, AND METHODS FOR PERFORMING SECONDARY AND / OR TERTIARY PROCESSING (57) Asystem,methodandapparatus for executing a bioinformatics analysis on genetic sequence data is pro- vided. Particularly, a genomics analysis platform for ex- ecuting a sequence analysis pipeline is provided. The genomicsanalysis platform includes oneormoreof a first integratedcircuit,whereeach first integratedcircuit forms a central processing unit (CPU) that is responsive to one or more software algorithms that are configured to in- struct the CPU to perform a first set of genomic proces- sing steps of the sequence analysis pipeline. Addition- ally, a second integrated circuit is also provided, where each second integrated circuit forming a field program- mable gate array (FPGA), the FPGA being configured by firmware to arrangeaset of hardwireddigital logic circuits that are interconnected by a plurality of physical inter- connects to perform a second set of genomic processing steps of the sequence analysis pipeline, the set of hard- wired digital logic circuits of each FPGA being arranged as a set of processing engines to perform the second set of genomic processing steps. A shared memory is also provided. EP 4 68 2 89 1 A 2 Processed by Luminess, 75001 PARIS (FR) (Cont. next page) 2 EP 4 682 891 A2 Description Related Application

[0001] This application claims priority to U.S. Provisional Patent Application Serial No. 62 / 414,637, filed onOctober 28, 2016, the contents of which are hereby fully incorporated by reference. Field of the Disclosure

[0002] The subject matter described herein relates to bioinformatics, and more particularly to systems, apparatuses, and methods for implementing bioinformatic protocols, such as performing one or more functions for analyzing genomic data on an integrated circuit, such as on a hardware processing platform. Background to the Disclosure

[0003] As described in detail herein, some major computational challenges for high-throughput DNA sequencing analysis is to address the explosive growth in available genomic data, the need for increased accuracy and sensitivity when gathering that data, and the need for fast, efficient, and accurate computational tools when performing analysis on a wide range of sequencing data sets derived from such genomic data.

[0004] Keeping pace with such increased sequencing throughput generated by Next Gen Sequencers has typically beenmanifestedasmultithreadedsoftware tools that havebeenexecutedonever greater numbersof faster processors in computer clusters with expensive high availability storage that requires substantial power and significant ITsupport costs. Importantly, future increases in sequencing throughput rates will translate into accelerating real dollar costs for these secondary processing solutions.

[0005] The devices, systems, andmethods of their use described herein are provided, at least in part, so as to address these and other such challenges. Summary of the Disclosure

[0006] The present disclosure is directed to devices, systems, andmethods for employing the same in the performance of oneormoregenomicsand / or bioinformaticsprotocolsondatagenerated throughaprimaryprocessingprocedure, such as on genetic sequence data. For instance, in various aspects, the devices, systems, and methods herein provided are configured for performing secondary and / or tertiary analysis protocols on genetic data, such as data generated by the sequencing of RNA and / or DNA, e.g., by a Next Gen Sequencer ("NGS"). In particular embodiments, one or more secondary and / or tertiary processingpipelines for processinggenetic sequencedata is provided.Specifically, oneormore tertiary processing pipelines for processing genetic sequence data is provided, such as where the pipelines, and / or individual elements thereof, deliver superior sensitivity and improved accuracy on awider rangeof sequence derived data than is currently available in the art.

[0007] For example, provided herein is a system, such as for executing one or more of a sequence and / or genomic analysis pipeline ongenetic sequencedata and / or other data derived therefrom. In various embodiments, the systemmay include one or more of an electronic data source that provides digital signals representing a plurality of reads of genetic and / or genomic data, such aswhere each of the plurality of reads of genomic data include a sequence of nucleotides. The systemmay further include amemory, e.g., a DRAM, or a cache, such as for storing one or more of the sequenced reads, one or a plurality of genetic reference sequences, and one or more indices of the one or more genetic reference sequences. The systemmay additionally include one ormore integrated circuits, such as aFPGA,ASIC, or sASIC, and / or a CPU and / or a GPU and / or Quantum Processing Units (QPUs), which integrated circuit, e.g., with respect to the FPGA, ASIC, or sASIC may be formed of a set of hardwired digital logic circuits that are interconnected by a plurality of physical electrical interconnects. The systemmay additionally include a quantumcomputing processing unit, for use in implement- ing one or more of the methods disclosed herein.

[0008] In various embodiments, one ormore of the plurality of electrical interconnectsmay include an input to the one or more integrated circuits that may be connected or connectable, e.g., directly, via a suitable wired connection, or indirectly such as via a wireless network connection (for instance, a cloud or hybrid cloud), with the electronic data source. Regardless of a connection with the sequencer, an integrated circuit of the disclosuremay be configured for receiving the plurality of reads of genomic data, e.g., directly from the sequencer or from an associated memory. The reads may be digitally encoded in a standardFASTQorBCL file format. Accordingly, the systemmay include an integrated circuit having oneormoreelectrical interconnects thatmaybeaphysical interconnect that includesamemory interface soas toallow the integrated circuit to access the memory.

[0009] Particularly, the hardwired digital logic circuit of the integrated circuit may be arranged as a set of processing 3 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 engines, such as where each processing engine may be formed of a subset of the hardwired digital logic circuits so as to perform one or more steps in the sequence, genomic, and / or tertiary analysis pipeline, as described herein below, on the plurality of reads of genetic data as well as on other data derived therefrom. For instance, each subset of the hardwired digital logic circuits may be in a wired configuration to perform the one or more steps in the analysis pipeline. Additionally, where the integrated circuit is anFPGA, such steps in the sequenceand / or further analysis processmay involve thepartial reconfiguration of the FPGA during the analysis process.

[0010] Particularly, the set of processing engines may include a mapping module, e.g., in a wired configuration, to access, according to at least some of the sequence of nucleotides in a read of the plurality of reads, the index of the one or more genetic reference sequences, from the memory via the memory interface, so as to map the read to one or more segments of the one ormore genetic reference sequences based on the index. Additionally, the set of processing engines may include an alignmentmodule in the wired configuration to access the one ormore genetic reference sequences from thememory via thememory interface to align the read, e.g., themapped read, to one ormore positions in the one ormore segments of theoneormoregenetic referencesequences, e.g., as received from themappingmoduleand / or stored in the memory.

[0011] Further, the set of processing enginesmay include a sortingmodule so as to sort each aligned read according to theoneormorepositions in theoneormoregenetic referencesequences.Furthermore, thesetof processingenginesmay include a variant call module, such as for processing the mapped, aligned, and / or sorted reads, such as with respect to a reference genome, to thereby produce an HMM readout and / or variant call file for use with and / or detailing the variations between the sequenced genetic data and the reference genomic reference data. In various instances, one or more of the plurality of physical electrical interconnectsmay includeanoutput from the integrated circuit for communicating result data from the mapping module and / or the alignment and / or sorting and / or variant call modules.

[0012] Particularly, with respect to the mapping module, in various embodiments, a system for executing a mapping analysis pipeline on a plurality of reads of genetic data using an index of genetic reference data is provided. In various instances, the genetic sequence, e.g., read, and / or the genetic reference data may be represented by a sequence of nucleotides, whichmay be stored in amemory of the system. Themappingmodule may be included within the integrated circuit and may be formed of a set of pre-configured and / or hardwired digital logic circuits that are interconnected by a plurality of physical electrical interconnects, which physical electrical interconnects may include a memory interface for allowing the integrated circuit to access the memory. In more particular embodiments, the hardwired digital logic circuits may be arranged as a set of processing engines, such as where each processing engine is formed of a subset of the hardwired digital logic circuits to perform one or more steps in the sequence analysis pipeline on the plurality of reads of genomic data.

[0013] For instance, in one embodiment, the set of processing engines may include a mapping module in a hardwired configuration, where the mapping module, and / or one or more processing engines thereof is configured for receiving a read of genomic data, such as via one ormore of a plurality of physical electrical interconnects, and for extracting a portion of the read in such a manner as to generate a seed therefrom. In such an instance, the read may be represented by a sequence of nucleotides, and the seed may represent a subset of the sequence of nucleotides represented by the read. Themappingmodulemay includeor beconnectable toamemory that includesoneormoreof the reads, oneormoreof the seeds of the reads, at least a portion of one ormore of the reference genomes, and / or one ormore indexes, such an index built from the one or more reference genomes. In certain instances, a processing engine of the mapping module employ the seed and the index to calculate an address within the index based on the seed.

[0014] Once an address has been calculated or otherwise derived and / or stored, such as in an onboard or offboard memory, the address may be accessed in the index in the memory so as to receive a record from the address, such as a record representing position information in the genetic reference sequence. This position informationmay then be used to determine one or more matching positions from the read to the genetic reference sequence based on the record. Then at least one of the matching positions may be output to the memory via the memory interface.

[0015] In another embodiment, a set of the processing engines may include an alignment module, such as in a pre- configured and / or hardwired configuration. In this instance, one or more of the processing engines may be configured to receive one or more of the mapped positions for the read data via one or more of the plurality of physical electrical interconnects. Then thememory (internal or external)may be accessed for eachmapped position to retrieve a segment of the reference sequence / genome corresponding to the mapped position. An alignment of the read to each retrieved reference segment may be calculated along with a score for the alignment. Once calculated, at least one best-scoring alignment of the read may be selected and output. In various instances, the alignment module may also implement a dynamic programming algorithm when calculating the alignment, such as one or more of a Smith-Waterman algorithm, e.g., with linear or affine gap scoring, a gapped alignment algorithm, and / or a gapless alignment algorithm. In particular instances, the calculating of the alignment may include first performing a gapless alignment to each reference segment, and based on the gapless alignment results, selecting reference segments with which to further perform gapped alignments.

[0016] In various embodiments, a variant call module may be provided for performing improved variant call functions 4 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 that when implemented in one or both of software and / or hardware configurations generate superior processing speed, better processed result accuracy, and enhanced overall efficiency than the methods, devices, and systems currently known in the art. Specifically, in one aspect, improvedmethods for performing variant call operations in software and / or in hardware, suchas forperformingoneormoreHMMoperationsongenetic sequencedata, areprovided. Inanotheraspect, novel devices including an integrated circuit for performing such improved variant call operations, where at least a portion of the variant call operation is implemented in hardware, are provided.

[0017] Accordingly, in various instances, the methods disclosed herein may include mapping, by a first subset of hardwired and / or quantum digital logic circuits, a plurality of reads to one or more segments of one or more genetic referencesequences.Additionally, themethodsmay includeaccessing, by the integratedand / or quantumcircuits, e.g., by one or more of the plurality of physical electrical interconnects, from the memory or a cache associated therewith, one or more of themapped reads and / or one ormore of the genetic reference sequences; and aligning, by a second subset of the hardwired and / or quantum digital logic circuits, the plurality of mapped reads to the one or more segments of the one or more genetic reference sequences.

[0018] In various embodiments, the method may additionally include accessing, by the integrated and / or quantum circuit, e.g., by one or more of the plurality of physical electrical interconnects from a memory or a cache associated therewith, the aligned plurality of reads. In such an instance the method may include sorting, by a third subset of the hardwired and / or quantumdigital logic circuits, the aligned plurality of reads according to their positions in the one ormore genetic reference sequences. In certain instances, the method may further include outputting, such as by one or more of the plurality of physical electrical interconnects of the integrated and / or quantum circuit, result data from the mapping and / or thealigningand / or thesorting, suchaswhere the result data includespositionsof themappedand / oralignedand / or sorted plurality of reads.

[0019] In some instances, the method may additionally include using the obtained result data, such as by a further subset of the hardwired and / or quantum digital logic circuits, for the purpose of determining how the mapped, aligned, and / or sorted data, derived from the subject’s sequenced genetic sample, differs from a reference sequence, so as to produce a variant call file delineating the genetic differences between the two samples. Accordingly, in various embodi- ments, the method may further include accessing, by the integrated and / or quantum circuit, e.g., by one or more of the plurality of physical electrical interconnects from a memory or a cache associated therewith, the mapped and / or aligned and / or sorted plurality of reads. In such an instance the method may include performing a variant call function, e.g., an HMMor paired HMMoperation, on the accessed reads, by a third or fourth subset of the hardwired and / or quantum digital logic circuits, soas toproduceavariant call filedetailinghow themapped,aligned,and / or sorted readsvary from thatof one or more reference, e.g., haplotype, sequences.

[0020] Accordingly, in accordance with particular aspects of the disclosure, presented herein is a compact hardware, e.g., chip based, or quantum accelerated platform for performing secondary and / or tertiary analyses on genetic and / or genomic sequencing data. Particularly, a platform or pipeline of hardwired and / or quantum digital logic circuits that have specifically been designed for performing secondary and / or tertiary genetic analysis, such as on sequenced genetic data, or genomic data derived therefrom, is provided. Particularly, a set of hardwired digital and / or quantum logic circuits, which may be arranged as a set of processing engines,may be provided, such aswhere the processing enginesmay be present in a preconfigured and / or hardwired and / or quantum configuration on a processing platformof the disclosure, andmay be specifically designed for performing secondary mapping and / or aligning and / or variant call operations related to genetic analysis on DNA and / or RNA data, and / or may be specifically designed for performing other tertiary processing on the results data.

[0021] In particular instances, the present devices, systems, andmethods of employing the same in the performance of oneormoregenomicsand / or bioinformatics secondary and / or tertiary processing protocols, havebeenoptimized soas to deliver an improvement in processing speed that is orders of magnitude faster than standard secondary processing pipelines that are implemented in software. Additionally, the pipelines and / or components thereof as set forth herein provide better sensitivity and accuracy on a wide range of sequence derived data sets for the purposes of genomics and bioinformatics processing. In various instances, one or more of these operations may be performed on by an integrated circuit that is part of or configured as a general purpose central processing unit and / or a graphics processing unit and / or a quantum processing unit.

[0022] For example, genomics and bioinformatics are fields concerned with the application of information technology and computer science to the field of genetics and / or molecular biology. In particular, bioinformatics techniques can be applied to process and analyze various genetic and / or genomic data, such as from an individual, so as to determine qualitativeandquantitative informationabout that data that can thenbeusedbyvariouspractitioners in thedevelopment of prophylactic, therapeutic, and / or diagnostic methods for preventing, treating, ameliorating, and / or at least identifying diseased states and / or their potential, and thus, improving the safety, quality, and effectiveness of health care on an individualized level. Hence, because of their focus on advancing personalized healthcare, genomics and bioinformatics fieldspromote individualizedhealthcare that isproactive, insteadof reactive,and thisgives thesubject inneedof treatment theopportunity tobecomemore involved in their ownwellness.Anadvantageof employing thegenetics, genomics, and / or 5 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 bioinformatics technologies disclosed herein is that the qualitative and / or quantitative analyses of molecular biological, e.g., genetic, datacanbeperformedonabroader rangeof samplesetsat amuchhigher rateof speedandoften timesmore accurately, thus expediting the emergence of a personalized healthcare system. Particularly, in various embodiments, the genomics and / or bioinformatics related tasks may form a genomics pipeline that includes one or more of a micro-array analysis pipeline, agenome,e.g.,wholegenomeanalysis pipeline, genotypinganalysis pipeline, exomeanalysis pipeline, epigenomeanalysis pipeline,metagenomeanalysis pipeline,microbiomeanalysis pipeline, genotyping analysis pipeline, including joint genotyping, variantsanalysis pipelines, includingstructural variants, somatic variants, andGATK,aswell as RNA sequencing and other genetic analyses pipelines.

[0023] Accordingly, to make use of these advantages there exists enhanced and more accurate software implementa- tions for performing one or a series of such bioinformatics based analytical techniques, such as for deployment by a general purpose CPU and / or GPU and / or may be implemented in one or more quantum circuits of a quantum processing platform. However, common characteristics of traditionally configured software based bioinformatics methods and systems is that they are labor intensive, take a long time to execute on such general purpose processors, and are prone to errors. Therefore, bioinformatics systems as implemented herein that could perform these algorithms, such as implemented in software by a CPU and / or GPU of quantum processing unit in a less labor and / or processing intensive manner with a greater percentage accuracy would be useful.

[0024] Such implementations have been developed and are presented herein, such as where the genomics and / or bioinformatics analyses are performed by optimized software run on a CPU and / or GPU and / or quantum computer in a system that makes use of the genetic sequence data derived by the processing units and / or integrated circuits of the disclosure. Further, it is to benoted that the cost of analyzing, storing, and sharing this rawdigital data has far outpaced the cost of producing it. Accordingly, also presented herein are "just in time" storage and / or retrievalmethods that optimize the storageof suchdata inamanner that substitutes thespeedof regenerating thedata inexchange for thecost of storingsuch data collectively.Hence, thedatageneration, analysis, and "just in time"or "JIT" storagemethodspresentedherein solvea key bottleneck that is a long felt but unmet obstacle standing between the ever-growing raw data generation and storage and the real medical insight being sought from it.

[0025] Presented herein, therefore, are systems, apparatuses, and methods for implementing genomics and / or bioinformatic protocols or portions thereof, such as for performing one or more functions for analyzing genomic data, for instance, on one or both of an integrated circuit, such as on a hardware processing platform, and a general purpose processor, such as for performing one ormore bioanalytic operations in software and / or on firmware. For example, as set forth herein below, in various implementations, an integrated circuit and / or quantum circuit is provided so as to accelerate one or more processes in a primary, secondary, and / or tertiary processing platform. In various instances, the integrated circuit may be employed in performing genetic analytic related tasks, such as mapping, aligning, variant calling, compressing, decompressing, and the like, in an accelerated manner, and as such the integrated circuit may include a hardware accelerated configuration. Additionally, in various instances, an integrated and / or quantum circuit may be providedsuchaswhere thecircuit ispart of aprocessingunit that is configured forperformingoneormoregenomicsand / or bioinformatics protocols on the generated mapped and / or aligned and / or variant called data.

[0026] Particularly, in a first embodiment, a first integrated circuitmay be formedof anFPGA,ASIC, and / or sASIC that is coupled to or otherwise attached to themotherboard and configured, or in the case of an FPGAmay be programmable by firmware to be configured, as a set of hardwired digital logic circuits that are adapted to perform at least a first set of sequence analysis functions in a genomics analysis pipeline, such as where the integrated circuit is configured as described herein above to include one ormore digital logic circuits that are arranged as a set of processing engines, which are adapted to perform one ormore steps in amapping, aligning, and / or variant calling operation on the genetic data so as to produce sequence analysis results data. The first integrated circuit may further include an output, e.g., formed of a plurality of physical electrical interconnects, such as for communicating the result data from the mapping and / or the alignment and / or other procedures to the memory.

[0027] Additionally, a second integrated and / or quantumcircuitmaybe included, coupled to or otherwiseattached to the motherboard, and in communication with the memory via a communications interface. The second integrated and / or quantum circuit may be formed as a central processing unit (CPU) or graphics processing unit (GPU) or quantum processing unit (QPU) that is configured for receiving themapped and / or aligned and / or variant called sequence analysis result dataandmaybeadapted tobe responsive tooneormoresoftwarealgorithms that are configured to instruct theCPU orGPU to performone ormore genomics and / or bioinformatics functions of the genomic analysis pipeline on themapped, aligned, and / or variant called sequenceanalysis result data.Specifically, thegenomicsand / or bioinformatics related tasks may formagenomics analysis pipeline that includes one ormore of amicro-array analysis, a genomepipeline, e.g., whole genome analysis pipeline, genotyping analysis pipeline, exome analysis pipeline, epigenome analysis pipeline, meta- genome analysis pipeline, microbiome analysis pipeline, genotyping analyses pipelines, including joint genotyping, variants analyses pipelines, including structural variants, somatic variants, and GATK, as well as RNA sequencing analysis pipeline and other genetic analyses pipelines.

[0028] For instance, in one embodiment, the CPU and / or GPUand / or QPU of the second integrated circuit may include 6 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 software that is configured for arranging the genome analysis pipeline for executing a whole genome analysis pipeline, such as a whole genome analysis pipeline that includes one or more of genome-wide variation analysis, whole-exome DNA analysis, whole transcriptome RNA analysis, gene function analysis, protein function analysis, protein binding analysis, quantitative gene analysis, and / or a gene assembly analysis. In certain instances, the whole genome analysis pipelinemaybeperformed for thepurposesofoneormoreofancestryanalysis, personalmedical historyanalysis, disease diagnostics, drug discovery, and / or protein profiling. In a particular instance, the whole genome analysis pipeline is performed for the purposes of oncology analysis. In various instances, the results of this datamay bemade available, e.g. globally, throughout the system.

[0029] In various instances, the CPU and / or GPU and / or a quantum processing unit (QPU) of the second integrated and / or quantum circuit may include software that is configured for arranging the genome analysis pipeline for executing a genotyping analysis, such as a genotyping analysis including joint genotyping. For instance, the joint genotyping analysis may be performed using a Bayesian probability calculation, such as a Bayesian probability calculation that results in an absolute probability that a given determined genotype is a true genotype. In other instances, the software may be configured for performing ametagenomeanalysis so as to producemetagenome result data thatmay in turn be employed in the performance of a microbiome analysis.

[0030] In certain instances, the first and / or second integratedcircuit and / or thememorymaybehousedonanexpansion card, such as a peripheral component interconnect (PCI) card. For instance, in various embodiments, one or more of the integrated circuits may be one or more chips coupled to a PCIe card or otherwise associated with the motherboard. In various instances, the integrated and / or quantum circuit(s) and / or chip(s) may be a component within a sequencer or computer, or server, such as part of a server farm. In particular embodiments, the integrated and / or quantum circuit(s) and / or expansion card(s) and / or computer(s) and / or server(s) may be accessible via the internet, e.g., cloud.

[0031] Further, in some instances, the memory may be a volatile random access memory (RAM), e.g., a direct access memory (DRAM). Particularly, in various embodiments, the memory may include at least two memories, such as a first memory that is anHMEM,e.g., for storing the referencehaplotype sequencedata, andasecondmemory that is anRMEM, e.g., for storing the read of genomic sequence data. In particular instances, each of the twomemoriesmay include awrite port and / or a read port, such aswhere thewrite port and the read port each accessing a separate clock. Additionally, each of the two memories may include a flip-flop configuration for storing a multiplicity of genetic sequence and / or processing result data.

[0032] The details of one or more variations of the subject matter described herein are set forth in the accompanying drawingsand thedescriptionbelow.Other featuresandadvantagesof thesubjectmatter describedhereinwill beapparent from thedescription anddrawings, and from the claims.While certain features of the currently disclosed subjectmatter are described for illustrativepurposes in relation toanenterprise resourcesoftwaresystemorotherbusinesssoftwaresolution or architecture, it should be readily understood that such features are not intended to be limiting. The claims that follow this disclosure are intended to define the scope of the protected subject matter. Brief Description of the Figures

[0033] The accompanying drawings, which are incorporated in and constitute a part of this specification, show certain aspects of the subject matter disclosed herein and, together with the description, help explain some of the principles associated with the disclosed implementations. FIG.1Adepictsasequencingplatformwithaplurality of genetic samples thereon, apluralityof exemplary tilesarealso depicted, as well as a three-dimensional representation of the sequenced reads. FIG. 1B depicts a representation of a flow cell with the various lanes represented. FIG. 1C depicts a lower corner of the flow cell platform of FIG. 1B, showing a constellation of sequenced reads. FIG. 1D depicts a virtual array of the results of the sequencing performed on the reads of FIGS. 1 and 2, where the reads are set forth in an output column by column order. FIG. 1E depicts the method by which the transposition of the outcome reads from column by column order to row by row read order may be implemented. FIG. 1F depicts the transposition of the outcome reads from column by column order, to row by row read order. FIG. 1G depicts the system components for performing the transposition. 7 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 FIG 1H depicts the transposition order. FIG. 1I depicts the architecture for electronically transposing the sequenced data. FIG. 2 depicts an HMM3-state basedmodel illustrating the transition probabilities of going from one state to another. FIG. 3A depicts a high-level view of an integrated circuit of the disclosure including a HMM interface structure. FIG. 3B depicts the integrated circuit of FIG. 3A, showing an HMM cluster features in greater detail. FIG. 4 depicts an overview of HMM related data flow throughout the system including both software and hardware interactions. FIG. 5 depicts exemplary HMM cluster collar connections. FIG. 6 depicts a high-level view of the major functional blocks within an exemplary HMM hardware accelerator. FIG. 7 depicts an exemplary HMM matrix structure and hardware processing flow. FIG. 8depicts anenlarged viewof a portionof FIG. 2 showing thedata flowanddependencies betweennearby cells in the HMM M, I, and D state computations within the matrix. FIG. 9 depicts exemplary computations useful for M, I, D state updates. FIG. 10 depicts M, I, and D state update circuits, including the effects of simplifying assumptions of FIG. 9 related to transition probabilities and the effect of sharing some M, I, D adder resources with the final sum operations. FIG. 11 depicts Log domain M, I, D state calculation details. FIG. 12A depicts an HMM state transition diagram showing the relation between GOP, GCP and transition probabilities. FIG. 12B depicts a particular embodiment of an exemplary HMM state transition diagram showing the relation between GOP, GCP and transition probabilities. FIG. 12C depicts a pileup of a region in the genome evidencing short tandem repeats (STR). FIG. 12D depicts an area under the curve graph expressing indels within a given region. FIG. 13 depicts an HMM Transprobs and Priors generation circuit to support the general state transition diagram of FIG. 12. FIG. 14 depicts a simplified HMM state transition diagram showing the relation between GOP, GCP and transition probabilities. FIG. 15 depicts a HMM Transprobs and Priors generation circuit to support the simplified state transition. FIG. 16 depicts an exemplary theoretical HMM matrix and illustrates how such an HMM matrix may be traversed. FIG. 17A presents a method for performing a multi-region joint detection pre-processing procedure. FIG. 17B presents an exemplary method for computing a connectionmatrix such as in the pre-processing procedure of FIG. 17A. FIG. 18A depicts an exemplary event between two homologous sequenced regions in a pileup of reads. FIG. 18B depicts the constructed reads of FIG. 18A, demarcating nucleotide difference between the two sequences. 8 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 FIG. 18C depicts various bubbles of a De Brujin graph that may be used in performing an accelerated variant call operation. FIG. 18D depicts a representation of a pruning the tree function as described herein. FIG. 18E depicts one of the bubbles of FIG. 18C. FIG. 19 is a graphical representation of the exemplary pileup pursuant to the connection matrix of FIG. 17. FIG. 20 is a processing matrix for performing the pre-processing procedure of FIGS. 17A and B. FIG. 21 is an example of a bubble formation in a De Brujin graph in accordance with the methods of FIG. 20. FIG. 22 is an example of a variant pathway through an exemplary De Brujin graph. FIG. 23 is a graphical representation of an exemplary sorting function. FIG. 24 is another example of a processing matrix for a pruned multi-region joint detection procedure. FIG. 25 illustrates a joint pileup of paired reads for two regions. FIG. 26 sets forth a probability table in accordance with the disclosed herein. FIG. 27 is a further example of a processing matrix for a multi-region joint detection procedure. FIG. 28 represents a selection of candidate solutions for the joint pile up of FIG. 25. FIG. 29 representsa further selectionof candidate solutions for thepile upofFIG.28, after a pruning functionhasbeen performed. FIG. 30 represents the final candidates of FIG. 28, and their associated probabilities, after the performanceof aMRJD function. FIG. 31 illustrates the ROC curves for MRJD and a conventional detector. FIG. 32 illustrates the same results of FIG. 31 displayed as a function of the sequence similarity of the references. FIG. 33A depicts an exemplary architecture illustrating a loose coupling between a CPU and an FPGA of the disclosure. FIG.33Bdepictsanexemplaryarchitecture illustratinga tight couplingbetweenaCPUandanFPGAof thedisclosure. FIG. 34A depicts a direct coupling of a CPU and a FPGA of the disclosure. FIG. 34B depicts an alternative embodiment of the direct coupling of a CPU and a FPGA of FIG. 34A. FIG. 35 depicts an embodiment of a package of a combinedCPUand FPGA,where the two devices share a common memory and / or cache. FIG. 36 illustrates a core of CPUs sharing one ormorememories and / or caches, wherein theCPUs are configured for communicating with one or more FPGAs that may also include a shared or common memory or caches. FIG. 37 illustrates an exemplary method of data transfer throughout the system. FIG. 38 depicts the embodiment of FIG. 36 in greater detail. FIG. 39 depicts an exemplary method for the processing of one or more jobs of a system of the disclosure. 9 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 FIG. 40Adepicts a block diagram for a genomic infrastructure for onsite and / or cloud basedgenomics processing and analysis. FIG. 40B depicts a block diagram of a cloud-based genomics processing platform for performing the BioITanalysis disclosed herein. FIG. 40C depicts a block diagram for an exemplary genomic processing and analysis pipeline. FIG. 40D depicts a block diagram for an exemplary genomic processing and analysis pipeline. FIG. 41A depicts a block diagram of a local and / or cloud based computing function of FIG. 40A for a genomic infrastructure for onsite and / or cloud based genomics processing and analysis. FIG. 41B depicts the block diagram of FIG. 41A illustrating greater detail regarding the computing function for a genomic infrastructure for onsite and / or cloud based genomics processing and analysis. FIG. 41Cdepicts the block diagramof FIG. 40 illustrating greater detail regarding the 3rd‑Party analytics function for a genomic infrastructure for onsite and / or cloud based genomics processing and analysis. FIG. 42A depicts a block diagram illustrating a hybrid cloud configuration. FIG. 42B depicts the block diagram of FIG. 42A in greater detail, illustrating a hybrid cloud configuration. FIG. 42C depicts the block diagram of FIG. 42A in greater detail, illustrating a hybrid cloud configuration. FIG. 43A depicts a block diagram illustrating a primary, secondary, and / or tertiary analysis pipeline as presented herein. FIG. 43Bprovides an exemplary tertiary processing epigenetics analysis for execution by themethods anddevices of the system herein. FIG. 43Cprovidesanexemplary tertiary processingmethylation analysis for execution by themethodsanddevices of the system herein. FIG. 43D provides an exemplary tertiary processing structural variants analysis for execution by the methods and devices of the system herein. FIG. 43E provides an exemplary tertiary cohort processing analysis for execution by the methods and devices of the system herein. FIG. 43F provides an exemplary joint genotyping tertiary processing analysis for execution by the methods and devices of the system herein. FIG. 44 depicts a flow diagram for an analysis pipeline of the disclosure. FIG. 45 is a block diagram of a hardware processor architecture in accordance with an implementation of the disclosure. FIG. 46 is a block diagram of a hardware processor architecture in accordance with another implementation. FIG. 47 is a block diagram of a hardware processor architecture in accordance with yet another implementation. FIG. 48 illustrates a genetic sequence analysis pipeline. FIG. 49 illustrates processing steps using a genetic sequence analysis hardware platform. FIG. 50A illustrates an apparatus in accordance with an implementation of the disclosure. 10 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 FIG. 50B illustrates another apparatus in accordance with an alternative implementation of the disclosure. FIG. 51 illustrates a genomics processing system in accordance with an implementation. Detailed Description of the Disclosure

[0034] As summarized above, the present disclosure is directed to devices, systems, and methods for employing the same in the performance of one or more genomics and / or bioinformatics protocols, such as a mapping, aligning, sorting, and / or variant call protocol ondatagenerated throughaprimaryprocessingprocedure, suchasongenetic sequencedata. For instance, in various aspects, the devices, systems, and methods herein provided are configured for performing secondary analysis protocols on genetic data, such as data generated by the sequencing of RNA and / or DNA, e.g., by a Next Gen Sequencer ("NGS"). In particular embodiments, one or more secondary processing pipelines for processing genetic sequence data is provided, such aswhere the pipelines, and / or individual elements thereof, may be implemented in software, hardware, or a combination thereof in a distributed and / or an optimized fashion so as to deliver superior sensitivity and improved accuracy on a wider range of sequence derived data than is currently available in the art. Additionally, as summarizedabove, the present disclosure is directed to devices, systems, andmethods for employing the same in the performance of one ormore genomics and / or bioinformatics tertiary protocols, such as amicro-array analysis protocol, a genome, e.g., whole genome analysis protocol, genotyping analysis protocol, exome analysis protocol, epigenome analysis protocol, metagenome analysis protocol, microbiome analysis protocol, genotyping analysis pro- tocol, including joint genotyping, variants analysis protocols, including structural variants, somatic variants, andGATK, as well asRNAsequencing protocols andother genetic analyses protocols suchasonmapped, aligned, and / or other genetic sequence data, such as employing one or more variant call files.

[0035] Accordingly, provided herein are software and / or hardware e.g., chip based, accelerated platform analysis technologies for performing secondary and / or tertiary analysis of DNA / RNA sequencing data. More particularly, a platform, or pipeline, of processing engines, such as in a software implemented and / or hardwired configuration, which havespecifically beendesigned for performing secondarygenetic analysis, e.g.,mapping, aligning, sorting, and / or variant calling; and / or may be specifically designed for performing tertiary genetic analysis, such as a micro-array analysis, a genome,e.g.,wholegenomeanalysis, genotypinganalysis, exomeanalysis, epigenomeanalysis,metagenomeanalysis, microbiome analysis, genotyping analysis, including joint genotyping analysis, variants analysis, including structural variants analysis, somatic variants analysis, and GATK analysis, as well as RNA sequencing analysis and other genetic analysis, such aswith respect to genetic based sequencing data, whichmay have been generated in an optimized format that delivers an improvement in processing speed that ismagnitudes faster than standard pipelines that are implemented in known software alone. Additionally, the pipelines presented herein provide better sensitivity and accuracy on a wide range of sequence derived data sets, such as on nucleic acid or protein derived sequences.

[0036] As indicated above, in various instances, it is a goal of bioinformatics processing to determine individual genomes and / or protein sequences of people, which determinations may be used in gene discovery protocols as well as for prophylaxis and / or therapeutic regimes tobetter enhance the livelihoodof eachparticular personandhumankindas awhole. Further, knowledge of an individual’s genomeand / or protein compellationmay beused such as in drug discovery and / or FDA trials to better predict with particularity which, if any, drugs will be likely to work on an individual and / or which would be likely to have deleterious side effects, such as by analyzing the individual’s genome and / or a protein profile derived therefrom and comparing the same with predicted biological response from such drug administration.

[0037] Such bioinformatics processing usually involves three well defined, but typically separate phases of information processing. The first phase, termed primary processing, involves DNA / RNA sequencing, where a subject’s DNA and / or RNA is obtained and subjected to various processes whereby the subject’s genetic code is converted to a machine- readable digital code, e.g., a FASTQ file. The second phase, termed secondary processing, involves using the subject’s generated digital genetic code for the determination of the individual’s genetic makeup, e.g., determining the individual’s genomic nucleotide sequence. And the third phase, termed tertiary processing, involves performingoneormoreanalyses on the subject’s genetic makeup so as to determine therapeutically useful information therefrom.

[0038] Accordingly, once a subject’s genetic code is sequenced, such as by a NextGen sequencer, so as to produce a machine readable digital representation of the subject’s genetic code, e.g., in a FASTQ and / or BCL file format, it may be useful to further process the digitally encoded genetic sequence data obtained from the sequencer and / or sequencing protocol, such as by subjecting digitally represented data to secondary processing. This secondary processing, for instance, can be used to map and / or align and / or otherwise assemble an entire genomic and / or protein profile of an individual, such as where the individual’s entire genetic makeup is determined, for instance, where each and every nucleotide of each and every chromosome is determined in sequential order such that the composition of the individual’s entire genome has been identified. In such processing, the genome of the individual may be assembled such as by comparison to a reference genome, such as a reference standard, e.g., one or more genomes obtained from the human genome project or the like, so as to determine how the individual’s geneticmakeup differs from that of the referent(s). This 11 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 process is commonly known as variant calling. As the difference between the DNA of any one person to another is 1 in 1,000 base pairs, such a variant calling process can be very labor and time intensive, requiringmany steps thatmay need to be performed one after the other and / or simultaneously, such as in a pipeline, so to analyze the subject’s genomic data and determine how that genetic sequence differs from a given reference.

[0039] Inperformingasecondaryanalysispipeline, suchas forgeneratingavariant call file for agivenquerysequenceof an individual subject; a genetic sample, e.g., DNA,RNA, protein sample, or the likemaybeobtained, form the subject. The subject’s DNA / RNA may then be sequenced, e.g., by a NextGen Sequencer (NGS) and / or a sequencer-on-a-chip technology, e.g., in a primary processing step, so as to produce a multiplicity of read sequence segments ("reads") covering all or a portion of the individual’s genome, such as in an oversampledmanner. The end product generated by the sequencing device may be a collection of short sequences, e.g., reads, that represent small segments of the subject’s genome,e.g., short genetic sequences representing the individual’sentire genome.As indicated, typically, the information represented by these readsmay be an image file or in a digital format, such as in FASTQ, BCL, or other similar file format.

[0040] Particularly, in a typical secondary processing protocol, a subject’s geneticmakeup is assembled by comparison toa referencegenome. This comparison involves the reconstruction of the individual’s genome frommillionsuponmillions of short read sequences and / or the comparison of the whole of the individual’s DNA to an exemplary DNA sequence model. In a typical secondary processing protocol an image, FASTQ, and / or BCL file is received from the sequencer containing the raw sequenced read data. In order to compare the subject’s genome to that of the standard reference genome, it needs to be determinedwhere each of these readsmap to the reference genome, such as how each is aligned with respect to one another, and / or how each read can also be sorted by chromosome order so as to determine at what position and inwhich chromosomeeach readbelongs.Oneormoreof these functionsmay takeplaceprior to performinga variant call function on the entire full-length sequence, e.g., once assembled. Specifically, once it is determined where in thegenomeeach readbelongs, the full lengthgenetic sequencemaybedetermined, and then thedifferencesbetween the subject’s genetic code and that of the referent can be assessed.

[0041] For instance, reference based assembly in a typical secondary processing assembly protocol involves the comparison of sequenced genomic DNA / RNA of a subject to that of one or more standards, e.g., known reference sequences. Various mapping, aligning, sorting, and / or variant calling algorithms have been developed to help expedite these processes. These algorithms, therefore, may include some variation of one or more of: mapping, aligning, and / or sorting the millions of reads received from the image, FASTQ, and / or BCL file communicated by the sequencer, to determine where on each chromosome each particular read is located. It is noted that these processes may be implemented in software or hardware, such as by the methods and / or devices described in U.S. Patent Nos. 9,014,989 and 9,235,680 both assigned to Edico Genome Corporation and incorporated by reference herein in their entireties.Oftena common feature behind the functioningof these various algorithmsand / or hardware implementations is their use of an index and / or an array to expedite their processing function.

[0042] For example, with respect to mapping, a large quantity, e.g., all, of the sequenced reads may be processed to determine thepossible locations in the referencegenome towhich those readscouldpossibly align.Onemethodology that canbeused for this purpose is todoadirect comparisonof the read to the referencegenomesoas to findall thepositionsof matching. Another methodology is to employ a prefix or suffix array, or to build out a prefix or suffix tree, for the purpose of mapping the reads to various positions in the reference genome. A typical algorithmuseful in performing such a function is a Burrows-Wheeler transform, which is used to map a selection of reads to a reference using a compression formula that compresses repeating sequences of data.

[0043] Additionally, analigning functionmaybeperformed todetermineout of all thepossible locationsagiven readmay map to on agenome, such as in those instanceswhere a readmaymap tomultiple positions in the genome,which is in fact the location fromwhich it actuallywasderived, suchasbybeing sequenced therefromby theoriginal sequencingprotocol. This function may be performed on a number of the reads, e.g., mapped reads, of the genome and a string of ordered nucleotide bases representing a portion or the entire genetic sequence of the subject’sDNA / RNAmay be obtained. Along with theorderedgenetic sequenceascoremaybegiven for eachnucleotide in agivenposition, representing the likelihood that for any given nucleotide position, the nucleotide, e.g., "A", "C", "G", "T" (or "U"), predicted to be in that position is in fact the nucleotide that belongs in that assigned position. Typical algorithms for performing alignment functions include Needleman-Wunsch and Smith-Waterman algorithms. In either case, these algorithms perform sequence alignments between a string of the subject’s query genomic sequence and a string of the reference genomic sequence whereby instead of comparing the entire genomic sequences, one with the other, segments of a selection of possible lengths are compared.

[0044] Once the reads have been assigned a position, such as relative to the reference genome, which may include identifying towhich chromosome the read belongs and / or its offset from the beginning of that chromosome, the readsmay be sorted by position. This may enable downstream analyses to take advantage of the oversampling procedures described herein. All of the reads that overlap a given position in the genome will be adjacent to each other after sorting and they can be organized into a pileup and readily examined to determine if themajority of them agreewith the reference value or not. If they do not, a variant can be flagged. 12 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55

[0045] For instance, in various embodiments, the methods of the disclosure may include generating a variant call file (VCF) identifying one or more, e.g., all, of the genetic variants in the individual who’s DNA / RNA were sequenced, e.g., relevant to one or more reference genomes. For instance, once the actual sample genome is known and compared to the referencegenome, thevariationsbetween the twocanbedetermined,anda list ofall thevariations / deviationsbetween the reference genome(s) and the sample genomemay be called out, e.g., a variant call file may be produced. Particularly, in one aspect, a variant call file containing all the variations of the subject’s genetic sequence to the reference sequence(s) may be generated.

[0046] Accordingly, a useful element of the methods and systems disclosed herein is a genomic reference from which mapping, aligning, variant calling, and other such processes of the systemmay be performed, such as in comparison to a referent. Typically, suchmaping, aligning, variant calling, and / or the likemay be performedwith respect to a single human reference, e.g., an "ideal reference" that is a composite of genetic code froma variety of different sources, and as such the typical reference genome doesn’t match any single person. Such secondary analysis leverages the fact that most people have a geneticmakeup that is very similar to the reference. Hence, although it is not perfect, the typical reference genome is useful in helping tomap and align reads to the right place in a person’s genome based on their general similarity with the reference.

[0047] The typical reference is also useful with respect to forming the pile ups, as discussed herein, of all the reads over their given mapped and / or aligned place(s) in the reference, which pileups thereby allow a greater amount of evidence to be considered when making a variant call at any given position. Particularly, the reference allows one to consider prior probabilities of what particular base at a particular position of a given read should be, as compared to the reference, when determiningwhat that baseactually iswith respect to the read.Hence, useof a referenceallows for theassumption that the identity of anybaseat anyposition in the reference iswhat is themost likely content of that baseof the read that is present in the human genome at that position. Accordingly, secondary analysis is usually performed in a manner so as to figure out how any given individual differs from the typical reference.

[0048] However, although employing a single reference is useful for determining the identity of any given base pair of a read of a subject, but in some instances, there may be significant differences between a given subject and a typical reference that is used when performing secondary processing of that particular subject’s DNA / RNA. Alternatively, there are someplaces in the typical reference that are problematic for amultiplicity of people, and in certain instances, there are significant differences from the reference that occur commonly in various parts of the population.

[0049] For instance, in some instances, theremaybe individual variants, e.g., single nucleotidepolymorphisms (SNPs), whichoccur in somesignificant portionof thepopulation, suchas3%or5%or10%of thepopulation, ormore tenpercent of the population. Particularly, in various instances, for any given individual there may be one or more segments of various subjects genome that has been replaced by another sequence of a similar or different length, with different content, of course. Further complicating matters is that this genetic re-arrangement may occur in a single copy of the chromosome. Hence, at one haplotype the subject’s DNA may be similar to that of the reference, while at the other haplotype, the subject’s DNA may be vastly different from that of the reference.

[0050] Consequently, in some places a subject’s DNA may be identical to the standard reference, and in some places dramatically different from the standard reference. In some instances, such genetic variations may occur in predictable positions in the genome, and in particular geographical places. In other instances, the variantsmay occur in amuch larger percent of thepopulation, suchas80%of thepopulation. In suchan instance, the referencegenomemayactually show the less common content at a given region of the genome. Hence, in certain instances, there may be large sections of the reference genome, e.g., which might be hundreds or thousands or even millions of basis long, that are significantly different from a large sample set of the population. The outcome, therefore, is that if only the standard reference is employed in performing a secondary analysis process then the accuracy of such secondary analysis, e.g., mapping, aligning, and / or variant calling may not be as accurate is it could be.

[0051] It will of course will be better for those whose genome most closely match that of the reference, versus those whose genome has significant variations therefrom. The accuracy of secondary, and consequently, tertiary processing may be improved, therefore, if a reference being employed in the analyses is better fitted to the subject’s whose DNA is being processed, such as more closely aligned with that of their family members, ancestry, and the like. There are a multiplicity of methods and / or strategies that may be employed so as to overcome these potential inefficiencies of performing secondary processing using a standard reference genome.

[0052] For instance, a first traditional standard, e.g., linear, reference genome may be employed for determining the genomic identity of one strand of a subject’s DNA, e.g., one haplotype, and a second traditional, or non-traditional, reference genomemay be employed for determining the genomic identity of the other strand of the subject’s DNA.Hence, there may be one reference sequence for chromosome one, and another reference sequence for chromosome two, where, in certain instances, the reference sequences may be generated and / or otherwise employed dynamically, e.g., based on auxiliary data, e.g., ancestry, of the subject or person. In such an instance, a first secondary processing procedure, e.g., mapping, alignment, and variant calling procedure, may be performed, e.g., based on a standard reference, and in a second processing procedure, a second reference genome, e.g., one that is ancestry specific, may be 13 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 employed in the secondary processing procedure(s).

[0053] This secondary processing procedure may be performed with respect to the entire genome of the subject, or for one or more identified regions thereof. For example, where a region by region secondary processing procedure is being performed, various genetic markers may be used to identify the regions for more careful processing. Particularly, once a region of variation is determined, for a subject’s genome, a given secondary reference may then be employed by the system for the performance of a secondary processing procedure with respect to that one or more segments.

[0054] In amanner such as this a plurality of referencesmay be used, where each reference is selected to enhance the accuracyand / or efficiencyof the secondary processingprocedurebeingperformed. In particular instances, therefore, one cultural reference, e.g., a European or African reference, may be employed for processing a given portion of a subject’s DNA, while another cultural reference, e.g., one or more Asian, Indian, South American reference, may be employed for processing another given portion of the subject’s DNA. Particularly, a database storing a multiplicity of references, e.g., specific to given populations and / or geographies, may be employed, such that at any given time the system may dynamically switch between what reference is to be employed for determining any given segment of a subject’s DNA.

[0055] Hence, in a particular use instance, a given long segment of a subject’s DNA, e.g., 1 million base pairs, may be analyzed employing a generated or standard European reference, and another long segment of DNA, e.g., 2million base pairs, may be analyzed by a generated or standard other, e.g., North American, reference. Particularly, a statistical analysis may be performed, e.g., to determine the percentage homology of any given portion of a genome to a particular reference standard, and which of a selection of reference standards is to be employed may be determined based on the outcome of the statistical analysis.

[0056] More particularly, the artificial intelligence module, discussed herein below, may be employed to determine the most relevant reference to use for performing secondary analysis on any given region of a subject’s DNA, e.g., so as to employ the reference that fits best. In various instances, a plurality of these standards may be mixed and matched in any logical order, so as to produce a combined, e.g., Chimeric, reference, which may be built of segments from a variety of sources.

[0057] As indicated, thismay be performed in a haploid or diploidmanner, such aswhere the reference is applied to only onecopy, e.g., strand, of theDNA,or toboth copiesof theDNA.Thismaybecomplicatedby the fact that different strandsof DNA of a subject from different sources may have different splicing patterns. These patterns, however, may be used as a map by which to differentially and / or dynamically employ reference genomes. Such splicing patterns may be based on ancestral genetic background. And, in some instances, based on these differential slice patterns, a chimeric reference genome, e.g., including different culturally relevant reference genomes.

[0058] Likewise, these references may then be used as a guide for the mapping, aligning, and / or variant calling procedures, such that use of a non-traditional, e.g., other standard or chimeric, reference genomewill give a closermatch to the actual genome of the user, and therefore will provide for amore accuratemapping, aligning, and / or variant calling of the subject’s genomic sequence. Thus, the overall analysis will have a higher probability of accuracy.

[0059] As indicated, the reference genome to be employed may be dynamic, and may be built, e.g., on the fly, to specifically, and more closely represent the genome of the subject. For instance, in various instances, the chimeric reference genomemay be assembled, such as in a De Bruijn graph format, where variations from the standard reference may be represented by bubbles in the graph, which bubbles may refer to various mapped coordinates in the standard reference. Particularly, a graph based reference may be generated, such that wherever a variation in the reference is to occur, the change from the standard may be represented as a bubble in a graph.

[0060] So, where a newly built, e.g., chimeric, reference is used, those regions where the chimeric reference matches with the standard reference, e.g., the backbone, may be represented as a straight line, but where the chimeric reference includes a differential segment, e.g., a branch, this difference may be represented as a bubble in the standard reference, e.g., where the bubble represents the different base pairs from the reference. The bubbles may be any length, and one regionof bubbles in the chimeric referenceneednot be the same lengthas theother.Hence, once the referencegenome is assembled, it may be backtraced and / or otherwise mapped to the traditional reference to track the manner in which the dynamic, e.g., chimeric, reference differs from the traditional reference.

[0061] In such an instance, a local assembly reference may be generated so as to accord with the specific ancestry and / or culture of the subject, such aswhere the bubble regions represent ancestral differences fromastandard reference. In a manner such as this, a dynamic reference may be generated where each reference employed is specific to the individual, and consequently, no two references will be alike.

[0062] Another manner in which a dynamic reference may be built and / or employed is to build a chimeric reference basedon knownpopulation variations, e.g., common for the detected ancestry and / or knowndifferent segment,where the standard reference is changed in various regions to include known segments of variations, e.g., known ancestral and / or cultural variations, which variations may then be annotated so as to build a map of the chimeric reference. In such an instance, when used, it may be known which reference segment from which source is being used when performing a mapping and / or aligning operation on the subject’s DNA and / or for determining how that DNA differs from the reference.

[0063] For instance, once it has been determined for a subject, which part of their DNA comes from which part of their 14 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 ancestry, a reference coherent with that ancestry, over the identified sequence length, may be employed as at least a segment of the chimeric reference. For example, a reference genome may be built based on known variations in human populations, such as based on geography, culture, ancestry, etc., such as where common alleles are known and may be used in producing a chimeric reference.

[0064] Specifically, where a sequence includes a plurality of SNPs in a row, a certain part of the population may have a certain order of the combination, and a certain other part of the population may have a different order of the combination, these variations may be represented as either annotations or as bubbles, such as in a De Bruijn graph format. These variations may be representative of different haplotypes of the population, where such common variations from the standard reference may be coded and represented as bubbles or annotations in the reference, e.g., graph, of variable lengths. In such an instance, typical variant callers will not distinguish between these differences, and will not be able to resolve this area of the genome.

[0065] However, using thedifferential referencegenomeof the system, these regionsmaybemoreaccurately resolved. In such an instance, it might be more useful to represent these variations as bubble instead of annotations as individual SNPs so that the difference is clear, since the SNPs are near each other or otherwise densely spaced. Hence, there are advantages to having bubbles, even longer bubbles to represent such variations. Consequently, the complete reference neednot benon-standard, in some instances, only various segments needbeswappedout, e.g., editedandannotated, so as to form the chimeric reference. Specifically, in certain instances, the differential segments need not be absolutely changed, the change may be made optional, e.g., variable, depending on how the system determines which reference code, traditional or variable, in which circumstances, any of which may be implemented in the hardwired configurations disclosed herein.

[0066] In suchamanner, for anyvariation in the reference, e.g., at anygivennucleotideposition, theremaybeavariation at that position between one nucleotide and another, the absolute determination of which variation depends on which reference is to be employed, and in certain instances, may be determined on the fly, such as during the analysis process. For instance, with respect to one reference genome, such as a dominant reference thatmatches a large percentage, e.g., 75% of the population, the dominant reference may indicate an "A" at a given position, while a sub-dominant reference, whichmaymatch a smaller percentageof the population, e.g. 25%,maydiffer from the dominant referenceby havinga "T" at that particular position. Hence, when employing only the dominant reference, a "non-match"may occur, but when using the non-dominant reference, a match may occur. Consequently, employing a plurality of references, or a chimeric reference, may lead to better accuracy.

[0067] In various instances, this known variability may simply be flagged or annotated as being a known variation in the population. Specifically, in various instances, these variables may be annotated by one or more flags so as to demarcate regions of variability within the reference. This is especially useful for determining one ormore SNPs. However, in various instances, thismay lead toanotherproblem, suchaswhere theremaybe threeSNPs ina row,suchasan "A" "A" "A",where at each position a known variable may be present and flagged, such as where the first "A" may alternatively be a "C", the second may be "T", and the third variable may alternatively be a "G". In such an instance, these three bases could be flagged as having three independent variables, but in some instances, each variation may not represent an independent SNP, butmay actually be amore common haplotype in the population. Hence, the first haplotypemay represent an "AAA" sequence, and the second haplotype may represent an "CTG," in which case these variations do not sort randomly, but rather collectively. That is part of the population has an "AAA" haplotype, while another part has a "CTG" haplotype.

[0068] In such instances, rather than flag each base individually as a variable, it may be useful to indicate the variation collectively as a "bubble" in the reference graph. Accordingly, in various instances, one or more segments of the genome may fromhaplotypes that are very similar or even identical to oneanother. As such, oneormore reads in thegenomeof the subject may correspond to these one or more haplotypes in a primary or secondary assembly. Using a typical reference, such a read covering oneof these haplotypes, in a conventional system,will not bemapped or aligned because itmatches to too many different positions.

[0069] Specifically, in various instances, a read from a subject could correspond ormatch to one particular haplotype or another or maymatch to the primary assembly. In various instances, the read maymatch in all these places substantially equally well. The typical mapper, however, may not be able to resolve this difference. This situation may be overcome by having the mapper simply choose one position over another, in which case the odds of being correct decline with the number of potential matching positions, or it maymap to any and all overlapping positions, but thismay lead to a decrease in resolution.

[0070] Not mapping the sequence leaves viable information unaccounted for. One way to overcome this dilemma is for the system to use alternative references that show regions of variable haplotypes to which known reads containing the variant haplotype configurations may be mapped and aligned. As indicated above, a graph based mapper may be employed to indicate known alternative haplotype variations.

[0071] Specifically, in such an instance, the system may be configured to perform an alt aware type analysis. For instance,where various reads fromasubject are identical or substantially identical, a branchedgraphof the referencemay be generated to indicate the presence of alternative haplotypes, such aswhere each haplotype forms a different branch in 15 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 the graph. The branch, or bubble, may be longer or shorter as is required to meet the length of the haplotype sequence. Additionally, the number of branchesmay vary based on the number of known variant haplotypes there are, which number may be in the tens, hundreds, or more.

[0072] The system, therefore, may be configured such that the mapper will understand that each branch represents a potential alternate haplotype, in comparison to the primary assembly backbone. Anotherway to overcome this dilemma is for the system to take such substantially identical read sequences and consider them as a "new" chromosome. Specifically, the system may be configured to treat an alternative haplotype, e.g., which alternative has significant difference from the traditionally employed reference, as an entirely new chromosome by which to examine potential candidate sequences.

[0073] Suchaconfiguration is useful because it reduces falsepositives byassuming readsand / or the seeds thereof that don’t match the primary referencemay in fact align to the alternative haplotype. Particularly, without access to a reference including alternative haplotypes, such sequences may be force fit into a primary reference where they do not actually fit, resulting in a false positive, e.g., for aSNP, being called.However, in various instances, analternate haplotypemayhavea sequence that is quite long, and in various instances,may have portions thatmatch the primary reference. Thismay result in a read that appears to match both the primary and the haplotype reference.

[0074] In suchasituation the readmaynotbeable tobemapped, or itmaysimplybe randomlyassigned toone reference or the other, in which case coverage is reduced by 50%, assuming it has an equal chance of matching either reference, resulting in a lowerMAPQ because the two references now become in competition for one another. However, themapper maybeconfigured soas tobeAlt-aware, suchasbyemploying agraphbasedbackbonebywhich to placeboth references so as to not be in competition with one another with respect to determining best fit. Consequently, the mapper may be adapted such that it understands that a branch in the chain backbone represents analternative sequence that is related to, e.g., branched off from, the graph of the primary assembly, as such the two references will not be in competition with one another.

[0075] One way to accomplish this functionality is to employ a hash table that is adapted so as to be populated with the substantially similar reads, such as in accordancewith the hash table basedmapper disclosed above, but in this instance, a virtual, e.g., chimeric, referencemaybeemployedas the index. For instance, known variations, suchas knownalternate haplotypesequences,maybe includedwithinand / oremployedas the index, andmaybeused in thepopulationof thehash table, such as where the identified alternate haplotypes are entered into the hash table, e.g., as a virtual reference index, for seed mapping purposes.

[0076] In a manner such as this, matches in those positions may be identified, so as to improve the sensitivity of the system, and allowing reads that would otherwise remain unresolved, e.g., due to alternate haplotypes, to be resolved. Thus, the relationship of substantially identical haplotype reads, whichmay otherwisemap to the primary assembly, but in actuality do not belong there, may be determined. Hence, the mapper may be configured to take responsibility for the sorting and finding of the bestmatch in analternate haplotype, and then remapping it to its identified, e.g., lift-over, position in the primary assembly graph. Hence, in various instances, the virtual reference may be employed as a graph and / or branchoffof the reference,e.g., built upfront into themapper configuration, andmappingmayoccurasdescribedabove for pre-fix and suffix tree mapping.

[0077] Thesemethods provide for enhanced sensitivity and increased accuracy of the systemoverall, e.g., with respect to mapping and / or aligning, such as by minimizing false positive when substantially identical reads are not mapped, randomly mapped, or mapped to a multiplicity or wrong positions. Accordingly, in various embodiments, as described herein, a dynamic reference based system may be configured so as to employ a multiple graph branch configuration to map multiple substantially identical sequences that often occur non-randomly in a population, such as by employing a population significant and / or chimeric reference genome. And, as population studies increase, and more and more population related data is employed to build chimeric reference genomes, the accuracy of this system will continue to improve. Changes in the building of such graphs and / or tables may be informed by the changes in these population data, such as by accommodating ever increasing branches or bubbles in the graph and / or the number of alternate haplotypes available for consideration.

[0078] In various embodiments, a super dynamic reference may be generated, such as where the reference is specializing particularly to a specific community or family or event to the individual subject themselves, such as based on the subject’s specific ancestry. Accordingly, in accordance with the methods disclosed herein, the system may be configured for performing a first analysis, employing a standard reference, andmay further be configured for performing a second analysis employing a non-standard or modified, e.g., specialized, reference.

[0079] For instance, a first passmay be performedwith regard to the standard reference, the subject’s ancestrymay be determined, or othermarkers, e.g., geneticmarkers, identified, haplotypic informationmaybe identified, and / or a chimeric reference, e.g., includinghaplotype information,maybeassembled,whichchimeric referencemay thenbeusedwithin the system for purposes of mapping and / or aligning, such as when building the hash table.

[0080] Specifically, thechimericassembly canbut neednot bebuilt fromscratch.Rather, identifiedhaplotypessimplybe insertedor otherwisesubstitutedwithin themain referencebackbonesuchaswhere their branchchainwould indicate they 16 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 be inserted, and this referencemay then be inserted into the hash table for hashing.Hence, the chimeric referencemaybe built it not by completely replacing segments but by substituting segments of specific ancestral references, e.g., lift-over sequences, and listing or flagging them as alternate haplotypes for substitution into the primary reference.

[0081] For example, whether a seed of a readmaps to a non-chimeric or chimeric, e.g., annotated, reference segment, this informationmay be included, such as by an appropriate annotation, within the hash table. Particularly, the information tobe includedwithin thehash tablemay indicate that the referenceand / or read / seed / Kmer is annotated, that the reference is primary, and / or that one or more alt. haplotypes are included and / or matching, and / or that one or more lift-over groups, e.g., a lift-over seed group, are included, and the like. The actual candidates, therefore, may be in a lift-over group, where each lift-over group may be assigned a score, e.g., of the best representative, and the primary alignment, e.g., MAPQ, of this group may be reported, with respect to the difference in score from the second best group.

[0082] Specifically, it is useful to determine how the best lift-over group scored, as well as the distance in score from the second best lift-over group, which if the distance in score is substantial indicates a higher confidence of a correct match, regardless of how close the MAPQ scores are with regard to the sequence(s) in question matching the primary and alt. references. The system, therefore,may be configured to keep track of all of the annotations, to build the hash table, and to implement the hash function, score the results, as well as tomap and align the best results, e.g., in a pipeline fashion, and thus, keeping the primary reference as a backbone in building a dynamic reference is an important feature for facilitating the extensive bookkeeping that allows the subsequent functions to work efficiently and with better accuracy.

[0083] In amanner such as this, two ormore seeds thatmatch each other reasonablywell, but do not necessarilymatch the primary reference, need not be discarded if they match an alt. reference segment. In such an instance, they may be grouped together as alt. seeds.

[0084] Accordingly, the hash table may employ one or more of these techniques to recognize the various possible organizational structures of the seeds as well as their positions corresponding to either the ALT haplotype or primary assemblies, may organize them as such (e.g., Alt, Alt, Primary, etc.), and annotate them, e.g., some as being from the alt. and some from the primary assembly, etc., in the organizational structure of the hash table so as to ensure any relevant information theycontain isnot lostbut isuseable. In various instances, this informationand / ororganizational structuremay be employed by the mapper and carried over to the aligner.

[0085] In manners such as these, one or more of SW, HMM, and / or variant calling may be performed against the primary / chimeric reference, without having to juggle between alternative references and / or competing coordinates thereof, resulting in a more normalized coverage, better sensitivity, and a clear MAPQ. Likewise, the output file may be in any suitable file format, such as a typical BAM and / or SAM file (e.g., an altBAM / SAM file), and / or may bemodified to indicate the reference was chimeric and / or which haplotype sequences were implemented in the reference, e.g., an indication may be made for which haplotype was included within the primary reference and where, what coordinates (‑ such as a lift-over map), and which sequences mapped to the haplotype as compared to the primary reference, and the like. In various instances, it may then be useful to include this seed group as a lift-over position in the chimeric reference.

[0086] Specifically, in the context of using a graph based, dynamic reference, as herein disclosed, a more sensitive mapping and / or aligning may be performed resulting in better accuracy, where the graph indicates how the dynamic referencewas stitched together and / or how the subject’s genetic sequencemapped thereto. Further, as indicated in detail above, this dynamic referencemaybe implemented in optimized software, suchasbyperformancebyaCPUand / orGPU, or may be implemented in hardware, such as by an integrated circuit, e.g., FPGA, ASIC, or the like, of the disclosure.

[0087] Hence, in particular embodiments, a platform of technologies for performing genetic analyses are provided where theplatformmay include theperformanceof oneormoreof:mapping, aligning, sorting, local realignment, duplicate marking, base quality score recalibration, variant calling, compression, and / or decompression functions. For instance, in various aspects a pipeline may be provided wherein the pipeline includes performing one or more analytic functions, as describedherein, onagenomicsequenceof oneormore individuals, suchasdataobtained inan imagefileand / oradigital, e.g., FASTQorBCL, file format fromanautomatedsequencer.A typical pipeline tobeexecutedmay includeoneormoreof sequencing genetic material, such as a portion or an entire genome, of one or more individual subjects, which genetic material may include DNA, ssDNA, RNA, rRNA, tRNA, and the like, and / or in some instances the genetic material may represent coding or noncoding regions, such as exomes and / or episomes of the DNA. The pipeline may include one or more of performing an image processing procedure, a base calling and / or error correction operation, such as on the digitized genetic data, and / ormay include oneormore of performing amapping, analignment, and / or a sorting function on the genetic data. In certain instances, the pipelinemay include performingoneormore of a realignment, a deduplication, a basequality or score recalibration, a reduction and / or compression, and / or a decompression on the digitized genetic data. In certain instances thepipelinemay includeperformingavariant callingoperation, suchasaHiddenMarkovModel, on the genetic data.

[0088] Accordingly, in certain instances, the implementation of oneormoreof theseplatform functions is for thepurpose of performing one or more of determining and / or reconstructing a subject’s consensus genomic sequence, comparing a subject’sgenomicsequence toa referent sequence, e.g., a referenceormodel genetic sequence, determining themanner in which the subject’s genomic DNA or RNA differs from a referent, e.g., variant calling, and / or for performing a tertiary 17 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 analysis on the subject’s genomic sequence, such as for genome-wide variation analysis, gene function analysis, protein function analysis, e.g., protein bindinganalysis, quantitative and / or assembly analysis of genomesand / or transcriptomes, as well as for various diagnostic, and / or a prophylactic and / or therapeutic evaluation analyses.

[0089] As indicated above, in one aspect one or more of these platform functions, e.g., mapping, aligning, sorting, realignment, duplicate marking, base quality score recalibration, variant calling, compression, and / or decompression functions is configured for implementation in software. In some aspects, one or more of these platform functions, e.g., mapping, aligning, sorting, local realignment, duplicatemarking, base quality score recalibration, decompression, variant calling, compression, and / or decompression functions is configured for implementation in hardware, e.g., firmware. In certain aspects, these genetic analysis technologies may employ improved algorithms that may be implemented by software that is run in a less processing intensive and / or less time-consuming manner and / or with greater percentage accuracy, e.g., the hardware implemented functionality is faster, less processing intensive, and more accurate.

[0090] In particular, where the algorithm is to be implemented in a software solution, the algorithm and / or its attendant processes, has been optimized so as to be performed faster and / or with better accuracy for execution by that media. Likewise, where the functions of the algorithm are to be implemented in a hardware solution, e.g., as firmware, the hardware has been designed to perform these functions and / or their attendant processes in an optimizedmanner so as to be performed faster and / or with better accuracy for execution by that media. Further, where the algorithm is to be implemented in a quantumprocessing solution, the algorithmand / or its attendant processes, has been optimized so as to be performed faster and / or with better accuracy for execution by that media. These methods, for instance, can be employed such as in an iterative mapping, aligning, sorting, variant calling, and / or tertiary processing procedure. In another instance, systems and methods are provided for implementing the functions of one or more algorithms for the performance of one ormore steps for analyzing genomic data in a bioinformatics protocol, as set forth herein, wherein the functions are implemented on a hardware and / or quantumaccelerator, whichmay ormaynot be coupledwith one ormore general purpose processors and / or super computers and / or quantum computers.

[0091] In one aspect, in various embodiments, once the subject’s genome has been reconstructed and / or a VCF has been generated, such datamay then be subjected to tertiary processing so as to interpret it, such as for determining what the data means with respect to identifying what diseases this personmay or may have the potential for suffer from and / or for determining what treatments or lifestyle changes this subject may want to employ so as to ameliorate and / or prevent a diseased state. For example, the subject’s genetic sequence and / or their variant call file may be analyzed to determine clinically relevant genetic markers that indicate the existence or potential for a diseased state and / or the efficacy of a proposed therapeutic or prophylactic regimenmay have on the subject. This datamay then be used to provide the subject with one or more therapeutic or prophylactic regimens so as to better the subject’s quality of life, such as treating and / or preventing a diseased state.

[0092] Particularly, once one or more of an individual’s genetic variations are determined, such variant call file information can be used to develop medically useful information, which in turn can be used to determine, e.g., using various known statistical analysis models, health related data and / or medical useful information, e.g., for diagnostic purposes, e.g., diagnosingadiseaseorpotential therefore, clinical interpretation (e.g., looking formarkers that represent a disease variant), whether the subject should be included or excluded in various clinical trials, and other such purposes. More particularly, in various instances, the generated genomics and / or bioinformatics processed results data may be employed in the performance of one or more genomics and / or bioinformatics tertiary protocols, such as a micro-array analysis protocol, a genome, e.g., whole genome analysis protocol, a genotyping analysis protocol, an exome analysis protocol, anepigenomeanalysis protocol, ametagenomeanalysis protocol, amicrobiomeanalysis protocol, a genotyping analysis protocol, including joint genotyping, variants analyses protocols, including structural variants, somatic variants, and GATK, as well as RNA sequencing protocols and other genetic analyses protocols.

[0093] As there are a finite number of diseased states that are caused by genetic malformations, in tertiary processing variants of a certain type, e.g., those known to be related to the onset of diseased states, can be queried for, such as by determining if oneormoregenetic baseddiseasedmarkersare included in thevariant call fileof thesubject.Consequently, in various instances, the methods herein disclosed may involve analyzing, e.g., scanning, the VCF and / or the generated sequence, against a known disease sequence variant, such as in a data base of genomic markers therefore, so as to identify the presence of the genetic marker in the VCF and / or the generated sequence, and if present to make a call as to the presence or potential for a genetically induced diseased state. Since there are a large number of known genetic variationsanda largenumberof individual’s suffering fromdiseasescausedbysuchvariations, in someembodiments, the methods disclosed herein may entail the generation of one or more databases linking sequenced data for an entire genome and / or a variant call file pertaining thereto, e.g., such as from an individual or a plurality of individuals, and a diseased state and / or searching the generated databases to determine if a particular subject has a genetic composition that would predispose them to having such diseased state. Such searching may involve a comparison of one entire genome with one or more others, or a fragment of a genome, such as a fragment containing only the variations, to one or more fragments of one or more other genomes such as in a database of reference genomes or fragments thereof.

[0094] Therefore, in various instances, a pipeline of the disclosure may include one or more modules, wherein the 18 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 modules are configured for performing one or more functions, such as an image processing or a base calling and / or error correction operation and / or amapping and / or an alignment, e.g., a gapped or gapless alignment, and / or a sorting function on genetic data, e.g., sequenced genetic data. And in various instances, the pipeline may include one or more modules, wherein themodules are configured for performing onemore of a local realignment, a deduplication, a base quality score recalibration, a variant calling, e.g., HMM, a reduction and / or compression, and / or a decompression on the genetic data. Additionally, the pipeline may include one or more modules, wherein the modules are configured for performing a tertiary analysis protocol, such as micro-array protocols, genome, e.g., whole genome protocols, genotyping protocols, exome protocols, epigenome protocols, metagenome protocols, microbiome protocols, genotyping protocols, including joint genotyping protocols, variants analysis protocols, including structural variants protocols, somatic variants protocols, and GATK protocols, as well as RNA sequencing protocols and other genetic analyses protocols.

[0095] Many of these modules may either be performed by software or on hardware, locally or remotely, e.g., via software or hardware, such as on the cloud, e.g., on a remote server and / or server bank, such as a quantum computing cluster. Additionally,many of thesemodules and / or steps of the pipeline are optional and / or can be arranged in any logical order and / or omitted entirely. For instance, the software and / or hardware disclosed herein may or may not include an imageprocessingand / or a base callingor sequencecorrection algorithm, suchaswhere theremaybeaconcern that such functionsmay result in a statistical bias. Consequently, the systemmay include ormay not include the base calling and / or sequence correction function, respectively, dependent on the level of accuracy and / or efficiency desired. And as indicated above, one ormore of the pipeline functionsmay be employed in the generation of a genomic sequence of a subject such as through a reference based genomic reconstruction. Also, as indicated above, in certain instances, the output from the secondaryprocessingpipelinemaybeavariant call file (VCF,gVCF) indicatingaportionor all thevariants in agenomeora portion thereof.

[0096] For instance, in various embodiments, a Next Generation sequencer, or a sequencer on a chip technology, may be configured to perform a sequencing operation on received genetic data. For instance, as can be seen with respect to FIG. 1A, the genetic data 6a may be coupled to a sequencing platform 6 for insertion into a Next Gen sequencer to be sequenced in an iterative fashion, such that each sequence will be grown by the stepwise addition of one nucleotide after another. Specifically, the sequencing platform 6 may include a number of template nucleotide sequences 6a from the subject that are arranged in a grid like fashion to form tiles 6b on the platform 6, which template sequences 6a are to be sequenced. The platform 6may be added to a flow cell 6c of the sequencer that is adapted for performing the sequencing reactions.

[0097] As the sequencing reactions take place, at each step a nucleotide having a fluorescent tag 6d is added to the platform6of the flowcell 6c. If a hybridizing reaction occurs, fluorescence is observed, an image is taken, the image is then processed, and an appropriate base call is made. This is repeated base by base until all of the template sequences, e.g., the entire genome, has been sequenced and converted into reads, thereby producing the read data of the system.Hence, once sequenced, the generated data, e.g., reads, need to be transferred from the sequencing platform into the secondary processing system. For instance, typically, this image data is converted into a BCL and / or FASTQ file that can then be transported into the system.

[0098] However, in various instances, this conversion and / or transfer processmay bemademore efficient. Specifically, presented herein are methods and architectures for expedited BCL conversion into files that can be rapidly processed within the secondary processing system. More specifically, in particular instances, instead of transmitting the raw BCL or FASTQ files, the images produced representing each tile of the sequencing operationmay be transferred directly into the system and prepared for mapping and aligning et al. For instance, the tiles may be streamed across a suitably configured PCIe and into the ASIC, FPGA, or QPU, wherein the read data may be extracted therefrom directly, and the reads advanced into the mapping and aligning and / or other processing engines.

[0099] Particularly, with respect to the transfer of the data from the tiles obtained by the sequencer to the FPGA / C- PU / GPU / QPU, as can be seenwith respect to FIG. 1A, the sequencing platform6may be imaged as a 3-D cube 6c, within which the growing sequences 6a are generated. Essentially, as can be seen with respect to FIG. 1B, the sequencing platform6maybe composedof 16 lanes, 8 in the front and 8 in the back,whichmaybe configured to formabout 96 tiles 6b. Within each tile 6b are a number of template sequences 6a to be sequenced thereby forming reads, where each read represents the nucleotide sequence for a given region of the genomeof a subject, each column represents one file, and as digitally encoded represents1byte foreveryfile,with8bitsper file, suchaswhere2bits represents thecalledbase, and the remaining 6 bits represents the quality score.

[0100] More particularly, with respect to Next Gen Sequencing, the sequencing is typically performed on glass plates 6 that form flow cells 6c that are entered into the automated sequencer for sequencing. As can be seen with respect to FIG. 1B, a flowcell 6c is aplatform6composedof 8 vertical columnsand8horizontal rows (front andback), togetherwhich form 16 lanes, where each lane is sufficient for the sequencing of an entire genome. TheDNA and / or RNA 6a of a subject to be sequenced is associated within designated positions in between fluidly isolated intersections of the columns and rows of the platform 6 so as to form the tiles 6b, where each tile includes template genetic material 6a to be sequenced. The sequencingplatform6, therefore, includesanumberof templatenucleotidesequences from thesubject,whichsequences 19 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 arearranged inagrid like fashionof tileson theplatform. (SeeFIG.1B.)Thegenetic data6 is thensequenced inan iterative fashionwhereeachsequence is grownby thestepwise introductionof onenucleotideafter another into theflowcell,where each iterative growth step represents a sequencing cycle.

[0101] As indicated, an image is captured after each step, and the growing sequence, e.g., of images, form the basis by which the BCL file is generated. As can be seen with respect to FIG. 1C, the reads from the sequencing procedure may form clusters, and it is these clusters that form the theoretical 3-D cube 6c. Accordingly, within this theoretical 3-D cube, each base of each growing nucleotide strand being sequenced will have an x dimension and a y dimension. The image data, or tiles 6b, from this 3-D cube 6cmay be extracted and compiled into a two-dimensionalmap, fromwhich amatrix, as seen in FIG.1ADmay be formed. Thematrix is formed of the sequencing cycles, which represent the horizontal axis, and the read identities,which represent the vertical axis. Accordingly, as canbe seenwith reference toFIG. 1C, the sequenced reads form clusters in the flow cell 6c, which clusters may be defined by a vertical and horizontal axis, cycle by cycle, and the base by base data from each cycle for each read may be inserted into the matrix of FIG. 1D, such as in a streaming and / or pipelined fashion.

[0102] Specifically, each cycle represents the potential growth of each read within the flow cell by the addition of one nucleotide, which when sequencing one or several human genomes, may represent the growth of about 1 billion or more reads per lane. The growth of each read, e.g., by the addition of a nucleotide base, is identified by the iterative capturing of images, of the tiles 6b, of the flow cell 6c in between the growth steps. From these images base calls aremade, and quality scores determined, and the virtualmatrix of FIG1D is formed.Accordingly, therewill be both abase call andaquality score entered into the matrix, where each tile from each cycle represents a separate file. It is to be noted that where the sequencing is performed on an integrated circuit, sensed electronic data may be substituted for the image data.

[0103] For instance, as can be seen with respect to FIG. 1D, the matrix itself will grow iteratively as the images are captured and processed, bases are called, and quality scores are determined for each read, cycle by cycle. This is repeated for eachbase in the read, for each tile of the flowcell. For example, the cluster of reads. 1Cmaybenumberedand entered into thematrix as the vertical axis. Likewise, the cycle numbermaybe entered as the horizontal axis, and the base call and quality scoremay then be entered so as to fill out thematrix columnby column, rowby row. Accordingly, each read will be represented by a number of bases, e.g., about 100 or 150 up to 1000 or more bases per read depending on the sequencer, and there may be up to 10million or more reads per tile. So, if there are about 100 tiles each having 10million reads, the matrix would contain about 1 billion reads, which need to be organized and streamed into the secondary processing apparatus.

[0104] Accordingly, such organization is fundamental to rapidly and efficiently processing the data. Hence, in one aspect, presented herein are methods for transposing the data represented by the virtual sequencing matrix in a manner so that the data may be more directly and efficiently streamed into the pipelines of the system herein disclosed. For instance, the generation of the sequencing data, as represented by the star cluster of FIG. 1C, is largely unorganized, which isproblematic fromadataprocessingstandpoint.Particularly, as thedata isgeneratedby thesequencingoperation, it is organizedasone file per cycle,whichmeans that by theendof the sequencingoperation there aremillionsandmillions of files generated, which files are represented in FIG. 1E, by the data in the columns, demarcated by the solid lines. However, for the purposes of secondary and / or tertiary processing, as disclosed herein, the file data needs to be re- organized into read data, demarcated by the dashed lines of FIG. 1E.

[0105] More particularly, in order to more efficiently stream the data generated by the sequencer into the secondary processingdata, thedata representedby thevirtualmatrix shouldbe transposed, suchasby reorganizing thefiledata from acolumnby columnbasis of tiles per cycle, to a rowby rowbasis identifying thebasesof eachof the reads. Specifically, the data structureof the generated files forming thematrix, as it is producedby the sequencer, is organizedonacycle by cycle, column by column, basis. By the processes disclosed herein, this data may be transposed, e.g., substantially simulta- neously, so as to be represented, as seen within the virtual matrix, on a read by read, row by row basis, where each row represents an individual read, and each read is represented by a sequential number of base calls and quality scores, thereby identifying both the sequence for each read and its confidence. Thus, in a transpose operation as herein described, the datawithin thememorymay be re-organized, e.g., within the virtualmatrix, froma columnby columnbasis, representing the input data order, to a row by row basis, representing the output data order, thereby transposing the data order fromavertical to ahorizontal organization. Further, although theprocessmaybe implementedefficiently in software, it may be made even more efficiently and faster, by being implemented in hardware and / or by a quantum processor.

[0106] For instance, in various instances, this transposition process may be accelerated by being implemented in hardware.Forexample, inone implementation, inafirst step, thehost software, e.g., of thesequencer,maywrite inputdata into thememory, associatedwith theFPGA,onacolumnbycolumnbasis, e.g., in the input order.Specifically, as thedata is generated and stored into an associated memory, the data may be organized into files, cycle by cycle, where the data is saved as separate individual files. This data may be represented by the 3-D cube of FIG. 1A. This generated column organized data may then be queued and / or streamed, e.g., in flight, into the hardware where dedicated processing engines will queue up the column organized data and transpose that data from a column by column, cycle order configuration, to a row by row, read order configuration, in a manner as described herein above, such as by converting 20 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 the 3-D tile data into a 2-Dmatrix,whereby the columndatamaybe reorganized into rowdata, e.g., on a read to readbasis. This transposed data may then be stored in the memory in a more strategic order.

[0107] For example, the host software may be configured to write input data into the memory associated with the chip, e.g., FPGA, such as in a column-wise input order, and likewise the hardware may be configured to queue the data in a manner so that it is red into thememory in a strategic manner, such as set forth in FIG. 1F. Specifically, the hardwaremay include an array of registers 8a intowhich the cycle filesmaybe dispersedand re-organized into individual read data, such as by writing one base from a column into registers that are organized into rows. More specifically, as can be seen with respect to FIG. 1G, the hardware device 1, including the transposition processing engine 8, may include a DRAM port 8a that may queue up the data to be transposed, where the port is operably coupled to a memory interface 8b that is associatedwithaplurality of registersand / oranexternalmemory8c,and is configured forhandlingan increasedamountof transactions per cycle, where the queued data is transmitted in bursts.

[0108] Particularly, this transposition may take place one data segment at a time, such as where thememory accesses are queued up in such amanner as to takemaximal advantage of theDDR transmission rate. For instance, with respect to DRAM, theminimal burst length of the DDRmay be, for example, 64 bytes. Accordingly, the column arranged data stored in the host memory may be accessed in a manner such that with each memory access a column worth of corresponding, e.g., 64, bytes of data is obtained. Hence, with one access of the memory a portion of a tile, e.g., representing a corresponding "64" cycles or files, may be accessed, on a column by column basis.

[0109] However, as can be seen with respect to FIG. 1F, although the data in the host memory is accessed as column data,when transmitted to thehardware, itmaybeuploaded into associated smallermemories, e.g., registers, in a different orderwhereby thedatamaybeconverted intobytes, e.g., 64bytes, of rowby row readdata, suchas in accordancewith the minimal burst rate of the DDR, so as to generate a corresponding "64" memory units or blocks per access. This is exemplified by the virtual matrix of FIG. 1Dwhere a number of reads, e.g., 64 reads, are accessed in blocks, and read into memory in segments, as represented by FIG. 1E, such as where each register, or flip-flop, accounts for a particular read, e.g., 64 cycles x 64 reads x 8 bits per read = 32K flip-flops. Specifically, thismay be accomplished in various different ways inhardware, suchaswhere the inputwiring isorganized tomatch thecolumnordering, and theoutputwiring isorganized to match the row order. Hence in this configuration, the hardware may be adapted so as to both read and / or write to "64" different addresses per cycle.

[0110] More particularly, the hardware may be associated with an array of registers such that each base of a read is directed and written into a single register (or multiple registers in a row) such that when each block is complete, the newly ordered rowdatamay be transmitted tomemory as an output, e.g., FASTQdata, in a rowby roworganization. TheFASTQ data may then be accessed by one or more further processing engines of the secondary processing system for further processing, such as by a mapping, aligning, and / or variant calling engine, as described herein. It is to be noted, as described herein, the transpose is performed in small blocks, however, the systemmay be adapted for the processing of larger blocks as well, as the case may be.

[0111] As indicated, once aBCL file has been converted into a FASTQ file, as described above, and / or a BCL or FASTQ file has otherwise been received by the secondary processing platform, a mapping operation may be performed on the receiveddata.Mapping, in general, involvesplotting the reads toall the locations in the referencegenome towhere there is a match. For example, dependent on the size of the read there may be one or a plurality of locations where the read substantially matches a corresponding sequence in the reference genome. Hence, the mapping and / or other functions disclosedhereinmaybe configured for determiningwhereout of all the possible locations oneormore readsmaymatch to in the reference genome is actually the true location to where they map.

[0112] The output returned from the performance of a mapping function may be a list of possibilities as to where one or more, e.g., each, readmaps to one or more reference genomes. For instance, the output for eachmapped readmay be a list of possible locations the read may be mapped to a matching sequence in the reference genome. In various embodiments, an exact match to the reference for at least a piece, e.g., a seed of the read, if not all of the read may be sought. Accordingly, in various instances, it is not necessary for all portions of all the reads to match exactly to all the portions of the reference genome.

[0113] More particularly, in various instances, amappingmodulemay be provided, such aswhere themappingmodule is configured to perform one or more mapping functions, such as in a hardwired configuration. Specifically, the hardwired mappingmodulemaybeconfigured toperformoneormore functions typically performedbyoneormorealgorithms runon a CPU, such as the functions that would typically be implemented in a software based algorithm that produces a prefix and / or suffix tree, aBurrows-Wheeler Transform, and / or runsahash function, for instance, a hash function thatmakesuse of, or otherwise relies on, a hash-table indexing, such as of a reference, e.g., a reference genome sequence. In such instances, the hash functionmay be structured so as to implement a strategy, such as an optimizedmapping strategy that may be configured tominimize the number of memory accesses, e.g., large-memory random accesses, being performed soas to therebymaximize theutility of theon-boardorotherwiseassociatedmemorybandwidth,whichmay fundamentally be constrained such as by space within the chip architecture.

[0114] It hasbeendeterminedwhereall thepossiblematchesare for theseedsagainst the referencegenome, itmust be 21 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 determinedwhich out of all the possible locations a given readmaymatch to is in fact the correct position towhich it aligns. Hence, after mapping there may be a multiplicity of positions that one or more reads appear to match in the reference genome. Consequently, theremay be a plurality of seeds that appear to be indicating the exact same thing, e.g., theymay match to the exact same position on the reference, if you take into account the position of the seed in the read. The actual alignment, therefore, must be determined for each given read. This determinationmay bemade in several different ways.

[0115] In one instance, all the reads may be evaluated so as to determine their correct alignment with respect to the reference genome based on the positions indicated by every seed from the read that returned position information during the mapping, e.g., hash lookup, process. However, in various instances, prior to performing an alignment, a seed chain filtering function may be performed on one or more of the seeds. For instance, in certain instances, the seeds associated with a given read that appear tomap to the samegeneral place asagainst the referencegenomemaybeaggregated into a single chain that references the samegeneral region.All of the seedsassociatedwith one readmaybegrouped intooneor more seed chains such that each seed is a member of only one chain. It is such chain(s) that then cause the read to be aligned to each indicated position in the reference genome.

[0116] Specifically, in various instances, all the seeds that have the same supporting evidence indicating that they all belong to the same general location(s) in the referencemay be gathered together to form one or more chains. The seeds that group together, therefore, or at least appear as they are going to be near one another in the reference genome, e.g., within a certain band, will be grouped into a chain of seeds, and those that are outside of this band will be made into a different chain of seeds.Once these various seeds have been aggregated into one ormore various seed chains, it may be determinedwhichof thechainsactually represents thecorrect chain tobealigned.Thismaybedone,at least inpart, byuse of a filtering algorithm that is a heuristic designed to eliminate weak seed chains which are highly unlikely to be the correct one.

[0117] The outcome from performing one or more of these mapping, filtering, and / or editing functions is a list of reads which includes for each read a list of all the possible locations to where the readmaymatchupwith the reference genome. Hence, amapping functionmaybe performed so as to quickly determinewhere the reads of the image file, BCL file, and / or FASTQ file obtained from the sequencer map to the reference genome, e.g., to where in the whole genome the various reads map. However, if there is an error in any of the reads or a genetic variation, you may not get an exact match to the reference and / or theremaybe several places oneormore reads appear tomatch. It, therefore,must bedeterminedwhere the various reads actually align with respect to the genome as a whole.

[0118] Accordingly, after mapping and / or filtering and / or editing, the location positions for a large number of reads have been determined, where for some of the individual reads a multiplicity of location positions have been determined, and it now needs to be determined which out of all the possible locations is in fact the true or most likely location to which the various reads align. Such aligning may be performed by one or more algorithms, such as a dynamic programming algorithm thatmatches themapped reads to the reference genomeand runs analignment function thereon. Anexemplary aligning function compares one or more, e.g., all of the reads, to the reference, such as by placing them in a graphical relation to one another, e.g., such as in a table, e.g., a virtual array or matrix, where the sequence of one of the reference genome or the mapped reads is placed on one dimension or axis, e.g., the horizontal axis, and the other is placed on the opposed dimensions or axis, such as the vertical axis. A conceptual scoringwave front is then passed over the array so as to determine the alignment of the readswith the reference genome, suchasby computing alignment scores for each cell in the matrix.

[0119] Thescoringwave front represents oneormore, e.g., all, the cells of amatrix, or a portionof those cells,whichmay be scored independently and / or simultaneously according to the rules of dynamic programming applicable in the alignment algorithm, such as Smith-Waterman, and / or Needleman-Wunsch, and / or related algorithms. Alignment scores may be computed sequentially or in other orders, such as by computing all the scores in the top row from left to right, followed by all the scores in the next row from left to right, etc. In this manner the diagonally sweeping diagonal wave front representsanoptimal sequenceofbatchesof scorescomputedsimultaneouslyor inparallel in aseriesofwave front steps.

[0120] For instance, in oneembodiment, awindowof the referencegenomecontaining the segment towhich a readwas mappedmaybeplacedon thehorizontal axis, and the readmaybepositionedon the vertical axis. In amanner suchas this an array ormatrix is generated, e.g., a virtualmatrix, whereby the nucleotide at each position in the readmaybe compared with the nucleotide at each position in the reference window. As the wave front passes over the array, all potential ways of aligning the read to the referencewindoware considered, including if changes toone sequencewouldbe required tomake the readmatch the reference sequence, such as by changing one or more nucleotides of the read to other nucleotides, or inserting one or more new nucleotides into one sequence, or deleting one or more nucleotides from one sequence.

[0121] Analignment score, representing theextentof thechanges thatwouldbe required tobemade toachieveanexact alignment, is generated,wherein this score and / or other associateddatamaybe stored in thegiven cells of the array.Each cell of thearraycorresponds to thepossibility that thenucleotideat its positionon the readaxis aligns to thenucleotideat its position on the reference axis, and the score generated for each cell represents the partial alignment terminating with the cell’s positions in the read and the reference window. The highest score generated in any cell represents the best overall alignment of the read to the reference window. In various instances, the alignment may be global, where the entire read 22 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 must be aligned to some portion of the reference window, such as using a Needleman-Wunsch or similar algorithm; or in other instances, the alignment may be local, where only a portion of the read may be aligned to a portion of the reference window, such as by using a Smith-Waterman or similar algorithm.

[0122] Accordingly, in various instances, analignment functionmaybeperformed, suchason thedataobtained from the mapping module. Hence, in various instances, an alignment function may form a module, such as an alignment module, that may form part of a system, e.g., a pipeline, that is used, such as in addition with a mapping module, in a process for determining the actual entire genomic sequence, or a portion thereof, of an individual. For instance, the output returned from theperformanceof themapping function, suchas fromamappingmodule, e.g., the list of possibilities as towhereone or more or all of the reads maps to one or more positions in one or more reference genomes, may be employed by the alignment function so as to determine the actual sequence alignment of the subject’s sequenced DNA.

[0123] Such an alignment function may at times be useful because, as described above, often times, for a variety of different reasons, the sequenced reads do not alwaysmatch exactly to the reference genome. For instance, theremay be anSNP (single nucleotide polymorphism) in one ormore of the reads, e.g., a substitution of one nucleotide for another at a single position; there may be an "indel," insertion or deletion of one or more bases along one or more of the read sequences, which insertion or deletion is not present in the reference genome; and / or there may be a sequencing error (e.g., errors in sample prep and / or sequencer read and / or sequencer output, etc.) causing one or more of these apparent variations. Accordingly, when a read varies from the reference, such as by an SNP or Indel, this may be because the reference differs from the trueDNAsequence sampled, or because the read differs from the trueDNAsequence sampled. The problem is to figure out how to correctly align the reads to the reference genome given the fact that in all likelihood the two sequences are going to vary from one another in a multiplicity of different ways.

[0124] As indicated, typically, an algorithm is used to perform such an alignment function. For example, a Smith- Waterman and / or a Needleman-Wunsch alignment algorithm may be employed to align two or more sequences against one another. In this instance, they may be employed in a manner so as to determine the probabilities that for any given position where the read maps to the reference genome that the mapping is in fact the position from where the read originated. Typically, these algorithms are configured so as to be performed by software, however, in various instances, such as herein presented, one or more of these algorithms can be configured so as to be executed in hardware, as described in greater detail herein below.

[0125] In particular, the alignment function operates, at least in part, to align one or more, e.g., all, of the reads to the reference genome despite the presence of one or more portions of mismatches, e.g., SNPs, insertions, deletions, structural artifacts, etc. so as to determinewhere the reads are likely to fit in the genome correctly. For instance, the one or more reads are compared against the reference genome, and the best possible fit for the read against the genome is determined,whileaccounting for substitutionsand / or Indelsand / or structural variants.However, tobetter determinewhich of themodified versions of the read best fits against the reference genome, the proposed changesmust be accounted for, and as such a scoring function may also be performed.

[0126] In view of the above, there are, therefore, at least two goals that may be achieved from performing an alignment function. One is a report of the best alignment, including position in the reference genome and a description of what changes are necessary to make the read match the reference segment at that position, and the other is the alignment quality score. For instance, in various instances, the output from the alignment module may be a Compact Idiosyncratic Gapped Alignment Report, e.g., a CIGAR string, wherein the CIGAR string output is a report detailing all the changes that were made to the reads so as to achieve their best fit alignment, e.g., detailed alignment instructions indicating how the queryactually alignswith the reference.SuchaCIGARstring readoutmaybeuseful in further stagesofprocessingsoas to better determine that for thegiven subject’s genomicnucleotide sequence, thepredicted variationsas comparedagainst a reference genome are in fact true variations, and not just due to machine, software, or human error.

[0127] One or more of such alignment procedures may be performed by any suitable alignment algorithm, such as a Needleman-Wunsch alignment algorithm and / or a Smith-Waterman alignment algorithm that may have beenmodified to accommodate the functionality herein described. In general both of these algorithms and those like them basically perform, in some instances, ina similarmanner. For instance, as set forth above, thesealignment algorithms typically build the virtual array in a similar manner such that, in various instances, the horizontal top boundary may be configured to represent the genomic reference sequence,whichmaybe laid out across the top rowof the array according to its basepair composition. Likewise, the vertical boundary may be configured to represent the sequenced and mapped query sequences that have been positioned in order, downwards along the first column, such that their nucleotide sequence order is generally matched to the nucleotide sequence of the reference to which they mapped. The intervening cells may then be populated with scores as to the probability that the relevant base of the query at a given position, is positioned at that location relative to the reference. In performing this function, a swath may be moved diagonally across the matrix populating scores within the intervening cells and the probability for each base of the query being in the indicated position may be determined.

[0128] With respect to a Needleman-Wunsch alignment function, which generates optimal global (or semi-global) alignments, aligning the entire read sequence to some segment of the reference genome, the wave front steeringmay be 23 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 configured such that it typically sweeps all theway from the top edge of the alignmentmatrix to the bottomedge.When the wave front sweep is complete, themaximumscoreon thebottomedgeof thealignmentmatrix (corresponding to theendof the read) is selected, and the alignment is back-traced to a cell on the top edge of the matrix (corresponding to the beginning of the read). In various of the instances disclosed herein, the reads can be any length long, can be any size, and there neednot beextensive readparameters as to how thealignment is performed, e.g., in various instances, the read can be as long as a chromosome. In such an instance, however, the memory size and chromosome length may be limiting factor.

[0129] With respect to a Smith-Waterman algorithm, which generates optimal local alignments, aligning the entire read sequence or part of the read sequence to some segment of the reference genome, this algorithm may be configured for finding thebest scoringpossiblebasedona full or partial alignmentof the read.Hence, in various instances, thewave front- scored band may not extend to the top and / or bottom edges of the alignment matrix, such as if a very long read had only seeds in itsmiddlemapping to the referencegenome, but commonly thewave frontmay still score from top tobottomof the matrix. Local alignment is typically achieved by two adjustments. First, alignment scores are never allowed to fall below zero (or some other floor), and if a cell score otherwise calculated would be negative, a zero score is substituted, representing the start of a new alignment. Second, the maximum alignment score produced in any cell in the matrix, not necessarily along the bottom edge, is used as the terminus of the alignment. The alignment is backtraced from this maximum score up and left through the matrix to a zero score, which is used as the start position of the local alignment, even if it is not on the top row of the matrix.

[0130] In view of the above, there are several different possible pathways through the virtual array. In various embodiments, the wave front starts from the upper left corner of the virtual array, and moves downwards towards identifiers of the maximum score. For instance, the results of all possible aligns can be gathered, processed, correlated, and scored to determine themaximumscore.When the endof a boundary or the endof the array has been reached and / or a computation leading to the highest score for all of the processed cells is determined (e.g., the overall highest score identified) then a backtracemay be performed so as to find the pathway that was taken to achieve that highest score. For example, a pathway that leads to a predicted maximum score may be identified, and once identified an audit may be performed so as to determine how thatmaximumscorewas derived, for instance, bymoving backwards following the best scorealignmentarrows retracing thepathway that led toachieving the identifiedmaximumscore, suchascalculatedby the wave front scoring cells.

[0131] Once it has been determined where each read is mapped, and further determined where each read is aligned, e.g., each relevant read has been given a position and a quality score reflecting the probability that the position is the correct alignment, such that the nucleotide sequence for the subject’s DNA is known, then the order of the various reads and / or genomic nucleic acid sequence of the subject may be verified, such as by performing a back trace functionmoving backwards up through the array so as to determine the identity of every nucleic acid in its proper order in the sample genomic sequence. Consequently, in some aspects, the present disclosure is directed to a backtrace function, such as is part of an alignmentmodule that performs both an alignment and a back trace function, such as amodule thatmay be part of a pipeline of modules, such as a pipeline that is directed at taking raw sequence read data, such as form a genomic sample form an individual, and mapping and / or aligning that data, which data may then be sorted.

[0132] In the case of affine gap scoring, scoring vector information may be extended, e.g. to 4 bits per scored cell. In addition to the e.g., 2-bit score-choice direction indicator, two 1-bit flags may be added, a vertical extend flag, and a horizontal extend flag. According to the methods of affine gap scoring extensions to Smith-Waterman or Needleman- Wunsch or similar alignment algorithms, for each cell, in addition to the primary alignment score representing the best- scoring alignment terminating in that cell, a ’vertical score’ should begenerated, corresponding to themaximumalignment score reaching that cell with a final vertical step, and a ’horizontal score’ should be generated, corresponding to the maximum alignment score reaching that cell with a final horizontal step; and when computing any of the three scores, a vertical step into the cellmay be computed either using the primary score from the cell aboveminus a gap-openpenalty, or using the vertical score from the cell aboveminus agap-extendpenalty, whichever is greater; andahorizontal step into the cell may be computed either using the primary score from the cell to the left minus a gap-open penalty, or using the horizontal score from the cell to the leftminus a gap-extendpenalty, whichever is greater. In caseswhere the vertical score minus agapextendpenalty is selected, the vertical extendflag in the scoring vector should be set, e.g., ’1’, andotherwise it should be unset, e.g., ’0’.

[0133] In cases when the horizontal score minus a gap extend penalty is selected, the horizontal extend flag in the scoring vector should be set, e.g. ’1’, and otherwise it should be unset, e.g. ’0’. During backtrace for affinegap scoring, any time backtrace takes a vertical step upward from a given cell, if that cell’s scoring vector’s vertical extend flag is set, the following backtrace step must also be vertical, regardless of the scoring vector for the cell above. Likewise, any time backtrace takes a horizontal step leftward from a given cell, if that cell’s scoring vector’s horizontal extend flag is set, the following backtrace stepmust also be horizontal, regardless of the scoring vector for the cell to the left. Accordingly, such a tableof scoring vectors, e.g. 129bits per row for 64cells using linear gapscoring, or 257bits per row for 64cells usingaffine gapscoring,with somenumberNRof rows, is adequate to support backtraceafter concludingalignment scoringwhere the 24 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 scoring wavefront took NR steps or fewer.

[0134] Hence, a method is given for performing incremental backtrace from partial alignment information, e.g., comprising partial scoring vector information for alignment matrix cells scored so far. From a currently completed alignment boundary, e.g., a particular scored wave front position, backtrace is initiated from all cell positions on the boundary. Such backtrace from all boundary cells may be performed sequentially, or advantageously, especially in a hardware implementation, all thebacktracesmaybeperformed together. It is not necessary toextract alignment notations, e.g., CIGAR strings, from thesemultiple backtraces; only to determinewhat alignmentmatrix positions they pass through during thebacktrace. In an implementation of simultaneous backtrace fromascoring boundary, a number of 1-bit registers may be utilized, corresponding to the number of alignment cells, initialized e.g., all to ’1’s, representing whether any of the backtraces pass through a corresponding position. For each step of simultaneous backtrace, scoring vectors correspond- ing to all the current ’1’s in these registers, e.g. from one row of the scoring vector table, can be examined, to determine a next backtrace step corresponding to each ’1’ in the registers, leading to a following position for each ’1’ in the registers, for the next simultaneous backtrace step.

[0135] Importantly, it is easily possible formultiple ’1’s in the registers tomerge into commonpositions, corresponding to multiple of the simultaneous backtraces merging together onto common backtrace paths. Once two or more of the simultaneous backtraces merge together, they remain merged indefinitely, because henceforth they will utilize scoring vector information from the same cell. It has been observed, empirically and for theoretical reasons, that with high probability, all of the simultaneous backtraces merge into a singular backtrace path, in a relatively small number of backtracesteps,whiche.g.maybeasmallmultiple, e.g. 8, times thenumberof scoringcells in thewavefront. Forexample, with a 64-cell wavefront, with high probability, all backtraces from a given wavefront boundary merge into a single backtrace path within 512 backtrace steps. Alternatively, it is also possible, and not uncommon, for all backtraces to terminate within the number, e.g. 512, of backtrace steps.

[0136] Accordingly, the multiple simultaneous backtraces may be performed from a scoring boundary, e.g. a scored wavefront position, far enough back that they all either terminate or merge into a single backtrace path, e.g. in 512 backtrace steps or fewer. If they all merge together into a singular backtrace path, then from the location in the scoring matrix where they merge, or any distance further back along the singular backtrace path, an incremental backtrace from partial alignment information is possible. Further backtrace from the merge point, or any distance further back, is commenced, by normal singular backtrace methods, including recording the corresponding alignment notation, e.g., a partial CIGAR string. This incremental backtrace, and e.g., partial CIGAR string, must be part of any possible final backtrace, and e.g., full CIGAR string, that would result after alignment completes, unless such final backtrace would terminate before reaching the scoring boundary where simultaneous backtrace began, because if it reaches the scoring boundary, it must follow one of the simultaneous backtrace paths, and merge into the singular backtrace path, now incrementally extracted.

[0137] Therefore, all scoring vectors for thematrix regions corresponding to the incrementally extractedbacktrace, e.g., in all table rows for wave front positions preceding the start of the extracted singular backtrace, may be safely discarded. When the final backtrace is performed from amaximum scoring cell, if it terminates before reaching the scoring boundary (or alternatively, if it terminates before reaching the start of the extracted singular backtrace), the incremental alignment notation, e.g. partial CIGAR string, may be discarded. If the final backtrace continues to the start of the extracted singular backtrace, its alignment notation, e.g., CIGAR string, may then be grafted onto the incremental alignment notation, e.g., partial CIGAR string. Furthermore, in a very long alignment, the process of performing a simultaneous backtrace from a scoringboundary,e.g., scoredwave front position,until all backtraces terminateormerge, followedbyasingular backtrace with alignment notation extraction, may be repeated multiple times, from various successive scoring boundaries. The incremental alignment notation, e.g. partial CIGAR string, from each successive incremental backtrace may then be grafted onto the accumulated previous alignment notations, unless the newsimultaneous backtrace or singular backtrace terminates early, inwhich caseaccumulatedpreviousalignment notationsmaybediscarded. Theeventual final backtrace likewise grafts its alignment notation onto the most recent accumulated alignment notations, for a complete backtrace description, e.g., CIGAR string.

[0138] Accordingly, in this manner, thememory to store scoring vectors may be kept bounded, assuming simultaneous backtraces always merge together in a bounded number of steps, e.g. 512 steps. In rare cases where simultaneous backtraces fail tomerge or terminate in the bounded number of steps, various exceptional actionsmay be taken, including failing the current alignment, or repeating it with a higher bound or with no bound, perhaps by a different or traditional method, such as storing all scoring vectors for the complete alignment, such as in external DRAM. In a variation, it may be reasonable to fail such an alignment, because it is extremely rare, and even rarer that such a failed alignment would have been a best-scoring alignment to be used in alignment reporting.

[0139] In various instances, the devices, systems, and theirmethods of use of the present disclosuremaybe configured for performingoneormoreof a full-readgapless and / or gappedalignments thatmay thenbescored soas todetermine the appropriate alignment for the reads in the dataset. For instance, in various instances, a gapless alignment proceduremay be performed on data to be processed, which gapless alignment procedure may then be followed by one or more of a 25 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 gapped alignment, and / or by a selective Smith-Waterman alignment procedure. For example, in a first step, a gapless alignment chain may be generated. As described herein, such gapless alignment functions may be performed quickly, such as without the need for accounting for gaps, which after a first step of performing a gapless alignment, may then be followed by then performing a gapped alignment.

[0140] For instance, analignment functionmaybeperformed in order to determinehowanygivennucleotide sequence, e.g., read, aligns to a reference sequencewithout the need for inserting gaps in one ormore of the reads and / or reference. An important part of performing such an alignment function is determining where and how there are mismatches in the sequence in question versus the sequence of the reference genome. However, because of the great homology within the human genome, in theory, any given nucleotide sequence is going to largelymatch a representative reference sequence. Where there are mismatches, these will likely be due to a single nucleotide polymorphism, which is relatively easy to detect, or they will be due to an insertion or deletion in the sequences in question, which are muchmore difficult to detect.

[0141] Consequently, in performing an alignment function, themajority of the time, the sequence in question is going to match the reference sequence, and where there is a mismatch due to an SNP, this will easily be determined. Hence, a relatively large amount of processing power is not required to perform such analysis. Difficulties arise, however, where there are insertions or deletions in the sequence in question with respect to the reference sequence, because such insertionsanddeletionsamount togaps in thealignment.Suchgaps requireamoreextensiveandcomplicatedprocessing platform so as to determine the correct alignment. Nevertheless, because there will only be a small percentage of indels, only a relatively smaller percentage of gapped alignment protocols need be performed as compared to the millions of gaplessalignments performed.Hence, only a small percentageof all of thegaplessalignment functions result in aneed for further processing due to the presence of an indel in the sequence, and therefore will need a gapped alignment.

[0142] When an indel is indicated in a gapless alignment procedure, only those sequences get passed on to an alignment engine for further processing, such as an alignment engine configured for performing an advanced alignment function, such as a Smith Waterman alignment (SWA). Thus, because either a gapless or a gapped alignment is to be performed, the devices and systems disclosed herein are a much more efficient use of resources. More particularly, in certain embodiments, both a gapless and a gapped alignment may be performed on a given selection of sequences, e.g., one right after the other, then the results are compared for each sequence, and the best result is chosen. Such an arrangementmay be implemented, for instance,where an enhancement in accuracy is desired, and an increased amount of time and resources for performing the required processing is acceptable.

[0143] Particularly, in various instances, a first alignment step may be performed without engaging a processing intensive Smith Waterman function. Hence, a plurality of gapless alignments may be performed in a less resource intensive, less time-consuming manner, and because less resources are needed less space need be dedicated for such processing on the chip. Thus, more processing may be performed, using less processing elements, requiring less time, therefore, more alignments can be done, and better accuracy can be achieved. More particularly, less chip resource- implementations for performingSmithWatermanalignmentsneedbededicatedusing less chiparea, as it doesnot require asmuch chip area for the processing elements required to perform gapless alignments as it does for performing a gapped alignment. As the chip resource requirements go down, themore processing can be performed in a shorter period of time, and with the more processing that can be performed, the better the accuracy can be achieved.

[0144] The output from the alignment module is a SAM (Text) or BAM (e.g., binary version of a SAM) file along with a mapping quality score (MAPA), which quality score reflects the confidence that the predicted and aligned location of the read to the reference is actually where the read is derived. Accordingly, once it has been determined where each read is mapped, and further determined where each read is aligned, e.g., each relevant read has been given a position and a quality score reflecting the probability that the position is the correct alignment, such that the nucleotide sequence for the subject’sDNA is knownaswell ashow thesubject’sDNAdiffers from that of the reference (e.g., theCIGARstringhasbeen determined), then the various reads representing the genomic nucleic acid sequence of the subject may be sorted by chromosome location, so that the exact location of the read on the chromosomes may be determined. Consequently, in some aspects, the present disclosure is directed to a sorting function, such as may be performed by a sorting module, which sortingmodulemay be part of a pipeline of modules, such as a pipeline that is directed at taking raw sequence read data, such as form a genomic sample form an individual, and mapping and / or aligning that data, which data may then be sorted.

[0145] Moreparticularly, once the readshavebeenassignedaposition, suchas relative to the referencegenome,which may include identifying to which chromosome the read belongs and / or its offset from the beginning of that chromosome, the readsmay be sorted by position. Sortingmay be useful, such as in downstreamanalyses, whereby all of the reads that overlap agivenposition in the genomemaybe formed into apile up soas to beadjacent to oneanother, suchasafter being processed through the sorting module, whereby it can be readily determined if the majority of the reads agree with the reference value or not. Hence, where the majority of reads do not agree with the reference value a variant call can be flagged. Sorting, therefore, may involve one ormore of sorting the reads that align to the relatively same position, such as the same chromosome position, so as to produce a pileup, such that all the reads that cover the same location are physically grouped together; andmay further involve analyzing the reads of the pileup to determine where the readsmay 26 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 indicate an actual variant in the genome, as compared to the reference genome, which variant may be distinguishable, such as by the consensus of the pileup, from an error, such as a machine read error or error an error in the sequencing methods which may be exhibited by a small minority of the reads.

[0146] Once the data has beenobtained there are oneormore othermodules thatmaybe run so as to clean up the data. For instance, onemodule that may be included, for example, in a sequence analysis pipeline, such as for determining the genomic sequence of an individual, may be a local realignment module. For example, it is often difficult to determine insertions and deletions that occur at the end of the read. This is because the Smith-Waterman or equivalent alignment process lacks enough context beyond the indel to allow the scoring to detect its presence. Consequently, the actual indel may be reported as one ormore SNPs. In such an instance, the accuracy of the predicted location for any given readmay be enhanced by performing a local realignment on the mapped and / or aligned and / or sorted read data.

[0147] In such instances, pileupsmay be used to help clarify the proper alignment, such aswhere a position in question isat theendof anygiven read, that sameposition is likely tobeat themiddleof someother read in thepileup.Accordingly, in performing a local realignment the various reads in a pileupmay be analyzed so as to determine if someof the reads in the pile up indicate that therewasan insertionor adeletionat a givenpositionwhereanother readdoesnot include the indel, or rather includes a substitution, at that position, then the indel may be inserted, such as into the reference, where it is not present, and the reads in the local pileup that overlap that region may be realigned to see if collectively a better score is achieved then when the insertion and / or deletion was not there. If there is an improvement, the whole set of reads in the pileupmay be reviewed and if the score of the overall set has improved then it is clear tomake the call that there really was an indel at that position. In amanner suchas this, the fact that there is not enough context tomoreaccurately align a readat the end of a chromosome, for any individual read,may be compensated for. Hence, when performing a local realignment, oneormorepileupswhereoneormore indelsmaybepositionedare examined, and it is determined if by addingan indel at any given position the overall alignment score may be enhanced.

[0148] Another module that may be included, for example, in a sequence analysis pipeline, such as for determining the genomic sequenceof an individual,maybeaduplicatemarkingmodule.For instance, aduplicatemarking functionmaybe performed so as to compensate for chemistry errors that may occur during the sequencing phase. For example, as described above, during some sequencing procedures nucleic acid sequences are attached to beads and built up from there using labeled nucleotide bases. Ideally there will be only one read per bead. However, sometimes multiple reads become attached to a single bead and this results in an excessive number of copies of the attached read. This phenomenon is known as read duplication.

[0149] Afteranalignment isperformedand the resultsobtained, and / orasorting function, local realignment, and / or ade- duplication is performed, a variant call functionmay be employed on the resultant data. For instance, a typical variant call function or parts thereofmay be configured so as to be implemented in a software and / or hardwired configuration, such as on an integrated circuit. Particularly, variant calling is a process that involves positioning all the reads that align to a given locationon the reference intogroupingssuch that all overlapping regions fromall the variousaligned reads forma "pile up." Then the pileup of reads covering a given region of the reference genome are analyzed to determine what the most likely actual content of the sampled individual’s DNA / RNA iswithin that region. This is then repeated, stepwise, for every region of the genome. The determined content generates a list of differences termed "variations" or "variants" from the reference genome, each with an associated confidence level along with other metadata.

[0150] Themost common variants are single nucleotide polymorphisms (SNPs), in which a single base differs from the reference. SNPs occur at about 1 in 1000 positions in a human genome. Next most common are insertions (into the reference) and deletions (from the reference), or "indels" collectively. These aremore common at shorter lengths, but can be of any length. Additional complications arise, however, because the collection of sequenced segments ("reads") is random, some regions will have deeper coverage than others. There are also more complex variants that include multi- base substitutions, and combinations of indels and substitutions that can be thought of as length-altering substitutions. Standard software based variant callers have difficulty identifying all of these, and with various limits on variant lengths. More specialized variant callers in both software and / or hardware are needed to identify longer variations, and many varieties of exotic "structural variants" involving large alterations of the chromosomes.

[0151] However, variant calling is adifficult procedure to implement in software, andworldsofmagnitudemoredifficult to deploy in hardware. In order to account for and / or detect these types of errors, typical variant callers may perform one or more of the following tasks. For instance, theymay come up with a set of hypothesis genotypes (content of the one or two chromosomes at a locus), use Bayesian calculations to estimate the posterior probability that each genotype is the truth given the observed evidence, and report the most likely genotype along with its confidence level. As such variant callers may be simple or complex. Simpler variant callers look only at the columnof bases in the aligned read pileup at the precise position of a call being made. More advanced variant callers are "haplotype based callers", which may be configured to take into account context, such as in a window, around the call being made.

[0152] A "haplotype" is particular DNA content (nucleotide sequence, list of variants, etc.) in a single common "strand", e.g. one of two diploid strands in a region, and a haplotype based caller considers the Bayesian implications of which differences are linked by appearing in the same read. Accordingly, a variant call protocol, as proposed herein, may 27 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 implement one or more improved functions such as those performed in a Genome Analysis Tool Kit (GATK) haplotype caller and / or usingaHiddenMarkovModel (HMM) tool and / or aDeBruijnGraph function, suchaswhereoneormore these functions typically employed by a GATK haplotype caller, and / or a HMM tool, and / or a De Bruijn Graph function may be implemented in software and / or in hardware.

[0153] More particularly, as implemented herein, various different variant call operationsmay be configured so as to be performed in software or hardware, andmay include one or more of the following steps. For instance, variant call function may includeanactive region identification, suchas for identifying placeswheremultiple readsdisagreewith the reference, and for generating a window around the identified active region, so that only these regions may be selected for further processing. Additionally, localizedhaplotypeassemblymay takeplace, suchaswhere, for eachgivenactive region, all the overlapping reads may be assembled into a "De Bruijn graph" (DBG) matrix. From this DBG, various paths through the matrix may be extracted, where each path constitutes a candidate haplotype, e.g., hypotheses, for what the true DNA sequence may be on at least one strand. Further, haplotype alignment may take place, such as where each extracted haplotype candidate may be aligned, e.g., Smith-Waterman aligned, back to the reference genome, so as to determine what variation(s) from the reference it implies.Furthermore,a read likelihoodcalculationmaybeperformed, suchaswhere each readmay be tested against each haplotype, or hypothesis, to estimate a probability of observing the read assuming the haplotype was the true original DNA sampled.

[0154] With respect to these processes, the read likelihood calculation will typically be the most resource intensive and time-consuming operation to be performed, often requiring a pair HMM evaluation. Additionally, the constructing of De Bruijn graphs for each pileup of reads, with associated operations of identifying locally and globally unique K-mers, as described below may also be resource intensive and / or time consuming. Accordingly, in various embodiments, one or more of the various calculations involved in performing one or more of these steps may be configured so as to be implemented in optimized software fashion or hardware, such as for being performed in an accelerated manner by an integrated circuit, as herein described.

[0155] As indicated above, in various embodiments, a Haplotype Caller of the disclosure, implemented in software and / or in hardware or a combination thereof may be configured to include one or more of the following operations: Active Region Identification, Localized Haplotype Assembly, Haplotype Alignment, Read Likelihood Calculation, and / or Geno- typing.For instance, thedevices, systems, and / ormethodsof thedisclosuremaybeconfigured toperformoneormoreof a mapping, aligning, and / or a sorting operation on data obtained from a subject’s sequenced DNA / RNA to generate mapped, aligned, and / or sorted results data. This results data may then be cleaned up, such as by performing a de duplication operation on it and / or that data may be communicated to one or more dedicated haplotype caller processing engines for performing a variant call operation, including one or more of the aforementioned steps, on that results data so as to generate a variant call file with respect thereto. Hence, all the reads that have been sequenced and / or beenmapped and / or aligned to particular positions in the reference genomemay be subjected to further processing so as to determine how the determined sequence differs from a reference sequence at any given point in the reference genome.

[0156] Accordingly, in various embodiments, a device, system, and / or method of its use, as herein disclosed, may include avariant or haplotype caller system that is implemented in a software and / or hardwired configuration to performan active region identification operation on the obtained results data. Active region identification involves identifying and determining places where multiple reads, e.g., in a pile up of reads, disagree with a reference, and further involves generating one ormore windows around the disagreements ("active regions") such that the region within the windowmay beselected for furtherprocessing.Forexample, duringamappingand / oraligningstep, identified readsaremappedand / or aligned to the regions in the reference genome where they are expected to have originated in the subject’s genetic sequence.

[0157] However, as the sequencing is performed in such amanner so as to create an oversampling of sequenced reads for any given region of the genome, at any given position in the reference sequencemay be seen a pile up of any and / all of the sequenced reads that line up and align with that region. All of these reads that align and / or overlap in a given region or pile up position may be input into the variant caller system. Hence, for any given read being analyzed, the read may be compared to the referenceat its suspected regionof overlap, and that readmaybecompared to the reference todetermine if it shows any difference in its sequence from the known sequence of the reference. If the read lines up to the reference, without any insertions or deletions and all the bases are the same, then the alignment is determined to be good.

[0158] Hence, for any givenmappedand / or aligned read, the readmay havebases that are different from the reference, e.g., the readmay include one or more SNPs, creating a position where a base is mismatched; and / or the readmay have one or more of an insertion and / or deletion, e.g., creating a gap in the alignment. Accordingly, in any of these instances, there will be one ormoremismatches that need to be accounted for by further processing. Nevertheless, to save time and increaseefficiency, such further processingshouldbe limited to those instanceswhereaperceivedmismatch isnon-trivial, e.g., a non-noise difference.

[0159] In determining the significance of a mismatch, places where multiple reads in a pile up disagree from the reference may be identified as an active region, a window around the active region may then be used to select a locus of disagreement that may then be subjected to further processing. The disagreement, however, should be non-trivial. This 28 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 maybedetermined inmanyways, for instance, thenon-referenceprobabilitymaybecalculated for each locus in question, such as by analyzing basematch vsmismatch quality scores, such as above a given threshold deemed to be a sufficiently significant amount of indication from those reads that disagree with the reference in a significant way.

[0160] For instance, if 30 of themapped and / or aligned reads all line up and / or overlap so as to form a pile up at a given position in the reference, e.g., an active region, and only 1 or 2 out of the 30 reads disagrees with the reference, then the minimal threshold for further processing may be deemed to not have been met, and the non-agreeing read(s) can be disregarded in view of the 28 or 29 reads that do agree. However, if 3 or 4, or 5, or 10, or more of the reads in the pile up disagree, then thedisagreementmaybe statistically significant enough towarrant further processing, andanactive region around the identified region(s) of difference might be determined. In such an instance, an active region window ascertaining the bases surrounding that difference may be taken to give enhanced context to the region surrounding the difference, and additional processing steps, such as performing a Gaussian distribution and sum of non-reference probabilities distributed across neighboring positions,may be taken to further investigate and process that region to figure out if and active region should be declared and if so what variances from the reference actually are present within that region if any. Therefore, the determining of an active region identifies those regions where extra processing may be needed to clearly determine if a true variance or a read error has occurred.

[0161] Particularly, because in many instances it is not desirable to subject every region in a pile up of sequences to further processing, an active region can be identified whereby it is only those regions where extra processing may be needed to clearly determine if a true variance or a read error has occurred that may be determined as needing of further processing. And, as indicated above, it may be the size of the supposed variance that determines the size of thewindowof the active region. For instance, in various instances, the boundsof the activewindowmayvary from1or 2or about 10or 20 or even about 25 or about 50 to about 200 or about 300, or about 500 or about 1000 bases long or more, where it is only within the bounds of the active window that further processing is taking place. Of course, the size of the active window can be any suitable length so long as it provides the context to determine the statistical importance of a difference.

[0162] Hence, if there are only one or two isolated differences, then the active window may only need to cover one or more to a few dozen bases in the active region so as to have enough context tomake a statistical call that an actual variant is present. However, if there is a cluster or a bunch of differences, or if there are indels present for which more context is desired, then thewindowmaybe configured soas to be larger. In either instance, itmay be desirable to analyze any andall the differences that might occur in clusters, so as to analyze them all in one or more active regions, because to do so can provide supporting informationabout each individual differenceandwill saveprocessing timebydecreasing thenumberof activewindowsengaged. In various instances, theactive regionboundariesmaybedeterminedbyactiveprobabilities that pass a given threshold, such as about 0.00001 or about 0.00001 or about 0.0001 or less to about 0.002 or about 0.02 or about 0.2 or more. And if the active region is longer than a given threshold, e.g., about 300 - 500 bases or 1000 bases or more, then the regioncanbebrokenup into sub-regions, suchasby sub-regionsdefinedby the locuswith the lowest active probability score.

[0163] In various instances, after an active region is identified, a localized haplotype assembly procedure may be performed. For instance, in each active region, all the piled up and / or overlapping reads may be assembled into a "De Bruijn Graph" (DBG). A DBGmay be a directed graph based on all the reads that overlapped the selected active region, which active region may be about 200 or about 300 to about 400 or about 500 bases long or more, within which active region the presence and / or identity of variants are to be determined. In various instances, as indicated above, the active region can be extended, e.g., by including another about 100 or about 200 or more bases in each direction of the locus in question so as to generate an extended active region, such as where additional context surrounding a difference may be desired.Accordingly, it is from theactive regionwindow,extendedor not, that all of the reads that haveportions that overlap theactive region arepiled up, e.g., to produceapileup, the overlappingportions are identified, and the read sequencesare threaded into the haplotype caller system and are thereby assembled together in the form of a De Bruin graph, much like the pieces of a puzzle.

[0164] For any given active window there will be reads that form a pile up such that en masse the pile up will include a sequence pathway throughwhich the overlapping regions of the various overlapping reads in the pile up covers the entire sequence within the active window. Hence, at any given locus in the active region, there will be a plurality of reads overlapping that locus, albeit any given read may not extend the entire active region. The result of this is that various regionsof various readswithinapileupareemployedby theDBG indeterminingwhethera variant actually is present ornot for any given locus in the sequencewithin the active region. As it iswithin the activewindow that this determination is being made, it is those portions of any given readwithin the borders of the activewindow that are considered, and those portions that are outside of the active window may be discarded.

[0165] As indicated, it is thosesectionsof the reads that overlap the referencewithin theactive region that are fed into the DBG system. The DBG system then assembles the reads like a puzzle into a graph, and then for each position in the sequence, it is determined based on the collection of overlapping reads for that position, whether there is a match or a mismatch for any given, and if there is a mismatch, what the probability of that mismatch is. For instance, where there are discreteplaceswheresegmentsof the reads in thepileupoverlapeachother, theymaybealigned tooneanotherbasedon 29 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 their areas of matching, and from stringing or stitching the matching reads together, as determined by their points of matching, it can be established for each position within that segment, whether and to what extent the reads at any given position match or mismatch each other. Hence, if two or more reads being compiled line up and match each other identically for a while, a graph having a single string will result; however, when the two or more reads come to a point of difference, a branch in the graph will form, and two or more divergent strings will result, until matching between the two or more reads resumes.

[0166] Hence, the pathways through the graph are often not a straight line. For instance, where the k-mers of a read varies from thek-mersof the referenceand / or thek-mers fromoneormoreoverlapping reads, e.g., in thepileup, a "bubble" will be formed in the graph at the point of difference resulting in two divergent strings that will continue along two different path linesuntilmatchingbetween the twosequences resumes.Eachvertexmaybegivenaweightedscore identifyinghow many times the respective k-mers overlap in all of the reads in the pileup.Particularly, each pathwayextending through the generated graph from one side to the other may be given a count. And where the same k-mers are generated from a multiplicity of reads, e.g., where each k-mer has the same sequence pattern, they may be accounted for in the graph by increasing thecount for thatpathwaywhere thek-meroverlapsanalreadyexistingk-merpathway.Hence,where thesame k-mer is generated fromamultiplicity of overlapping readshaving the samesequence, the pattern of the pathwaybetween thegraphwill be repeatedoverandoveragainand thecount for traversing thispathway through thegraphwill be increased incrementally in correspondence therewith. In such an instance, the pattern is only recorded for the first instance of the k- mer, and the count is incrementally increased for each k-mer that repeats that pattern. In thismode the various reads in the pile up can be harvested to determine what variations occur and where.

[0167] In a manner such as this, a graph matrix may be formed by taking all possible N base k-mers, e.g., 10 base k- mers, which can be generated from each given read by sequentially walking the length of the read in ten base segments, where the beginning of each new ten base segment is offset by one base from the last generated 10 base segment. This proceduremay then be repeated by doing the same for every read in the pile upwithin the activewindow. The generated k- mers may then be aligned with one another such that areas of identical matching between the generated k-mers are matched to the areaswhere they overlap, so as to build up a data structure, e.g., graph, thatmay then be scanned and the percentage ofmatching andmismatchingmaybe determined. Particularly, the reference and any previously processed k- mers aligned therewith may be scanned with respect to the next generated k-mer to determine if the instant generated k- mer matches and / or overlaps any portion of a previously generated k-mer, and where it is found to match the instant generated k-mer can then be inserted into the graph at the appropriate position.

[0168] Once built, the graph can be scanned and itmay be determined based on thismatchingwhether any givenSNPs and / or indels in the readswith respect to the referenceare likely tobeanactual variation in the subject’s genetic codeor the result of aprocessingorothererror. For instance, if all or asignificant portionof thek-mers, of all or a significant portionofall of the reads, in a given region include the same SNP and / or indel mismatch, but differ from the reference in the same manner, then it may be determined that there is an actually SNP and / or indel variation in the subject’s genome as compared to the reference genome. However, if only a limited number of k-mers from a limited number of reads evidence the artifact, it is likely to be caused bymachine and / or processing and / or other error and not indicative of a true variation at the position in question.

[0169] As indicated,where there is a suspected variance, a bubblewill be formedwithin thegraph.Specifically,whereall of thek-merswithinall of agiven regionof readsallmatch the reference, theywill lineup in suchamanneras to forma linear graph. However, where there is a difference between the bases at a given locus, at that locus of difference that graph will branch. This branchingmay be at any position within the k-mer, and consequently at that point of difference the 10 base k- mer, including that difference,will diverge from the rest of the k-mers in thegraph. In suchan instance, a newnode, forming a different pathway through the graph will be formed.

[0170] Hence, where everything may have been agreeing, e.g., the sequence in the given new k-mer being graphed is matching the sequence towhich it aligns in the graph, up to the point of difference the pathway for that k-merwill match the pathway for the graph generally and will be linear, but post the point of difference, a new pathway through the graph will emerge to accommodate the difference represented in the sequence of the newly graphed k-mer. This divergence being represented by a new nodewithin the graph. In such an instance, any new k-mers to be added to the graph thatmatch the newly divergent pathwaywill increase the count at that node. Hence, for every read that supports the arc, the count will be increased incrementally.

[0171] In various of such instances, the k-mer and / or the read it represents will once again start matching, e.g., after the point of divergence, such that there is now a point of convergence where the k-mer begins matching the main pathway through the graph represented by the k-mers of the reference sequence. For instance, naturally after a while the read(s) that support the branched node should rejoin the graph over time. Thus, over time, the k-mers for that read will rejoin the main pathway again.More particularly, for anSNPat a given locuswithin a read, the k-mer starting at that SNPwill diverge from the main graph and will stay separate for about 10 nodes, because there are 10 bases per k-mer that overlap that locus ofmismatching between the readand the reference.Hence, for anSNP, at the 11th position, the k-mers covering that locuswithin the readwill rejoin themain pathway as exactmatching is resumed. Consequently, it will take ten shifts for the 30 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 k-mers of a read having an SNP at a given locus to rejoin the main graph represented by the reference sequence.

[0172] As indicatedabove, there is typically onemainpathor lineorbackbone that is the referencepath, andwhere there is a divergence a bubble is formed at a node where there is a difference between a read and the backbone graph. Thus, there are some reads that diverge from the backbone and form a bubble, which divergence may be indicative of the presence of a variant. As the graph is processed, bubbles within bubbles within bubbles may be formed along the reference backbone, so that they are stacked up and a plurality of pathways through the graphmay be created. In such an instance, there may be a main path represented by the reference backbone, one path of a first divergence, and a further path of a second divergence within the first divergence, all within a given window, each pathway through the graph may represent an actual variation or may be an artifact such as caused by sequencing error, and / or PCR error, and / or a processing error, and the like.

[0173] Once such a graph has been produced, it must be determined which pathways through the graph represent actual variations present within the sample genome and which are mere artifacts. Albeit, it is expected that reads containing handling or machine errors will not be supported by the majority of reads in the sample pileup, however, this is not always the case. For instance, errors in PCR processing may typically be the result of a cloning mistake that occurs when preparing the DNA sample, such mistakes tend to result in an insertion and / or a deletion being added to the cloned sequence. Such indel errors may be more consistent among reads, and can wind up with generating multiple reads that have the sameerror from thismistake inPCRcloning.Consequently, a higher count line for suchapoint of divergencemay result because of such errors.

[0174] Hence, once a graph matrix has been formed, with many paths through the graph, the next stage is to traverse and thereby extract all of the paths through the graph, e.g., left to right, e.g., so as to derive one or more candidate haplotypes therefrom. One path will be the reference backbone, but there will be other paths that follow various bubbles along the way. All paths must be traversed and their count tabulated. For instance, if the graph includes a pathway with a two-level bubble in one spot and a three-level bubble in another spot, there will be (2 x 3)6 paths through that graph. So, each of the paths will individually need to be extracted, which extracted paths are termed as candidate haplotypes. Such candidate haplotypes represent theories for what could really be representative of the subject’s actual DNA that was sequenced, and the following processing steps, including one ormore of haplotype alignment, read likelihood calculation, and / or genotyping may be employed to test these theories so as to find out the probabilities that anyone and / or each of these theories is correct. The implementation of a De Bruijn graph reconstruction therefore represents a way to reliably extract a good set of hypotheses to test.

[0175] For instance, in performing a variant call function, as disclosed herein, an active region identification operation may be implemented, such as for identifying places wheremultiple reads in a pile up within a given region disagree with a reference, e.g., a standard or chimeric reference, and for generating a window around the identified active region, so that only these regions may be selected for further processing. Additionally, localized haplotype assembly may take place, such as where, for each given active region, all the overlapping reads in the pile up may be assembled into a "De Bruijn graph" (DBG) matrix. From this DBG, various paths through the matrix may be extracted, where each path constitutes a candidate haplotype, e.g., hypotheses, for what the true DNA sequence may be on at least one strand.

[0176] Further, haplotype alignment may take place, such as where each extracted haplotype candidate may be aligned, e.g., Smith-Waterman aligned, back to the reference genome, so as to determine what variation(s) from the reference it implies. Furthermore, a read likelihood calculationmay beperformed, suchaswhere each readmay be tested against each haplotype, to estimate a probability of observing the read assuming the haplotype was the true original DNA sampled. Finally, a genotyping operation may be implement, and a variant call file produced.

[0177] As indicated above, any or all of these operations may be configured so as to be implemented in an optimized manner in software and / or in hardware, and in various instances, because of the resource intensive and time consuming nature of building aDBGmatrix and extracting candidate haplotypes therefrom, and / or because of the resource intensive and time consuming nature of performing a haplotype alignment and / or a read likelihood calculation, which may include the engagement of an Hidden Markov Model (HMM) evaluation, these operations (e.g., localized haplotype assembly, and / or haplotypealignment, and / or read likelihood calculation) or a portion thereofmaybe configured so as to have oneor more functions of their operation implemented in a hardwired form, such as for being performed in an acceleratedmanner byan integrated circuit asdescribedherein. In various instances, these tasksmaybeconfigured tobe implementedbyone or more quantum circuits such as in a quantum computing device.

[0178] Accordingly, in various instances, thedevices, systems,andmethods forperforming thesamemaybeconfigured so as to perform a haplotype alignment and / or a read likelihood calculation. For instance, as indicated, each extracted haplotype may be aligned, such as Smith-Waterman aligned, back to the reference genome, so as to determine what variation(s) from the reference it implies. In various exemplary instances, scoring may take place, such as in accordance with the followingexemplary scoringparameters: amatch=20.0; amismatch= ‑15.0; agapopen ‑26.0; andagapextend= ‑1.1, other scoring parameters may be used. Accordingly, in this manner, a CIGAR strand may be generated and associatedwith thehaplotype toproduceanassembledhaplotype,whichassembledhaplotypemayeventually beused to identify variants. Accordingly, in a manner such as this, the likelihood of a given read being associated with a given 31 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 haplotypemaybecalculated forall read / haplotypecombinations. In such instances, the likelihoodmaybecalculatedusing a Hidden Markov Model (HMM).

[0179] For instance, thevariousassembledhaplotypesmaybealigned inaccordancewithadynamicprogramingmodel similar to a SWalignment. In such an instance, a virtual matrix may be generated such as where the candidate haplotype, e.g., generated by the DBG, may be positioned on one axis of a virtual array, and the readmay be positioned on the other axis. Thematrix may then be filled out with the scores generated by traversing the extracted paths through the graph and calculating the probabilities that any given path is the true path.

[0180] Hence, in such an instance, a difference in this alignment protocol from a typical SW alignment protocol is that with respect to finding the most likely path through the array, a maximum likelihood calculation may be used, such as a calculation performed by an HMMmodel that is configured to provide the total probability for alignment of the reads to the haplotype. Hence, an actual CIGAR strand alignment, in this instance, need not be produced. Rather all possible alignments are considered and their possibilities are summed. The pair HMM evaluation is resource and time intensive, and thus, implementing its operations within a hardwired configuration within an integrated circuit or via quantum circuits on a quantum computing platform is very advantageous.

[0181] For example, each read may be tested against each candidate haplotype, so as to estimate a probability of observing the read assuming the haplotype is the true representative of the original DNA sampled. In various instances, this calculationmay be performed by evaluating a "pair hiddenMarkovmodel" (HMM), whichmay be configured tomodel the various possible ways the haplotype candidatemight have beenmodified, such as by PCR or sequencing errors, and the like, and a variation introduced into the read observed. In such instances, the HMMevaluationmay employ a dynamic programmingmethod tocalculate the total probability of anyseriesofMarkovstate transitionsarrivingat theobserved read in view of the possibility that any divergence in the read may be the result of an error model. Accordingly, such HMM calculations may be configured to analyze all the possible SNPs and Indels that could have been introduced into one or more of the reads, such as by amplification and / or sequencing artifacts.

[0182] Particularly, paired HMM considers in a virtual matrix all the possible alignments of the read to the reference candidate haplotypes alongwith a probability associatedwith each of them,where all probabilities are added up. The sum of all of the probabilities of all the variants along a given path is added up to get one overarching probability for each read. This process is then performed for every pair, for every haplotype, read pair. For example, if there is a six pile up cluster overlapping a given region, e.g., a region of six haplotype candidates, and if the pile up includes about one hundred reads, 600HMMoperationswill thenneed tobeperformed.Moreparticularly, if thereare6haplotypes then therearegoing tobe6 branches through the path and the probability that each one is the correct pathway that matches the subject’s actual genetic code for that region must be calculated. Consequently, each pathway for all of the reads may be considered, and the probability for each read that you would arrive at this given haplotype is to be calculated.

[0183] The pair Hidden Markov Model is an approximate model for how a true haplotype in the sampled DNA may transform into a possible different detected read. It has been observed that these types of transformations are a combination of SNPs and Indels that have been introduced into the genetic sample set by the PCR process, by one or more of the other sample preparation steps, and / or by an error caused by the sequencing process, and the like. As can be seen with respect to FIG. 2, to account for these types of errors, an underlying 3-state base model may be employed, such as where: (M = alignment match, I = insertion, D = deletion), further where any transition is possible except I <-> D.

[0184] Ascanbeseenwith respect toFIG.2, the3-statebasemodel transitionsarenot ina timesequence, but rather are in a sequence of progression through the candidate haplotype and read sequences, beginning at position 0 in each sequence,where the first base is position 1. A transition toM implies position +1 in both sequences; a transition to I implies position+1 in the readsequenceonly; anda transition toD impliesposition+1 in thehaplotypesequenceonly. Thesame3- state model may be configured to underlie the Smith-Waterman and / or Needleman-Wunsch alignments, as herein described, as well. Accordingly, such a 3-state model, as set forth herein, may be employed in a SWand / or NW process thereby allowing for affine gap (indel) scoring, in which gap opening (entering the I or D state) is assumed to be less likely than gap extension (remaining in the I or D state). Hence, in this instance, the pair HMM can be seen as alignment, and a CIGAR string may be produced to encode a sequence of the various state transitions.

[0185] In various instances, the3-statebasemodelmaybecomplicatedbyallowing the transitionprobabilities to varyby position. For instance, the probabilities of all M transitionsmay bemultiplied by the prior probabilities of observing the next read base given its base quality score, and the corresponding next haplotype base. In such an instance, the base quality scoresmay translate toaprobability of a sequencingSNPerror.When the twobasesmatch, theprior probability is takenas oneminus this error probability, and when theymismatch, it is taken as the error probability divided by 3, since there are 3 possible SNP results.

[0186] Theabove discussion is regarding an abstract "Markovish"model. In various instances, themaximum-likelihood transition sequence may also be determined, which is termed herein as an alignment, and may be performed using a Needleman-Wunsch or other dynamic programming algorithm. But, in various instances, in performing a variant calling function, as disclosed herein, the maximum likelihood alignment, or any particular alignment, need not be a primary concern. Rather, the total probability may be computed, for instance, by computing the total probability of observing the 32 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 read given the haplotype, which is the sum of the probabilities of all possible transition paths through the graph, from read position zero at any haplotype position, to the read end position, at any haplotype position, each component path probability being simply the product of the various constituent transition probabilities.

[0187] Finding thesumofpathwayprobabilitiesmayalsobeperformedbyemployingavirtual arrayandusingadynamic programming algorithm, as described herein, such that in each cell of a (0 ... N) x (0 ... M)matrix, there are three probability values calculated, corresponding toM, D, and I transition states. (Or equivalently, there are 3matrices.) The top row (read position zero) of thematrix may be initialized to probability 1.0 in the D states, and 0.0 in the I andM states; and the rest of the left column (haplotype position zero) may be initialized to all zeros. (In software, the initial D probabilities may be set near the double-precision max value, e.g. 2^1020, so as to avoid underflow, but this factor may be normalized out later.)

[0188] This3-to‑1computationdependency restricts theorder that cellsmaybecomputed.Theycanbecomputed left to right in each row, progressing through rows from top to bottom, or top to bottom in each column, progressing rightward. Additionally, they may be computed in anti-diagonal wavefronts, where the next step is to compute all cells (n,m) where n+mequals the incremented stepnumber. Thiswavefront order has theadvantage that all cells in theanti-diagonalmaybe computed independently of each other. The bottom row of the matrix then, at the final read position, may be configured to represent the completed alignments. In such an instance, the Haplotype Caller will work by summing the I and M probabilities of all bottom row cells. In various embodiments, the system may be set up so that no D transitions are permitted within the bottom row, or a D transition probability of 0.0 may be used there, so as to avoid double counting.

[0189] As described herein, in various instances, each HMM evaluation may operate on a sequence pair, such as on a candidate haplotype and a read pair. For instance, within a given active region, each of a set of haplotypes may be HMM- evaluated vs. each of a set of reads. In such an instance, the software and / or hardware input bandwidth may be reduced and / orminimized by transferring the set of reads and the set of haplotypes once, and letting the software and / or hardware generate the NxM pair operations. In certain instances, a Smith-Waterman evaluator may be configured to queue up individualHMMoperations, eachwith itsowncopyof readandhaplotypedata.ASmith-Waterman (SW)alignmentmodule may be configured to run the pair HMMcalculation in linear space ormay operate in log probability space. This is useful to keep precision across the huge range of probability values with fixed-point values. However, in other instances, floating point operations may be used.

[0190] There are three parallel multiplications (e.g., additions in log space), then two serial additions (~5‑6 stage approximation pipelines), then an additional multiplication. In such an instance, the full pipeline may be about L = 12‑16 cycles long. The I&Dcalculationsmaybeabout half the length. Thepipelinemaybe fedamultiplicity of input probabilities, suchas2or 3or 5or 7ormore input probabilities eachcycle, suchas fromoneormorealready computedneighboringcells (Mand / orD from the left,Mand / or I fromabove, and / orMand / or I and / orD fromabove-left). Itmayalso includeoneormore haplotype bases, and / or one or more read bases such as with associated parameters, e.g., preprocessed parameters, each cycle. It outputs the M & I & D result set for one cell each cycle, after fall-through latency.

[0191] As indicated above, in performing a variant call function, as disclosed herein, a De Bruijn Graph may be formulated, andwhenall of the reads inapile upare identical, theDBGwill be linear.However,where therearedifferences, thegraphwill form "bubbles" that are indicativeof regionsof differences resulting inmultiple pathsdiverging frommatching the reference alignment and then later re-joining in matching alignment. From this DBG, various paths may be extracted, which form candidate haplotypes, e.g., hypotheses for what the true DNA sequencemay be on at least one strand, which hypotheses may be tested by performing an HMM, or modified HMM, operation on the data. Further still, a genotyping function may be employed such as where the possible diploid combinations of the candidate haplotypes may be formed, and for each of them, a conditional probability of observing the entire read pileup may be calculated. These results may then be fed into a Bayesian formula module to calculate an absolute probability that each genotype is the truth, given the entire read pileup observed.

[0192] Hence, inaccordancewith thedevices, systems,andmethodsof their usedescribedherein, in various instances, a genotyping operationmay be performed, which genotyping operationmay be configured so as to be implemented in an optimized manner in software and / or in hardware and / or by a quantum processing unit. For instance, the possible diploid combinations of the candidate haplotypesmaybe formed, and for each combination, a conditional probability of observing the entire read pileupmay be calculated, such as by using the constituent probabilities of observing each read given each haplotype from the pair HMMevaluation. The results of these calculations feed into a Bayesian formula so as to calculate an absolute probability that each genotype is the truth, given the entire read pileup observed.

[0193] Accordingly, in various aspects, the present disclosure is directed to a system for performing a haplotype or variant call operation on generated and / or supplied data so as to produce a variant call file with respect thereto. Specifically, as described herein above, in particular instances, a variant call file may be a digital or other such file that encodes the difference between one sequence and another, such as the difference between a sample sequence and a reference sequence. Specifically, in various instances, the variant call file may be a text file that sets forth or otherwise details the genetic and / or structural variations in a person’s genetic makeup as compared to one or more reference genomes.

[0194] For instance, a haplotype is a set of genetic, e.g., DNA and / or RNA, variations, such as polymorphisms that 33 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 reside in a person’s chromosomes and as suchmay be passed on to offspring and thereby inherited together. Particularly, a haplotype can refer to a combination of alleles, e.g., one of a plurality of alternative forms of a gene such asmay arise by mutation, which allelic variations are typically found at the same place on a chromosome. Hence, in determining the identity of a person’s genome it is important to know which form of various different possible alleles a specific person’s genetic sequence codes for. In particular instances, a haplotype may refer to one or more, e.g., a set, of nucleotide polymorphisms (e.g., SNPs) that may be found at the same position on the same chromosome.

[0195] Typically, in various embodiments, in order to determine the genotype, e.g., allelic haplotypes, for a subject, as described herein and above, a software based algorithm may be engaged, such as an algorithm employing a haplotype call program, e.g., GATK, for simultaneously determining SNPs and / or insertions and / or deletions, e.g., indels, in an individual’s genetic sequence. In particular, the algorithmmay involve one ormore haplotype assembly protocols such as for local de-novo assembly of a haplotype in one or more active regions of the genetic sequence being processed. Such processing typically involves the deployment of a processing function called a Hidden Markov Model (HMM) that is a stochastic and / or statistical model used to exemplify randomly changing systems such as where it is assumed that future states within the system depend only on the present state and not on the sequence of events that precedes it.

[0196] In such instances, the system being modeled bears the characteristics or is otherwise assumed to be a Markov process with unobserved (hidden) states. In particular instances, the model may involve a simple dynamic Bayesian network. Particularly, with respect to determining genetic variation, in its simplest form, there is one of four possibilities for the identity of any given base in a sequence being processed, such as when comparing a segment of a reference sequence, e.g., a hypothetical haplotype, and that of a subject’s DNA or RNA, e.g., a read derived from a sequencer. However, in order to determinesuch variation, in a first instance, a subject’sDNA / RNAmust be sequenced, e.g., via aNext Gen Sequencer ("NGS"), to produce a readout or "reads" that identify the subject’s genetic code.

[0197] Next, once the subject’s genome has been sequenced to produce one or more reads, the various reads, representative of the subject’s DNA and / or RNA need to be mapped and / or aligned, as herein described above in great detail. The next step in the process then is to determine how the genes of the subject that have just been determined, e.g., having been mapped and / or aligned, vary from that of a prototypical reference sequence. In performing such analysis, therefore, it is assumed that the read potentially representing a given gene of a subject is a representation of the prototypical haplotype albeit with various SNPs and / or indels that are to presently be determined.

[0198] Specifically, in particular aspects, devices, systems, and / or methods for practicing the same, such as for performing a haplotype and / or variant call function, such as deploying an HMM function, for instance, in an accelerated haplotype caller is provided. In various instances, in order to overcome these and other such various problems known in the art, the HMM accelerator herein presented may be configured to be operated in a manner so as to be implemented in software, implemented inhardware, or a combinationofbeing implementedand / orotherwisecontrolled inpart by software and / or in part by hardware and / or may include quantum computing implementations. For instance, in a particular aspect, the disclosure is directed to a method by which data pertaining to the DNA and / or RNA sequence identity of a subject and / or how the subject’s genetic information may differ from that of a reference genome may be determined.

[0199] In such an instance, themethodmay be performed by the implementation of a haplotype or variant call function, suchasemployinganHMMprotocol. Particularly, theHMMfunctionmaybeperformed inhardware, software, or via oneor more quantum circuits, such as on an accelerated device, in accordance with a method described herein. In such an instance, theHMMacceleratormaybe configured to receive andprocess the sequenced,mapped, and / or aligned data, to process the same, e.g., to produce a variant call file, aswell as to transmit the processeddata back throughout the system. Accordingly, the method may include deploying a system where data may be sent from a processor, such as a software- controlled CPU or GPU or even a QPU, to a haplotype caller implementing an accelerated HMM, which haplotype caller may be deployed on amicroprocessor chip, such as an FPGA, ASIC, or structured ASIC or implemented by one or more quantum circuits. The method may further include the steps for processing the data to produce HMM result data, which results may then be fed back to the CPU and / or GPU and / or QPU.

[0200] Particularly, in one embodiment, as can be seen with respect to FIG. 3A, a bioinformatics pipeline system including an HMM accelerator is provided. For instance, in one instance, the bioinformatics pipeline system may be configured as a variant call system 1. The system is illustrated as being implemented in hardware, but may also be implemented via one ormore quantum circuits, such as of a quantum computing platform. Specifically, FIG. 3A provides a high-level view of an HMM interface structure. In particular embodiments, the variant call system 1 is configured to accelerate at least a portion of a variant call operation, such as an HMMoperation. Hence, in various instances, the HMM systemmay be referenced herein as a part of the VC system1. The system1 includes a server having one ormore central processing units (CPU / GPU / QPU) 1000 configured for performing one or more routines related to the sequencing and / or processing of genetic information, such as for comparing a sequenced genetic sequence to one or more reference sequences.

[0201] Additionally, the system1 includes a peripheral device 2, such as an expansion card, that includes amicrochip 7, such as an FPGA, ASIC, or sASIC. In some instances, one or more quantum circuits may be provided and configured for performing the various operations set forth herein. It is also to be noted that the termASICmay refer equally to a structured 34 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 ASIC (sASIC), where appropriate. The peripheral device 2 includes an interconnect 3 and a bus interface 4, such as a parallel or serial bus, which connects the CPU / GPU / QPU 1000with the chip 7. For instance, the device 2may comprise a peripheral component interconnect, such as a PCI, PCI-X, PCIe, or QPI (quick path interconnect), andmay include a bus interface 4, that is adapted to operably and / or communicably connect theCPU / GPU / QPU1000 to the peripheral device 2, such as for low latency, high data transfer rates. Accordingly, in particular instances, the interface may be a peripheral component interconnect express (PCIe) 4 that is associated with the microchip 7, which microchip includes an HMM accelerator 8. For example, in particular instances, the HMM accelerator 8 is configured for performing an accelerated HMM function, such as where the HMM function, in certain embodiments, may at least partially be implemented in the hardware of the FPGA, AISC, or sASIC or via one or more suitably configured quantum circuits.

[0202] Specifically, FIG. 3A presents a high-level figure of an HMM accelerator 8 having an exemplary organization of one or more engines 13, such as a plurality of processing engines 13a - 13m+1, for performing one ormore processes of a variant call function, such as including an HMM task. Accordingly, the HMM accelerator 8 may be composed of a data distributor 9, e.g., CentCom, and one or a multiplicity of processing clusters 11 - 11n+1 that may be organized as or otherwise includeoneormore instances13, suchaswhereeach instancemaybeconfiguredasaprocessingengine, such as a small engine 13a - 13m+1. For instance, the distributor 9 may be configured for receiving data, such as from the CPU / GPU / QPU 1000, and distributing or otherwise transferring that data to one or more of the multiplicity of HMM processing clusters 11.

[0203] Particularly, in certain embodiments, the distributor 9 may be positioned logically between the on-board PCIe interface 4 and theHMMaccelerator module 8, such aswhere the interface 4 communicateswith the distributor 9 such as over an interconnect or other suitably configured bus 5, e.g., PCIe bus. The distributor module 9 may be adapted for communicatingwithoneormoreHMMaccelerator clusters11suchasoveroneormorecluster buses10.For instance, the HMMacceleratormodule 8may be configured as or otherwise include an array of clusters 11a‑11n+1, such aswhere each HMMcluster 11maybe configured as or otherwise includes a cluster hub 11 and / ormay include oneormore instances 13, which instancemaybeconfiguredasaprocessingengine13 that is adapted for performingoneormoreoperationsondata received thereby. Accordingly, in various embodiments, each cluster 11 may be formed as or otherwise include a cluster hub 11a‑11n+1, where each of the hubs may be operably associated with multiple HMM accelerator engine instances 13a‑13m+1, such aswhere each cluster hub 11maybe configured for directing data to a plurality of the processing engines 13a - 13m+1 within the cluster 11.

[0204] In various instances, the HMM accelerator 8 is configured for comparing each base of a subject’s sequenced genetic code, such as in read format, with the various known or generated candidate haplotypes of a reference sequence and determining the probability that any given base at a position being considered either matches or doesn’t match the relevant haplotype, e.g., the read includes anSNP, an insertion, or a deletion, thereby resulting in a variation of the base at the position being considered. Particularly, in various embodiments, the HMM accelerator 8 is configured to assign transition probabilities for the sequence of the bases of the read going between each of these states, Match ("M"), Insert ("I"), or Delete ("D") as set forth in FIG. 2 and as described in greater detail herein below.

[0205] More particularly, dependent on the configuration, the HMMacceleration function may be implemented in either software, such as by the CPU / GPU / QPU 1000 and / or microchip 7, and / or may be implemented in hardware and may be present within the microchip 7, such as positioned on the peripheral expansion card or board 2. In various embodiments, this functionality may be implemented partially as software, e.g., run by the CPU / GPU / QPU 1000, and partially as hardware, implemented on the chip 7 or via one or more quantum processing circuits. Accordingly, in various embodi- ments, thechip7maybepresent on themotherboardof theCPU / GPU / QPU1000,or itmaybepart of theperipheral device 2, or both. Consequently, the HMM accelerator module 8may include or otherwise be associated with various interfaces, e.g., 3, 5, 10, and / or 12 so as to allow the efficient transfer of data to and from the processing engines 13.

[0206] Accordingly, as can be seenwith respect to FIGS. 2 and 3, in various embodiments, amicrochip 7 configured for performing a variant, e.g., haplotype, call function is provided. Themicrochip 7may be associated with a CPU / GPU / QPU 1000 such as directly coupled therewith, e.g., included on the motherboard of a computer, or indirectly coupled thereto, such as being included as part of a peripheral device 2 that is operably coupled to the CPU / GPU / QPU 1000, such as via oneormore interconnects, e.g., 3, 4, 5, 10, and / or 12. In this instance, themicrochip7 ispresent on theperipheral device2. It is to be understood that although configured as a microchip, the accelerator could also be configured as one or more quantum circuits of a quantum processing unit, wherein the quantum circuits are configured as one or more processing engines for performing one or more of the functions disclosed herein.

[0207] Hence, the peripheral device 2 may include a parallel or serial expansion bus 4 such as for connecting the peripheral device 2 to the central processing unit (CPU / GPU / QPU) 1000 of a computer and / or server, such as via an interface 3, e.g., DMA. In particular instances, the peripheral device 2 and / or serial expansion bus 4 may be a Peripheral Component Interconnect express (PCIe) that is configured tocommunicatewithorotherwise include themicrochip7, such as via connection 5. As described herein, themicrochip 7may at least partially be configured as or may otherwise include an HMMaccelerator 8. The HMMaccelerator 8may be configured as part of themicrochip 7, e.g., as hardwired and / or as code to be run in association therewith, and is configured for performing a variant call function, such as for performing one 35 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 or more operations of a Hidden Markov Model, on data supplied to the microchip 7 by the CPU / GPU / QPU 1000, such as over the PCIe interface 4. Likewise, once one ormore variant call functions have been performed, e.g., one ormore HMM operations run, the results thereof may be transferred from the HMM accelerator 8 of the chip 7 over the bus 4 to the CPU / GPU / QPU 1000, such as via connection 3.

[0208] For instance, in particular instances, a CPU / GPU / QPU 1000 for processing and / or transferring information and / orexecuting instructions isprovidedalongwithamicrochip7 that isat least partially configuredasanHMMaccelerator 8. The CPU / GPU / QPU 1000 communicates with the microchip 7 over an interface 5 that is adapted to facilitate the communication between the CPU / GPU / QPU 1000 and the HMM accelerator 8 of the microchip 7 and therefore may communicably connect theCPU / GPU / QPU1000 to theHMMaccelerator 8 that ispart of themicrochip7.To facilitate these functions, themicrochip 7 includesadistributormodule9,whichmaybeaCentCom, that is configured for transferringdata to amultiplicity of HMMengines 13, e.g., via one ormore clusters 11, where each engine 13 is configured for receiving and processing the data, such as by running an HMM protocol thereon, computing final values, outputting the results thereof, and repeating the same. In various instances, the performance of an HMMprotocol may include determining one ormore transition probabilities, as described herein below. Particularly, each HMM engine 13 may be configured for performing a job suchas includingoneormoreof thegeneratingand / or evaluatingof anHMMvirtualmatrix to produceandoutput a final sum value with respect thereto, which final sum expresses the probable likelihood that the called base matches or is different from a corresponding base in a hypothetical haplotype sequence, as described herein below.

[0209] FIG. 3B presents a detailed depiction of the HMM cluster 11 of FIG. 3A. In various embodiments, each HMM cluster 11 includes one or more HMM instances 13. One or a number of clusters may be provided, such as desired in accordance with the amount of resources provided, such as on the chip or quantum computing processor. Particularly, a HMM cluster may be provided, where the cluster is configured as a cluster hub 11. The cluster hub 11 takes the data pertaining to one or more jobs 20 from the distributor 9, and is further communicably connected to one or more, e.g., a plurality of, HMM instances13, suchas via oneormoreHMM instancebusses 12, towhich the cluster hub 11 transmits the job data 20.

[0210] The bandwidth for the transfer of data throughout the systemmay be relatively low bandwidth process, and once a job 20 is received, the system 1 may be configured for completing the job, such as without having to go off chip 7 for memory. In variousembodiments, one job20a is sent to oneprocessingengine13aat anygiven time, but several jobs20a- n may be distributed by the cluster hub 11 to several different processing engines 13a‑13m+1, such as where each of the processing engines 13will beworking on a single job 20, e.g., a single comparison between one ormore reads and one or more haplotype sequences, in parallel and at high speeds.

[0211] As described below, the performance of such a job 20 may typically involve the generation of a virtual matrix whereby the subject’s "read" sequences may be compared to one or more, e.g., two, hypothetical haplotype sequences, so as to determine the differences there between. In such instances, a single job 20may involve the processing of one or morematrices having amultiplicity of cells therein that need to be processed for each comparison beingmade, such as on abasebybasebasis. As thehumangenome is about 3 billion basepairs, theremaybeon theorder of 1 to 2billion different jobs tobeperformedwhenanalyzinga30Xoversamplingof ahumangenome (which isequitable toabout 20 trillion cells in the matrices of all associated HMM jobs).

[0212] Accordingly, as described herein, each HMM instance 13 may be adapted so as to perform an HMM protocol, e.g., the generating and processing of an HMM matrix, on sequence data, such as data received thereby from the CPU / GPU / QPU1000. For example, asexplainedabove, in sequencinga subject’s geneticmaterial, suchasDNAorRNA, the DNA / RNA is broken down into segments, such as up to about 100 bases in length. The identity of these 100 base segments are then determined, such as by an automated sequencer, and "read" into a FASTQ text based file or other format that stores both each base identity of the read along with a Phred quality score (e.g., typically a number between 0 and 63 in log scale, where a score of 0 indicates the least amount of confidence that the called base is correct, with scores between 20 to 45 generally being acceptable as relatively accurate).

[0213] Particularly, as indicated above, a Phred quality score is a quality indicator that measures the quality of the identification of the nucleobase identities generated by the sequencing processor, e.g., by the automated DNA / RNA sequencer. Hence, each readbase includes its ownquality, e.g., Phred, scorebasedonwhat the sequencer evaluated the quality of that specific identification to be. ThePhred represents the confidencewith which the sequencer estimates that it got the called base identity correct. This Phred score is then used by the implemented HMM module 8, as described in detail below, to further determine the accuracy of each called base in the read as compared to the haplotype towhich it has beenmapped and / or aligned, such as by determining its Match, Insertion, and / or Deletion transition probabilities, e.g., in and out of theMatch state. It is to be noted that in various embodiments, the system 1maymodify or otherwise adjust the initial Phred score prior to the performance of an HMM protocol thereon, such as by taking into account neighboring bases / scores and / or fragments of neighboring DNA and allowing such factors to influence the Phred score of the base, e.g., cell, under examination.

[0214] In such instances, as can be seen with respect to FIG. 3A and 3B, the system 1, e.g., computer / quantum software, may determine and identify various active regions 500n within the sequenced genome that may be explored 36 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 and / or otherwise subjected to further processingasherein described,whichmaybebrokendown into jobs20n thatmaybe parallelized amongst the various cores and available threads 1007 throughout the system 1. For instance, such active regions 500may be identified as being sources of variation between the sequenced and reference genomes. Particularly, the CPU / GPU / QPU 1000 may have multiple threads 1007 running, identifying active regions 500a, 500b, and 500c, compiling and aggregating various different jobs 20 n to be worked on, e.g., via a suitably configured aggregator 1008, basedon theactive region(s) 500a-c currently being examined. Any suitable number of threads 1007maybeemployed so as to allow the system 1 to run at maximum efficiency, e.g., the more threads present the less active time spent waiting.

[0215] Once identified, compiled, and / or aggregated, the threads 1007 / 1008 will then transfer the active jobs 20 to the data distributor 9, e.g., CentCom, of theHMMmodule 8, such as via PCIe interface 4, e.g., in a fire and forgetmanner, and will thenmove on to a different processwhile waiting for theHMM8 to send the output data back so as to bematched back up to the corresponding active region 500 to which it maps and / or aligns. The data distributor 9 will then distribute the jobs 20 to the various different HMMclusters 11, such as on a job-by-jobmanner. If everything is running efficiently, thismay be on a first in first out format, but suchdoes not need to be the case. For instance, in various embodiments, raw jobs data and processed job results data may be sent through and across the system as they become available.

[0216] Particularly, as can be seen with respect to FIGS. 2, 3, and 4, the various job data 20may be aggregated into 4K byte pages of data, whichmay be sent via the PCIe 4 to and through the CentCom 9 and on to the processing engines 13, e.g., via the clusters 11. The amount of data being sent may bemore or less than 4K bytes, but will typically include about 100 HMM jobs per 4K (e.g., 1024) page of data. Particularly, these data then get digested by the data distributor 9 and are fed to each cluster 11, such as where one 4K page is sent to one cluster 11. However, such need not be the case as any given job 20 may be sent to any given cluster 11, based on the clusters that become available and when.

[0217] Accordingly, the cluster 11 approach as presented here efficiently distributes incoming data to the processing engines 13 at high-speed. Specifically, as data arrives at the PCIe interface 4 from the CPU / GPU / QPU 1000, e.g., over DMAconnection3, the receiveddatamay thenbesent over thePCIebus5 to theCentComdistributor 9of thevariant caller microchip 7. The distributor 9 then sends the data to one or more HMM processing clusters 11, such as over one or more cluster dedicated buses 10, which cluster 11may then transmit the data to one or more processing instances 13, e.g., via one or more instance buses 12, such as for processing. In this instance, the PCIe interface 4 is adapted to provide data through theperipheral expansionbus5, distributor 9, and / or cluster 10and / or instance12bussesat a rapid rate, suchasat a rate that can keep one or more, e.g., all, of the HMM accelerator instances 13a‑(m+1) within one or more, e.g., all, of the HMMclusters11a‑(n+1) busy, suchasoveraprolongedperiodof time,e.g., full time,during theperiodoverwhich thesystem 1 is being run, the jobs 20 are being processed, andwhilst also keeping upwith the output of the processedHMMdata that is to be sent back to one or more CPUs 1000, over the PCIe interface 4.

[0218] For instance, any inefficiency in the interfaces3, 5, 10, and / or 12 that leads to idle time for oneormoreof theHMM accelerator instances 13 may directly add to the overall processing time of the system 1. Particularly, when analyzing a human genome, theremay be on the order of two ormore billion different jobs 20 that need to be distributed to the various HMM clusters 11 and processed over the course of a time period, such as under 1 hour, under 45 minutes, under 30 minutes, under 20 minutes including 15 minutes, 10 minutes, 5 minutes, or less.

[0219] Accordingly, FIG. 4 sets forth an overview of an exemplary data flow throughout the software and / or hardware of the system1, as described generally above. As can be seenwith respect to FIG. 4, the system1may be configured in part to transfer data, such as between the PCIe interface 4 and the distributor 9, e.g., CentCom, such as over the PCIe bus 5. Additionally, the system1may further be configured in part to transfer the received data, such as between the distributor 9 and the oneormoreHMMclusters 11, such as over the oneormore cluster buses 10. Hence, in various embodiments, the HMMaccelerator 8may include one ormore clusters 11, such as one ormore clusters 11 configured for performing one or more processes of anHMM function. In such an instance, there is an interface, such as a cluster bus 10, that connects the CentCom 9 to the HMM cluster 11.

[0220] For instance,FIG. 5 is ahigh-level diagramdepicting the interface in to andout of theHMMmodule8, suchas into andoutof aclustermodule.Ascanbeseenwith respect toFIG.6, eachHMMcluster11maybeconfigured tocommunicate with, e.g., receive data from and / or send final result data, e.g., sum data, to the CentCom data distributor 9 through a dedicatedcluster bus10.Particularly, anysuitable interfaceorbus5maybeprovidedso longas it allows thePCIe interface 4 to communicate with the data distributor 9. More particularly, the bus 5 may be an interconnect that includes the interpretation logic useful in talking to the data distributor 9, which interpretation logicmay be configured to accommodate any protocol employed to provide this functionality. Specifically, in various instances, the interconnect may be configured as a PCIe bus 5.

[0221] Additionally, the cluster 11 may be configured such that single or multiple clock domains may be employed therein, andhence, oneormore clocksmaybepresentwithin the cluster 11. Inparticular instances,multiple clockdomains may be provided. For example, a slower clockmay be provided, such as for communications, e.g., to and from the cluster 11. Additionally, a faster, e.g., a high speed, clock may be provided which may be employed by the HMM instances 13 for use in performing the various state calculations described herein.

[0222] Particularly, in various embodiments, as canbe seenwith respect toFIG. 6, the system1maybe set up such that, 37 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 in a first instance, as the data distributor 9 leverages the existing CentCom IP, a collar, such as a gasket, may be provided, where the gasket is configured for translating signals to and from the CentCom interface 5 from and to the HMM cluster interface or bus 10. For instance, anHMMcluster bus 10maycommunicably and / or operably connect theCPU / GPU1000 to the various clusters 11 of the HMMacceleratormodule 8. Hence, as can be seenwith respect to FIG. 6, structuredwrite and / or read data for each haplotype and / or for each read may be sent throughout the system 1.

[0223] Following a job 20 being input into theHMMengine, anHMMengine 13may typically start either: a) immediately, if it is IDLE, orb)after it hascompleted its currently assigned task. It is tobenoted that eachHMMaccelerator engine13can handle ping and pong inputs (e.g., can be working on one data set while the other is being loaded), thus minimizing downtime between jobs. Additionally, the HMM cluster collar 11 may be configured to automatically take the input job 20 sent by the data distributor 9 and assign it to one of the HMM engine instances 13 in the cluster 11 that can receive a new job.Thereneednot beacontrol on thesoftware side that canselect a specificHMMengine instance13 fora specific job20. However, in various instances, the software can be configured to control such instances.

[0224] Accordingly, in view of the above, the system1may be streamlinedwhen transferring the results data back to the CPU / GPU / QPU, and because of this efficiency there is not much data that needs to go back to the CPU / GPU / QPU to achieve the usefulness of the results. This allows the system to achieve about a 30 minute or less, such as about a 25 or about a 20minute or less, for instance, about a 18 or about a 15minute or less, including about a 10 or about a 7minute or less, even about a 5 or about a 3 minute or less variant call operation, dependent on the system configuration.

[0225] FIG. 6 presents a high-level view of various functional blocks within an exemplary HMM engine 13 within a hardware accelerator 8, on the FPGA or ASIC 7. Specifically, within the hardware HMM accelerator 8 there are multiple clusters 11, and within each cluster 11 there are multiple engines 13. FIG. 6 presents a single instance of an HMM engine 13. As can be seenwith respect to FIG. 6, the engine 13may include an instance bus interface 12, a plurality ofmemories, e.g., an HMEM 16 and an RMEM 18, various other components 17, HMM control logic 15, as well as a result output interface 19. Particularly, on the engine side, the HMM instance bus 12 is operably connected to thememories, HMEM16 andRMEM18, andmay include interface logic that communicateswith the cluster hub 11,which hub is in communications with the distributor 9, which in turn is communicating with the PCIe interface 4 that communicates with the variant call software being run by the CPU / GPU and / or server 1000. The HMM instance bus 12, therefore, receives the data from the CPU 1000 and loads it into one or more of the memories, e.g., the HMEM and RMEM. This configuration may also be implemented in one or more quantum circuits and adapted accordingly.

[0226] In these instances, enoughmemory space should be allocated such that at least one or two ormore haplotypes, e.g., twohaplotypes,maybe loaded, e.g., in theHMEM16, per given readsequence that is loaded, e.g., into theRMEM18, which when multiple haplotypes are loaded results in an easing of the burden on the PCIe bus 5 bandwidth. In particular instances, two haplotypes and two read sequencesmay be loaded into their respective memories, which would allow the four sequences to be processed together in all relevant combinations. In other instances four, or eight, or sixteen sequences, e.g., pairs of sequences, may be loaded, and in like manner be processed in combination, such as to further ease the bandwidth when desired.

[0227] Additionally, enoughmemorymaybe reservedsuch that aping-pongstructuremaybe implemented therein such that once thememoriesare loadedwithanew job20a, suchason thepingsideof thememory, anew jobsignal is indicated, and the control logic 15 may begin processing the new job 20a, such as by generating the matrix and performing the requisite calculations, as describedherein andbelow.Accordingly, this leaves thepongsideof thememoryavailable soas to be loaded up with another job 20b, which may be loaded therein while the first job 20a is being processed, such that as the first job 20a is finished, the second job 20b may immediately begin to be processed by the control logic 15.

[0228] In suchan instance, thematrix for job20bmaybepreprocessedso that there is virtually nodown time, e.g., oneor two clock cycles, from the ending of processing of the first job 20a, and the beginning of processing of the second job 20b. Hence, when utilizing both the ping and pong side of thememory structures, theHMEM16may typically store 4 haplotype sequences, e.g., two a piece, and the RMEM 18 may typically store 2 read sequences. This ping-pong configuration is useful because it simply requires a little extra memory space, but allows for a doubling of the throughput of the engine 13.

[0229] During and / or after processing the memories 16, 18 feed into the transition probabilities calculator and lookup table (LUT) block 17a, which is configured for calculating various information related to "Priors" data, as explained below, which in turn feeds the Prior results data into theM, I, and D state calculator block 17b, for use when calculating transition probabilities. One or more scratch RAMs 17c may also be included, such as for holding the M, I, and D states at the boundary of the swath, e.g., the valuesof the bottom rowof the processing swath,which as indicated, in various instances, may be any suitable amount of cells, e.g., about 10 cells, in length so as to be commensurate with the length of the swath 35.

[0230] Additionally, a separate resultsoutput interfaceblock19maybe includedso thatwhen thesumsarefinished they, e.g., a 432-bitword, can immediately be transmittedback to the variant call softwareof theCPU / GPU / QPU1000. It is tobe noted that this configurationmay be adapted so that the system 1, specifically theM, I, and D calculator 17b is not held up waiting for the output interface 19 to clear, e.g., so long as it does not take as long to clear the results as it does to perform the job20.Hence, in this configuration, theremaybe threepipeline steps functioning in concert tomakeanoverall systems 38 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 pipeline, such as loading the memory, performing the MID calculations, and outputting the results. Further, it is noted that any givenHMMengine 13 is one ofmanywith their own output interface 19, however theymay share a common interface 10 back to the data distributor 9. Hence, the cluster hub 11 will include management capabilities to manage the transfer ("xfer") of information through the HMM accelerator 8 so as to avoid collisions.

[0231] Accordingly, the followingdetails the processesbeing performedwithin eachmoduleof theHMMengines13as it receives thehaplotypeand readsequencedata, processes it, andoutputs results datapertaining to thesame,asgenerally outlined above. Specifically, the high-bandwidth computations in the HMM engine 13, within the HMM cluster 11, are directed to computing and / or updating the match (M), insert (I), and delete (D) state values, which are employed in determining whether the particular read being examined matches the haplotype reference as well as the extent of the same, as described above.

[0232] Particularly, the read along with the Phred score and GOP value for each base in the read is transmitted to the cluster 11 from thedistributor 9 and is therebyassigned toaparticular processingengine13 for processing. Thesedataare then used by the M, I, and D calculator 17 of the processing engine 13 to determine whether the called base in the read is more or less likely to be correct and / or to be a match to its respective base in the haplotype, or to be the product of a variation, e.g., an insert or deletion; and / or if there is avariation,whether suchvariation is the likely result of a truevariability in the haplotype or rather an artifact of an error in the sequence generating and / or mapping and / or aligning systems.

[0233] As indicated above, a part of suchanalysis includes theMIDcalculator 17 determining the transition probabilities fromonebase to another in the read going fromoneM, I, orD state to another in comparison to the reference, such as from amatching state to anothermatching state, or amatching state to either an insertion state or to a deletion state. Inmaking such determinations each of the associated transition probabilities is determined and considered when evaluating whether any observed variation between the read and the reference is a true variation and not just some machine or processing error. For these purposes, the Phred score for each base being considered is useful in determining the transition probabilities in and out of the match state, such as going from a match state to an insert or deletion, e.g., a gapped, state in the comparison. Likewise, the transition probabilities of continuing a gapped state or going fromagapped state, e.g., an insert or deletion state, back to amatch state are also determined. In particular instances, the probabilities in or out of thedelete or insert state, e.g., exitingagapcontinuation state,maybeafixedvalue, andmaybe referencedherein as the gap continuation probability or penalty. Nevertheless, in various instances, such gap continuation penaltiesmay be floating and therefore subject to change dependent on the accuracy demands of the system configuration.

[0234] Accordingly, as depictedwith respect to FIGS. 7 and 8 each of theM, I, andD state values are computed for each possible read and haplotype base pairing. In such an instance, a virtual matrix 30 of cells containing the read sequence beingevaluatedononeaxis of thematrix and theassociatedhaplotypesequenceon theotheraxismaybe formed, suchas where each cell in the matrix represents a base position in the read and haplotype reference. Hence, if the read and haplotypesequencesareeach100bases in length, thematrix 30will include100by100cells, agivenportionofwhichmay need to be processed in order to determine the likelihood and / or extent to which this particular read matches up with this particular reference. Hence, once virtually formed, the matrix 30 may then be used to determine the various state transitions that take placewhenmoving fromone base in the read sequence to another and comparing the same to that of thehaplotypesequence, suchasdepicted inFIGS.7and8.Specifically, theprocessingengine13 is configuredsuch that a multiplicity of cellsmay be processed in parallel and / or sequential fashionwhen traversing thematrix with the control logic 15. For instance, as depicted inFIG. 7, a virtual processing swath 35 is propagated andmovesacross anddown thematrix 30, such as from left to right, processing the individual cells of the matrix 30 down the right to left diagonal.

[0235] Morespecifically, as canbeseenwith respect toFIG.7, each individual virtual cellwithin thematrix 30 includesan M, I, and D state value that needs to be calculated so as to assess the nature of the identity of the called base, and as depicted in FIG. 7 the data dependencies for each cell in this processmay clearly be seen. Hence, for determining a given M state of a present cell being processed, theMatch, Insert, andDelete states of the cell diagonally above the present cell need to bepushed into thepresent cell and used in the calculation of theMstate of the cell presently being calculated (e.g., thus, the diagonal downwards, forwards progression through the matrix is indicative of matching).

[0236] However, for determining the I state, only the Match and Insert states for the cell directly above the present cell need be pushed into the present cell being processed (thus, the vertical downwards "gapped" progression when continuing in an insertion state). Likewise, for determining the D state, only the Match and Delete states for the cell directly left of the present cell need be pushed into the present cell (thus, the horizontal cross-wards "gapped" progression whencontinuing inadeletion state). Ascanbeseenwith respect toFIG. 7, after computationof cell 1 (the shadedcell in the top most row) begins, the processing of cell 2 (the shaded cell in the second row) can also begin, without waiting for any results fromcell 1, because there is nodata dependencies between this cell in row2and the cell of row1where processing begins. This forms a reverse diagonal 35 where processing proceeds downwards and to the left, as shown by the arrow. This reverse diagonal 35 processing approach increases the processing efficiency and throughput of the overall system. Likewise, thedatagenerated in cell 1, can immediatelybepushed forward to thecell downand forward to the rightof the top most cell 1, thereby advancing the swath 35 forward.

[0237] For instance, FIG. 7 depicts an exemplary HMMmatrix structure 35 showing the hardware processing flow. The 39 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 matrix 35 includes the haplotype base index, e.g., containing 36 bases, positioned to run along the top edge of the horizontal axis, and further includes thebase read index, e.g., 10bases, positioned to fall along thesideedgeof thevertical axis in such a manner to from a structure of cells where a selection of the cells may be populated with an M, I, and D probability state, and the transition probabilities of transitioning from the present state to a neighboring state. In such an instance, as described in greater detail above, a move from a match state to a match state results in a forwards diagonal progression through the matrix 30, while moving from a match state to an insertion state results in a vertical downwards progressing gap, and a move from a match state to a deletion state results in a horizontal progressing gap. Hence, as depicted in FIG. 8, for a given cell, when determining the match, insert, and delete states for each cell, the match, insert, and delete probabilities of its three adjoining cells are employed.

[0238] Thedownwardsarrow inFIG. 7 represents theparallel and sequential nature of theprocessing engine(s) that are configuredsoas toproduceaprocessingswathorwave35 thatmovesprogressivelyalong thevirtualmatrix inaccordance with thedata dependencies, seeFIGS. 7and8, for determining theM, I, andDstates for eachparticular cell in the structure 30.Accordingly, in certain instances, itmaybedesirable tocalculate the identitiesof eachcell in adownwardsanddiagonal manner, as explained above, rather than simply calculating each cell along a vertical or horizontal axis exclusively, although this can be done if desired. This is due to the increased wait time, e.g., latency, that would be required when processing the virtual cells of thematrix 35 individually and sequentially along the vertical or horizontal axis alone, such as via the hardware configuration.

[0239] For instance, in suchan instance,whenmoving linearlyandsequentially through thevirtualmatrix30, suchas ina row by row or column by columnmanner, in order to process each new cell the state computations of each preceding cell would have to be completed, thereby increasing latency time overall. However, when propagating theM, I, D probabilities of each new cell in a downwards and diagonal fashion, the system 1 does not have to wait for the processing of its preceding cell, e.g., of row one, to complete before beginning the processing of an adjoining cell in row two of the matrix. Thisallows for parallel andsequential processingof cells inadiagonal arrangement tooccur, and further allows thevarious computational delays of the pipeline associated with the M, I, and D state calculations to be hidden. Accordingly, as the swath 35 moves across the matrix 30 from left to right, the computational processing moves diagonally downwards, e.g., towards the left (as shown by the arrow in FIG. 7). This configuration may be particularly useful for hardware and / or quantum circuit implementations, such as where the memory and / or clock-by-clock latency are a primary concern.

[0240] In these configurations, the actual value output fromeach cell of anHMMengine 13, e.g., after having calculated the entire matrix 30, may be a bottom row (e.g., Row 35 of FIG. 16) containingM, I, and D states, where theM and I states may be summed (the D states may be ignored at this point having already fulfilled their function in processing the calculations above), so as to produce a final sum value that may be a single probability that estimates, for each read and haplotype index, the probability of observing the read, e.g., assuming the haplotype was the true original DNA sampled.

[0241] Particularly, the outcomeof the processingof thematrix 30, e.g., of FIG. 7,maybeasingle value representing the probability that the read is an actual representation of that haplotype. This probability is a value between 0 and 1 and is formed by summing all of the M and I states from the bottom row of cells in the HMMmatrix 30. Essentially, what is being assessed is the possibility that something could have gone wrong in the sequencer, or associated DNA preparation methodsprior to sequencing, soas to incorrectly produceamismatch, insertion, or deletion into the read that is not actually present within the subject’s genetic sequence. In such an instance, the read is not a true reflection of the subject’s actual DNA.

[0242] Hence, accounting for such production errors, it can be determinedwhat any given read actually represents with respect to the haplotype, and thereby allows the system to better determine how the subject’s genetic sequence, e.g., en masse, may differ from that of a reference sequence. For instance, many haplotypes may be run against many read sequences, generating scores for all of them, and determining based on which matches have the best scores, what the actual genomic sequence identity of the individual is and / or how it truly varies from a reference genome.

[0243] Moreparticularly, FIG. 8depicts anenlargedviewof aportionof theHMMstatematrix 30 fromFIG.7.Asshown in FIG. 8, given the internal composition of each cell in thematrix 30, aswell as the structure of thematrix as awhole, theM, I, and D state probability for any given "new" cell being calculated is dependent on the M, I, and D states of several of its surrounding neighbors that have already been calculated. Particularly, as shown in greater detail with respect to FIGS. 1 and 16, in an exemplary configuration, there may be an approximately a .9998 probability of going from a match state to anothermatch state, and theremay be only a .0001 probability (gap open penalty) of going from amatch state to either an insertionor adeletion, e.g., gapped, state. Further,when ineither agapped insertionor gappeddeletionstate theremaybe only a0.1probability (gapextensionor continuationpenalty) of staying in that gappedstate,while there is a .9 probability of returning to a match state. It is to be noted that according to this model, all of the probabilities in to or out of a given state should sum to one. Particularly, the processing of the matrix 30 revolves around calculating the transition probabilities, accounting for the various gap open or gap continuation penalties and a final sum is calculated.

[0244] Hence, these calculated state transition probabilities are derived mainly from the directly adjoining cells in the matrix 30, such as from the cells that are immediately to the left of, the top of, and diagonally up and left of that given cell presently being calculated, as seen in FIGS. 8 and 16. Additionally, the state transition probabilitiesmay in part be derived 40 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 from the "Phred" quality score that accompanies each read base. These transition probabilities, therefore, are useful in computing theM, I, andDstate values for that particular cell, and likewise for any associatednewcell being calculated. It is to benoted that as describedherein, the gapopenandgap continuationpenaltiesmaybefixedvalues, however, in various instances, the gap open and gap continuation penalties may be variable and therefore programmable within the system, albeit by employing additional hardware resources dedicated to determining such variable transition probability calcula- tions. Such instancesmay be useful where greater accuracy is desired. Nevertheless, when such values are assumed to be constant, smaller resource usage and / or chip sizemay be achieved, leading to greater processing speed, as explained below.

[0245] Accordingly, there is a multiplicity of calculations and / or other mathematical computations, such as multi- plications and / or additions, which are involved in deriving each new M, I, and D state value. In such an instance, such as for calculating maximum throughput, the primitive mathematical computations involved in each M, I, and D transition state calculation may be pipelined. Such pipelining may be configured in a way that the corresponding clock frequencies are high, but where the pipeline depthmay benon-trivial. Further, such a pipelinemay be configured to have a finite depth, and in such instances it may take more than one clock cycle to complete the operations.

[0246] For instance, these computations may be run at high speeds inside the processor 7, such as at about 300MHz. This may be achieved such as by pipelining the FPGA or ASIC heavily with registers so little mathematical computation occurs between each flip-flop. This pipeline structure results in multiple cycles of latency in going from the input of the matchstate to theoutput, butgiven the reversediagonal computingstructure, set forth inFIG.7above, these latenciesmay be hidden over the entire HMM matrix 30, such as where each cell represents one clock cycle.

[0247] Hence, the number ofM, I, andDstate calculationsmaybe limited. In suchan instance, the processing engine13 may be configured in such a manner that a grouping, e.g., swath 35, of cells in a number of rows of the matrix 30 may be processed as a group (such as in a down-and-left-diagonal fashion as illustrated by the arrow in FIG. 7) before proceeding to the processing of a second swath below, e.g., where the second swath contains the same number of cells in rows to be processed as the first. In amanner such as this, a hardware implementation of an accelerator 8, as described herein, may be adapted so as to make the overall system more efficient, as described above.

[0248] Particularly, FIG. 9 sets forth an exemplary computational structure for performing the various state processing calculations herein described.More particularly, FIG. 9 sets forth three dedicated logic blocks 17 of the processing engine 13 for computing the state computations involved in generating each M, I, and D state value for each particular cell, or groupingof cells, beingprocessed in theHMMmatrix 30.These logic blocksmaybe implemented inhardware, but in some instances, may be implemented in software, such as for being performed by one or more quantum circuits.

[0249] As can be seen with respect to FIG. 9, thematch state computation 15a is more involved than either of the insert 15bor deletion 15c computations, this is because in calculating thematch state 15aof the present cell being processed, all of the previousmatch, insert, and delete states of the adjoining cells alongwith various other, e.g., prior, data are included in the present match computation, whereas only the match and either the insert and delete states are included in their respective calculations. Hence, as can be seen with respect to FIG. 9, in calculating amatch state, three statemultipliers, as well as two adders, and a final multiplier, which accounts for the prior, e.g., Phred, data are included. However, for calculating the I orDstate, only twomultipliers andoneadder are included. It is noted that in hardware,multipliers aremore resource intensive than adders.

[0250] Accordingly, to various extents, theM, I, andD state values for processing each new cell in the HMMmatrix uses the knowledge or pre-computation of the following values, such as the "previous"M, I, andD state values from left, above, and / or diagonally left and above of the currently-being-computed cell in the HMM matrix. Additionally, such values representing the prior information, or "priors", may at least in part be based on the "Phred" quality score, and whether the read base and the reference base at a given cell in the matrix 30 match or are different. Such information is particularly useful when determining a match state. Specifically, as can be seen with respect to FIG. 9, in such instances, there are basically seven "transition probabilities" (M-to-M, I-to-M, D-to-M, I-to-I, M-to-I, D-to-D, and M-to-D) that indicate and / or estimate the probability of seeing a gap open, e.g., of seeing a transition from a match state to an insert or delete state; seeing a gap close; e.g., going from an insert or delete state back to amatch state; and seeing the next state continuing in the same state as the previous state, e.g., Match-to-Match, Insert-to-Insert, Delete-to-Delete.

[0251] Thestate values (e.g., in any cell to beprocessed in theHMMmatrix 30), Priors, and transitionprobabilities areall values in the range of [0,1]. Additionally, there are also known starting conditions for cells that are on the left or top edge of the HMMmatrix. As can be seen from the logic 15a of FIG. 9, there are four multiplication and two addition computations that may be employed in the particular M state calculation being determined for any given cell being processed. Likewise, as can be seen from the logic of 15b and 15c there are two multiplications and one addition involved for each I state and eachDstatecalculation, respectively.Collectively, alongwith thepriorsmultiplier this sums toa total ofeightmultiplications and four addition operations for theM, I, andDstate calculations associatedwith each single cell in theHMMmatrix 8 to be processed.

[0252] The final sum output of the computation of the matrix, e.g., for a single job of comparing one read to one or two haplotypes, is the summation of the final M and I states across the entire bottom row of the matrix, which is the final sum 41 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 value that is output from theHMMaccelerator 8 and delivered to theCPU / GPU / QPU. This final summed value represents howwell the readmatches the haplotype(s). The value is a probability, e.g., of less than one, for a single job that may then be compared to the output resulting from another job such as form the same active region 500. It is noted that there are on the order of 20 trillion HMM cells to evaluate in a "typical" human genome at 30X coverage, where these 20 trillion HMM cells are spread across about 1 to 2 billion HMM matrices of all associated HMM jobs.

[0253] The results of such calculations may then be compared one against the other so as to determine, in a more precise manner, how the genetic sequence of a subject differs, e.g., on a base by base comparison, from that of one or more reference genomes. For the final sum calculation, the adders already employed for calculating the M, I, and / or D states of the individual cells may be re-deployed so as to compute the final sum value, such as by including a mux into a selection of the re-deployed adders thereby including one last additional row, e.g., with respect to calculation time, to the matrix soas to calculate this final sum,which if the read length is 100basesamounts to about a1%overhead. In alternative embodiments, dedicatedhardware resources canbeused for performing such calculations. In various instances, the logic for the adders for theMandD state calculationsmaybe deployed for calculating the final sum,whichDstate addermaybe efficiently deployed since it is not otherwise being used in the final processing leading to the summing values.

[0254] In certain instances, these calculations and relevant processes may be configured so as to correspond to the output of agivensequencingplatform, suchas includinganensembleof sequencers,whichasacollectivemaybecapable of outputting (on average) a new human genome at 30x coverage every 28 minutes (though they come out of the sequencer ensemble in groups of about 150 genomes every three days). In such an instance, when the presentmapping, aligning, and variant calling operations are configured to fit within such a sequencing platformof processing technologies, a portion of the 28minutes (e.g., about 10minutes) it takes for the sequencing cluster to sequenceagenome,maybeused by a suitably configuredmapper and / or aligner, as herein described, so as to take the image / BCL / FASTQ file results from the sequencer, such as streaming real-time, e.g., on the fly, and perform the steps ofmapping and / or aligning the genome, e.g., post-sequencer processing.

[0255] This leaves about 18 minutes of the sequencing time period for performing the variant calling step, of which the HMM operation is the main computational component, such as prior to the nucleotide sequencer sequencing the next genome, such as over the next 28minutes,where during the sequencing process, generated datamaybe streamed, such as substantially real-time into the present system, such as via the cloud, for instance, for processing to begin on the fly. Accordingly, in such instances, 18 minutes may be budgeted to computing the 20 trillion HMM cells that need to be processed in accordancewith theprocessing of a genome, suchaswhere eachof theHMMcells to beprocessed includes about twelve mathematical operations (e.g., eight multiplications and / or four addition operations). Such a throughput allows for the following computational dynamics (20 trillion HMM cells) x (12 math ops per cell) / (18 minutes x 60 seconds / minute), which is about 222 billion operations per second of sustained throughput.

[0256] FIG. 10 sets forth the logic blocks 17 of the processing engine of FIG. 9 including exemplary M, I, and D state update circuits that present a simplification of the circuit provided in FIG. 9. The systemmay be configured so as to not be memory-limited, so a single HMM engine instance 13 (e.g., that computes all of the single cells in the HMMmatrix 30 at a rate of onecell per clock cycle, onaverage, plusoverheads)maybe replicatedmultiple times (at least 65~70 times tomake the throughputefficient, asdescribedabove).Nevertheless, tominimize thesizeof thehardware, e.g., thesizeof thechip2 and / or its associated resource usage, and / or in a further effort to include asmany HMMengine instances 13 on the chip 2 asdesirable and / or possible, simplificationsmaybemadewith regard to the logic blocks15a’-c’ of the processing instance 13 for computing one or more of the transition probabilities to be calculated.

[0257] In particular, it may be assumed that the gap open penalty (GOP) and gap continuation penalty (GCP), as describedabove, suchas for inserts anddeletesare the sameandare knownprior to chip configuration. This simplification implies that the I-to-M and D-to-M transition probabilities are identical. In such an instance, one or more of themultipliers, e.g., set forth inFIG. 9,maybeeliminated, suchasbypre-adding I andDstatesbeforemultiplying bya common Indel-to-M transition probability. For instance, in various instances, if the I andD state calculations are assumed to be the same, then thestate calculationsper cell canbesimplifiedaspresented inFIG.10.Particularly, if the I andDstatevaluesare thesame, then the I state and the D state may be added and then that sum may be multiplied by a single value, thereby saving a multiply. This may be done because, as seen with respect to FIG. 10, the gap continuation and / or close penalties for the I and D states are the same. However, as indicated above, the system can be configured to calculate different values for both the I and D transition state probabilities, and in such an instance, this simplification would not be employed.

[0258] Additionally, in a further simplification, rather than dedicate chip or other computing resources configured specifically to perform the final sum operation at the bottom of the HMM matrix, the present HMM accelerator 8 may be configured so as to effectively append one or more additional rows to the HMMmatrix 30, with respect to computational time, e.g., overhead, it takes to perform the calculation, and may also be configured to "borrow" one or more adders from the M-state 15a and D-state 15c computation logic such as by MUXing in the final sum values to the existing adders as needed, soas to perform theactual final summingcalculation. In suchan instance, the final logic, including theM logic 15a, I logic 15b, and D logic 15c blocks, which blocks together form part of the HMMMID instance 17,may include 7multipliers and 4 adders along with the various MUXing involved. 42 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55

[0259] Accordingly, FIG. 10 sets forth the M, I, and D state update circuits 15a’, 15b’, and 15c’ including the effects of simplifying assumptions related to transition probabilities, aswell as the effect of sharing variousM, I, and / or D resources, e.g., adder resources, for the final sum operations. A delay block may also be added to the M-state path in the M-state computation block, as shown in FIG. 10. This delay may be added to compensate for delays in the actual hardware implementations of the multiply and addition operations, and / or to simplify the control logic, e.g., 15.

[0260] As shown in FIGS. 9 and 10, these respective multipliers and / or adders may be floating point multipliers and adders. However, in various instances, as can be seen with respect to FIG. 11, a log domain configuration may be implemented where in such configuration all of themultiplies turn into adds. FIG. 11 presents what log domain calculation would look like if all the multipliers turned into adders, e.g., 15a", 15b", and 15c", such as occurs when employing a log domain computational configuration. Particularly, all of themultiplier logic turns into an adder, but the adder itself turns into or otherwise includes a function where the function such as: f(a,b) = max(a,b) - log2(1+2^(‑[a-b]), such as where the log portion of the equation may be maintained within a LUT whose depth and physical size is determined by the precision required.

[0261] Given the typical read and haplotype sequence lengths as well as the values typically seen for read quality (Phred) scores and for the related transition probabilities, the dynamic range requirements on the internal HMM state values may be quite severe. For instance, when implementing the HMMmodule in software, various of the HMM jobs 20 may result in underruns, such aswhen implemented on single-precision (32-bit) floating-point state values. This implies a dynamic range that is greater than 80 powers of 10, thereby requiring the variant call software to bump up to double- precision (64-bit) floating point state values. However, full 64-bit double-precision floating-point representation may, in various instances, have somenegative implications, such as if compact, high-speed hardware is to be implemented, both storage and compute pipeline resource requirements will need to be increased, thereby occupying greater chip space, and / or slowing timing. In such instances, a fixed-point-only linear-domain number representation may be implemented. Nevertheless, the dynamic rangedemands on the state values, in this embodiment,make thebit widths involved in certain circumstances less than desirable. Accordingly, in such instances, fixed-point-only log-domain number representation may be implemented, as described herein.

[0262] In suchascheme,ascanbeseenwith respect toFIG.11, insteadof representing theactual statevalue inmemory and computations, the -log-base‑2 of the number may be represented. This may have several advantages, including employing multiply operations in linear space that translate into add operations in log space; and / or this log domain representationof numbers inherently supportswiderdynamic rangewithonly small increases in thenumberof integer bits. These log-domain M, I, D state update calculations are set forth in FIGS. 11 and 12.

[0263] Ascanbeseenwhencomparing the logic 17 configurationofFIG. 11with that of FIG.9, themultiplyoperationsgo away in the log-domain. Rather, they are replaced by add operations, and the add operations aremorphed into a function that can be expressed as a max operation followed by a correction factor addition, e.g., via a LUT, where the correction factor is a function of the difference between the two values being summed in the log-domain. Such a correction factor can be either computed or generated from the look-up-table. Whether a correction factor computation or look-up-table implementation is more efficient to be used depends on the required precision (bit width) on the difference between the sum values. In particular instances, therefore, the number of log-domain bits for state representation can be in the neighborhood of 8 to 12 integer bits plus 6 to 24 fractional bits, depending on the level of quality desired for any given implementation.This implies somewherebetween14and36bits total for log-domainstatevalue representation.Further, it has been determined that there are log-domain fixed-point representations that can provide acceptable quality and acceptable hardware size and speed.

[0264] In various instances, one read sequence is typically processed for each HMM job 20, which as indicated may include a comparison against one or two haplotype sequences, ormore. And like above for the haplotypememory, a ping- pong structure may also be used in the read sequence memory 18 to allow various software implemented functions the ability towritenewHMMjob information20bwhileacurrent job20a isstill beingprocessedby theHMMengine instance13. Hence, a read sequence storage requirement may be for a single 1024x32 two-port memory (such as one port for write, one port for read, and / or separate clocks for write and read ports).

[0265] Particularly, as described above, in various instances, the architecture employed by the system 1 is configured such that in determining whether a given base in a sequenced sample genome matches that of a corresponding base in one or more reference genomes, a virtual matrix is formed, wherein the reference genome is theoretically set across a horizontal axis, while the sequenced reads, representing the sample genome, is theoretically set in descending fashion down the vertical axis. Consequently, in performing an HMM calculation, the HMM processing engine 13, as herein described, is configured to traverse this virtual HMMmatrix. Such processing can be depicted as in FIG. 7, as a swath 35 moving diagonally down and across the virtual array performing the various HMM calculations for each cell of the virtual array, as seen in FIG. 8.

[0266] Moreparticularly, this theoretical traversal involvesprocessingafirst groupingof rowsof cells 35a from thematrix 30 in its entirety, suchas for all haplotypeand readbaseswithin thegrouping, before proceedingdown to thenext grouping of rows 35b (e.g., the next group of read bases). In such an instance, the M, I, and D state values for the first grouping are 43 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 stored at the bottom edge of that initial grouping of rows so that theseM, I, andD state values can then be used to feed the top rowof thenext grouping (swath) down in thematrix 30. In various instances, thesystem1maybeconfigured toallowup to 1008 length haplotypes and / or reads in theHMMaccelerator 8, and since the numerical representation employsW-bits for each state, this implies a 1008word x W-bit memory for M, I, and D state storage.

[0267] Accordingly, as indicated, such memory could be either a single-port or double-port memory. Additionally, a cluster-level, scratchpadmemory, e.g., for storing the results of theswathboundary,mayalsobeprovided.For instance, in accordancewith thedisclosureabove, thememoriesdiscussedalreadyareconfigured for aper-engine-instance13basis. InparticularHMMimplementations,multipleengine instances13a‑(n+1)maybegrouped intoacluster11 that is servicedby a single connection, e.g., PCIe bus 5, to the PCIe interface 4 andDMA3 via CentCom9.Multiple clusters 11a‑(n+1) can be instantiated so as to more efficiently utilize PCIe bandwidth using the existing CentCom 9 functionality.

[0268] Hence, in a typical configuration, somewhere between 16 and 64 engines 13m are instantiated within a cluster 11n, and one to four clusters might be instantiated in a typical FPGA / ASIC implementation of the HMM 8 (e.g., depending on whether it is a dedicated HMM FPGA image or whether the HMM has to share FPGA real estate with the sequencer / mapper / aligner and / or other modules, as herein disclosed). In particular instances, there may be a small amount ofmemory used at the cluster-level 11 in the HMMhardware. Thismemorymay be used as an elastic First In First Out ("FIFO") to captureoutput data from theHMMengine instances13 in the cluster andpass it on toCentCom9 for further transmittal back to the software of theCPU1000 via theDMA3andPCIe 4. In theory, this FIFO could be very small (on the order of two 32-bit words), as data are typically passed on to CentCom 9 almost immediately after arriving in the FIFO. However, to absorb potential disrupts in the output data path, the size of this FIFO may be made parametrizable. In particular instances, theFIFOmaybeusedwithadepthof 512words.Thus, thecluster-level storage requirementsmaybe a single 512x32 two-port memory (separate read and write ports, same clock domain).

[0269] FIG.12Asets forth thevariousHMMstate transitions17bdepicting the relationshipbetweenGapOpenPenalties (GOP),GapClosePenalties (GCP), and transition probabilities involved in determiningwhether andhowwell a given read sequencematches a particular haplotype sequence. In performing such an analysis, theHMMengine 13 includes at least three logic blocks 17b, such as a logic block for determining amatch state 15a, a logic block for determining an insert state 15b, and a logic block for determining a delete state 15c. These M, I, and D state calculation logic 17 when appropriately configured function efficiently to avoid high-bandwidth bottlenecks, such as of the HMM computational flow. However, once the M, I, D core computation architecture is determined, other system enhancements may also be configured and implemented so as to avoid the development of other bottlenecks within the system.

[0270] Particularly, the system1maybe configured soas tomaximize the process of efficiently feeding information from the computing core 1000 to the variant caller module 2 and back again, so as not to produce other bottlenecks that would limit overall throughput. One such block that feeds the HMM core M, I, D state computation logic 17 is the transition probabilities andpriors calculation block. For instance, as canbe seenwith respect to FIG. 9, each clock cycle employs the presentationof seven transition probabilities andonePrior at the input to theM, I,Dstate computationblock15a.However, after the simplifications that result in the architecture of FIG. 10, only four unique transition probabilities and one Prior are employed for each clock cycle at the input of the M, I, D state computation block. Accordingly, in various instances, these calculations may be simplified and the resulting values generated. Thus, increasing throughput, efficiency, and reducing the possibility of a bottleneck forming at this stage in the process.

[0271] Additionally, as described above, the Priors are values generated via the read quality, e.g., Phred score, of the particular base being investigated and whether, or not, that base matches the hypothesis haplotype base for the current cell being evaluated in the virtual HMM matrix 30. The relationship can be described via the equations bellow: First, the readPhred in questionmaybe expressed as a probability = 10^(‑(readPhred / 10)). Then thePrior can be computed based onwhether the readbasematches the hypothesis haplotypebase: If the readbase andhypothesis haplotype basematch: Prior = 1 - read Phred expressed as a probability. Otherwise: Prior = (read Phred expressed as probability) / 3. The divide- by-threeoperation in this last equation reflects the fact that thereareonly fourpossiblebases (A,C,G,T).Hence, if the read and haplotype base did notmatch, then it must be one of the three remaining possible bases that doesmatch, and each of the three possibilities is modeled as being equally likely.

[0272] Theper-read-basePhredscoresaredelivered to theHMMhardwareaccelerator 8as6-bit values.Theequations to derive the Priors, then, have 64 possible outcomes for the "match" case and an additional 64 possible outcomes for the "don’t match" case. This may be efficiently implemented in the hardware as a 128 word look-up-table, where the address into the look-up-table is a 7-bit quantity formedby concatenating thePhred valuewith a single bit that indicateswhether, or not, the read base matches the hypothesis haplotype base.

[0273] Further, with respect to determining the match to insert and / or match to delete probabilities, in various implementations of the architecture for the HMM hardware accelerator 8, separate gap open penalties (GOP) can be specified for the Match-to-Insert state transition, and the Match-to-Delete state transition, as indicated above. This equates to the M2I and M2D values in the state transition diagram of FIG. 12A being different. As the GOP values are delivered to the HMM hardware accelerator 8 as 6-bit Phred-like values, the gap open transition probabilities can be computed in accordance with the following equations: M2I transition probability = 10^(‑(read GOP(I) / 10)) and M2D 44 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 transitionprobability=10^(‑(readGOP(D) / 10)).Similar to thePriorsderivation inhardware, asimple64word look-up-table can be used to derive the M2I and M2D values. If GOP(I) and GOP(D) are inputted to the HMM hardware 8 as potentially different values, then two such look-up-tables (or one resource-shared look-up-table, potentially clocked at twice the frequency of the rest of the circuit) may be utilized.

[0274] Furthermore,with respect to determiningmatch tomatch transitionprobabilities, in various instances, thematch- to-match transition probability may be calculated as: M2M transition probability = 1 - (M2I transition probability + M2D transition probability). If theM2I andM2D transition probabilities can be configured to be less than or equal to a value of½, then in various embodiments the equation above can be implemented in hardware in a manner so as to increase overall efficiency and throughput, such as by reworking the equation to be: M2M transition probability = (0.5 - M2I transition probability) + (0.5 - M2D transition probability). This rewriting of the equation allows M2M to be derived using two 64 element look-up-tables followed by an adder, where the look-up-tables store the results.

[0275] Further still, with respect to determining the Insert to Insert and / or Delete to Delete transition probabilities, the I2I and D2D transition probabilities are functions of the gap continuation probability (GCP) values inputted to the HMM hardware accelerator 8. In various instances, theseGCP valuesmay be 6-bit Phred-like values given on a per-read-base basis. The I2I andD2D valuesmay then be derived as shown: I2I transition probability = 10^(‑(readGCP(I) / 10)), andD2D transition probability = 10^(‑(read GCP(D) / 10)). Similar to some of the other transition probabilities discussed above, the I2I and D2D values may be efficiently implemented in hardware, and may include two look-up-tables (or one resource- shared look-up-table), such as having the same form and contents as the Match-to-Indel look-up-tables discussed previously. That is, each look-up-table may have 64 words.

[0276] Additionally, with respect to determining the Inset and / or Delete to Match probabilities, the I2M and D2M transition probabilities are functions of the gap continuation probability (GCP) values and may be computed as: I2M transition probability = 1 - I2I transition probability, andD2M transition probability = 1 -D2D transition probability, where the I2I and D2D transition probabilities may be derived as discussed above. A simple subtract operation to implement the equations abovemay bemore expensive in hardware resources than simply implementing another 64word look-up-table and using two copies of it to implement the I2M and D2M derivations. In such instances, each look-up-table may have 64 words. Of course, in all relevant embodiments, simple or complex subtract operations may be formed with the suitably configured hardware.

[0277] FIG. 13 provides the circuitry 17a for a simplified calculation for HMM transition probabilities and Priors, as described above, which supports the general state transition diagramof FIG. 12A. As can be seenwith respect to FIG. 13, in various instances, a simple HMM hardware accelerator architecture 17a is presented, which accelerator may be configured to include separateGOPvalues for Insert andDelete transitions, and / or theremay be separateGCPvalues for Insert and Delete transitions. In such an instance, the cost of generating the seven unique transition probabilities and one Prior each clock cyclemaybe configuredas set forth below: eight 64word look-up-tables, one128word look-up-table, and one adder.

[0278] Further, in various instances, the hardware 2, as presented herein, may be configured so as to fit asmany HMM engine instances 13 as possible onto the given chip target (such as on anFPGA, sASIC, or ASIC). In such an instance, the cost to implement the transition probabilities and priors generation logic 17a can be substantially reduced relative to the costs as provided by the below configurations. Firstly, rather than supporting a more general version of the state transitions, such as set forth in FIG. 13, e.g., where there may be separate values for GOP(I) and GOP(D), rather, in various instances, it may be assumed that theGOP values for insert and delete transitions are the same for a given base. This results in several simplifications to the hardware, as indicated above.

[0279] In such instances, only one 64 word look-up-table may be employed so as to generate a single M2Indel value, replacing both the M2I and M2D transition probability values, whereas two tables are typically employed in the more general case. Likewise, only one 64 word look-up-table may be used to generate the M2M transition probability value, whereas two tables and an add may typically be employed in the general case, as M2M may now be calculated as 1‑2xM2Indel.

[0280] Secondly, the assumptionmay bemade that the sequencer-dependent GCP value for both insert and delete are the same AND that this value does not change over the course of an HMM job 20. This means that: a single Indel2Indel transition probabilitymaybecalculated insteadof separate I2I andD2Dvalues, usingone64word look-up-table insteadof two tables; and single Indel2Match transition probabilitymaybecalculated insteadof separate I2MandD2Mvalues, using one 64 word look-up-table instead of two tables.

[0281] Additionally, a further simplifying assumption can bemade that assumes the Inset2Insert andDelete2Delete (I2I and D2D) and Insert2Match and Delete2Match (I2M and D2M) values are not only identical between insert and delete transitions, but may be static for the particular HMM job 20. Thus, the four look-up-tables associated in the more general architecturewith I2I,D2D, I2M,andD2Mtransitionprobabilitiescanbeeliminatedaltogether. In variousof these instances, thestatic Indel2Indel and Indel2Matchprobabilities couldbemade tobeenteredvia softwareor viaanRTLparameter (and so would be bitstream programmable in an FPGA). In certain instances, these values may be made bitstream-program- mable, and in certain instances, a trainingmodemaybe implementedemployinga trainingsequencesoas to further refine 45 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 transition probability accuracy for a given sequencer run or genome analysis.

[0282] FIG. 14 sets forth what the new state transition 17b diagram may look like when implementing these various simplifying assumptions. Specifically, FIG. 14 sets forth the simplified HMM state transition diagram depicting the relationship between GOP, GCP, and transition probabilities with the simplifications set forth above.

[0283] Likewise, FIG. 15 sets forth the circuitry 17a,b for the HMM transition probabilities and priors generation, which supports the simplified state transition diagramof FIG. 14. As seenwith respect to FIG. 15, a circuit realization of that state transition diagram is provided. Thus, in various instances, for the HMMhardware accelerator 8, the cost of generating the transition probabilities and onePrior each clock cycle reduces to: Two64word look-up-tables, andOne128word look-up- table.

[0284] Accordingly, as can be seen with reference to the above discussion as well as FIGS 12B - 12D, one of the challenges in variant calling is distinguishing indel errors from true variants. To do so, a variant callermay be configured to employaHiddenMarkovModel (HMM), asdisclosedherein,whichmodels thestatistical behavior of indel errors, aspart of the probability calculation.As canbe seenwith respect to FIG. 12B, theHMMmayhave input parametersGOPins,GCPins, GOPdel, GCPdel,where GOP and GCP stand for the Gap Open Penalty and Gap Continuation Penalty, respectively, and the subscripts indicate insertion anddeletion. FIG. 12B, illustrates that theHMMparametersmaydependon the context of the readand / or the haplotype being processed, this is because indel errors aremore likely in the presence of short tandem repeats (STRs), and in such an instance, the error probability may depend on both the period and the length of the STR. The error process may differ significantly from one dataset to another, depending on factors such as PCR amplification, and / or other sources of error. For accurate detection, it is useful to use HMM parameters that accurately model the error process. However, where the variant caller is configured to use fixed parameters or predetermined functions, this may fail to accurately model the error process, resulting in poor detection performance.

[0285] Accordingly, in such an instance, such errors may be corrected for, such as through an auto-calibration process disclosed herein. Particularly, presented herein is an HMM Auto-Calibration addresses such problems, for instance, by estimating the PCR parameters directly from the dataset being processed. This operation may be performed after mapping & alignment and prior to variant calling, with or without knowledge of the ground truth and with or without using external databases of known mutations. In such an instance, the parameters depend on both the STR period and the repeat length.

[0286] For a givenSTRperiod and length, a set of N loci with the desired period and length, the pileups of readsmapped to those loci may be examined, counting the indels observed at each locus to estimate the parameters of interest. Particularly, the HMM parameters to be estimated include one or more ofGOPins, GCPins, GOPdel,GCPdel as well as the variant probabilities and , which represent the probability of an indel variant of length l, where positive values of l indicate insertions of l bases and negative values indicate deletions of |l| bases, and the superscript indicates whether the variant is heterozygous or homozygous. In various instances, it may be assumed that the underlying organism is diploid, but it is noted that this can be generalized to non-diploid organisms. Note also that for a single locus with limited coverage depth, it is often difficult to determinewhether the indels are due to errors or a true variant, and such pileupsmay not be particularly helpful for estimating the HMM parameters.

[0287] For example, as can be seen with respect to FIG. 12C, a pileup is presented wherein 11 out of 38 reads contain deletions.Specifically, FIG. 12Cpresents anSTR locuswithmultiple deletions in thepileup. In this instance, theSTRhasa periodof 1baseanda lengthof14bases. It is difficult to determine from thispileupalonewhether thesedeletionsareerrors or evidence of a true variant. However, by considering a sufficient number of loci, it’s possible to accurately estimate the parameters of interest. Thismay be done by finding the parameters thatmaximize the probability of producing the set of N observed pileups. Even pileups that seem completely unhelpful in isolation, in this instance, can play an important role when analyzed in conjunction with other pileups.

[0288] A straightforward way to estimate the parameters of interest is to use the HMM module to calculate the joint probability of theobservedpileups, sweeping theHMMparameters and choosing those thatmaximize the total probability. However, the computational complexity of doing so may be prohibitive, both because of the complexity of the HMM operation and because of the number of independent parameters to sweep. Accordingly, presented herein is a simplified method based on counting the number of indels of each length at each locus, without need using HMM. In such an instance, a qualifying read may be defined as one with a high-confidence alignment spanning the STR with a minimum number of flanking bases on each side. Accordingly, the appropriate calculations may be set forth as follows.

[0289] Let kl,i be the number of qualifying reads containing an indel of length l bases (relative to the reference) aligned at locus i, wherepositive valuesof l indicate insertionsandnegative values indicate deletions, and l=0 indicates theabsence of an indel. LetΨ be an approximation of the probability of making the observations (ni, kl,i),i = 1□N given the parameters GOPins, GCPins, GOPdel, GCPdel, and : 46 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 where: and λ is the STR length measured in bases. In general, our HMM auto-calibration procedure consists of tabulating the values of kl,i and then finding the values of GOPins, GCPins, GOPdel , GCPdel, and that maximize Ψ. This operation is performed for each STR period and length.

[0290] In practice, the number of independent parameters above can be problematic, both because there may be insufficient data to train a large number of parameters, and because searching over a large number of dimensions can be difficult or impractical. Fortunately, it is easy to reduce the number of independent parameters and still get good performance.

[0291] In one embodiment, the following assumptions may be made: This reduces the number of independent variables to 4.

[0292] In another embodiment, these calculations may further be simplified by disregarding the length of the indels. In this embodiment, ki represents the number of qualifying reads with an indel (of any length) aligned at locus i , and ni indicates the total number of qualifying reads at locus i. It may be assumed thatGCP is user-specified (by default,GCP = 10 / ω , where ω is the period of the STR), and αhet and αhom indicate the probability of indel variants of any non-zero length. The calculation may then be defined as: where ω is the period of the STR. 47 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 This reduces thenumber of independent variable to 2,which canbe easily performedby exhaustive search. It is noted that the expression forΨmaybeanapproximate expression that discounts or ignores the possibility that a locusmay contain a mixture of indel variants and indel errors (which may cancel an indel variant). This approximation may be employed in instances where it has little impact on the accuracy of the result.

[0293] In various instances, STRswith a period ranging from 1 to 8 and lengths ranging from 1 to 20 whole periodsmay be considered. In such an instance, each STR in the genomemay be classified according to the period for which it has the greatest repeat length, breaking ties toward shorter periods. A target quantity of 2K to 4K STR loci of each period / length combination may be sampled pseudo-randomly from the genomic regions covered by the aligned reads.

[0294] When fewer than 4K STR loci are available in a given period / length class, all covered STRsmay be considered, even though this quantity is much smaller than 2K for combinations of long period and high repeat length. In such an instance, each STR period / length class failing to meet aminimum sample count ofN ≥ 50may bemerged with other STR classes (e.g., merging with STRs with the same period but smaller repeat length) prior to maximum-likelihood parameter estimation. For each period and repeat length, a maximum-likelihood parameter estimation may be performed as described above, sweeping the parameters GOP and αhet over a 2-dimensional grid of integers on a phred scale. For each period, start with the lowest repeat length, where theGOP should be monotonically non-increasing with increasing repeat length, An increase inGOPmay be an indication of insufficient data. If an increase inGOP is observed, the class may be merged with the previous (shorter repeat-length) class.

[0295] Thismethodof indel errormodel estimation isapplicable todiploid germlineDNA-seq,givenasamplecoveringat the equivalent of human whole-exome (tens of millions of locus nucleotides) at substantial coverage depth (say 10x or deeper). Modification for other ploidy is straightforward. Substantially smaller samples, such as amplicon panels, lack enough STR loci to calibrate themodel across important period / length combinations; but variant calling on small samples could use amodel estimated froma larger dataset with similar PCRand sequencing protocols. Thismethod remains valid for whole-exome or whole-genome tumor samples, because although somatic variants violate the 50% / 100% allele frequency assumptions, there are too few real ones to disturbmodel parameter estimation. It also should be applicable to RNA-seq data, provided a sensitive spliced aligner is employed, and STR loci interrupted by alignment introns may be ignored.

[0296] FIG. 12D shows the indel ROC for a dataset SRA056922 (a human whole genome dataset). It can be seen that this HMM auto-calibration provides a large gain in indel sensitivity. For this dataset, the best f-measure increases from 0.9113 to 0.9319.

[0297] As set forth above, the engine control logic 15 is configured for generating the virtual matrix and / or traversing the matrix so as to reach the edge of the swath, e.g., via high-level engine state machines, where result data may be finally summed, e.g., via final sum control logic 19, and stored, e.g., via put / get logic. Accordingly, as can be seenwith respect to FIG. 16, in variousembodiments, amethod for producingand / or traversinganHMMcellmatrix 30 is provided.Specifically, FIG. 16 sets forth an example of how the HMM accelerator control logic 15 goes about traversing the virtual cells in the HMM matrix. For instance, assuming for exemplary purposes, a 5 clock cycle latency for each multiply and each add operation, theworst-case latency through theM, I, D state update calculationswould be the 20 clock cycles it would take to propagate through the M update calculation. There are half as many operations in the I and D state update calculations, implying a 10 clock cycle latency for those operations.

[0298] These latency implications of the M, I, and D compute operations can be understood with respect to FIG. 16, which sets forth variousexamplesof the cell-to-cell datadependencies. In such instances, theMandDstate informationof agivencell feed theDstatecomputationsof thecell in theHMMmatrix that is immediately to the right (e.g., having thesame read base as the given cell, but having the next haplotype base). Likewise, the M and I state information for the given cell feed the I state computationsof the cell in theHMMmatrix that is immediately below (e.g., having the samehaplotypebase as the give cell, but having the next read base). So, in particular instances, theM, I, and D states of a given cell feed the D and I state computations of cells in the next diagonal of the HMM cell matrix, as described above.

[0299] Similarly, the M, I, and D states of a given cell feed the M state computation of the cell that is to the right one and downone (e.g., havingboth thenext haplotypebaseANDthenext readbase). This cell is actually twodiagonalsaway from thecell that feeds it (whereas, the I andDstate calculations relyon states fromacell that is onediagonal away). This quality of the I andDstatecalculations relyingoncells onediagonal away,while theMstatecalculations relyoncells twodiagonals away, has a beneficial result for hardware design. 48 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55

[0300] Particularly, given these configurations, I and D state calculations may be adapted to take half as long (e.g., 10 cycles) as theMstate calculations (e.g., 20 cycles).Hence, ifMstate calculations are started10 cycles before I andDstate calculations for the samecell, then theM, I, andDstate computations for a cell in theHMMmatrix 30will all complete at the same time. Additionally, if thematrix 30 is traversed in a diagonal fashion, such ashaving a swath 35of about 10 cells each within it (e.g., that spans ten readbases), then:TheMandDstatesproducedbyagivencell at (hap, rd) coordinates (i, j) can be used by cell (i+1, j) D state calculations as soon as they are all the way through the compute pipeline of the cell at (i, j).

[0301] The M and I states produced by a given cell at (hap, rd) coordinates (i, j) can be used by cell (i, j+1) I state calculations one clock cycle after they are all theway through the compute pipeline of the cell at (i, j). Likewise, theM, I and D states produced by a given cell at (hap, rd) coordinates (i, j) can be used by cell (i+1, j+1) M state calculations one clock cycle after they are all the way through the compute pipeline of the cell at (i, j). Taken together, the above points establish that very little dedicated storage is needed for the M, I, and D states along the diagonal of the swath path that spans the swath length, e.g., of ten reads. In suchan instance, just the registers required to delay cell (i, j)M, I, andDstate valuesone clock cycle for use in cell (i+1, j+1) M calculations and cell (i, j+1) I calculations by one clock cycle). Moreover, there is somewhat of a virtuous cycle here as theMstate computations for a given cell are begun 10 clock cycles before the I andD state calculations for that same cell, natively outputting the new M, I, and D states for any given cell simultaneously.

[0302] In view of the above, and as can be seen with respect to FIG. 16, the HMM accelerator control logic 15 may be configured to process the data within each of the cells of the virtual matrix 30 in a manner so as to traverse the matrix. Particularly, in various embodiments, operations start at cell (0,0), with M state calculations beginning 10 clock cycles before I and D state calculations begin. The next cell to traverse should be cell (1,0). However, there is a ten cycle latency after thestart of I andDcalculationsbefore the results fromcell (0,0)will beavailable. Thehardware, therefore, inserts nine "dead" cycles into the compute pipeline. These are shown as the cells with haplotype index less than zero in FIG. 16.

[0303] After completing the dead cycle that has an effective cell position in the matrix of (‑9,‑9), the M, I, and D state values for cell (0,0)areavailable. These (e.g., theMandDstateoutputsof cell (0,0))maynowbeusedstraight away tostart theDstate computationsof cell (0,1).Oneclock cycle later, theM, I, andDstate values fromcell (0,0)maybeused to begin the I state computations of cell (0,1) and the M state computations of cell (1,1).

[0304] The next cell to be traversed may be cell (2,0). However, there is a ten cycle latency after the start of I and D calculations before the results from cell (1,0) will be available. The hardware, therefore, inserts eight dead cycles into the compute pipeline. Theseare shownas the cellswith haplotype index less than zero, as in FIG. 16 along the samediagonal ascells (1,0)and (0,1).After completing thedeadcycle thathasaneffectivecell position in thematrixof (‑8, ‑9), theM, I, and Dstate values for cell (1,0) are available. These (e.g., theMandD state outputs of cell (1,0)) are nowused straight away to start the D state computations of cell (2,0).

[0305] One clock cycle later, theM, I, andD state values from cell (1,0)may be used to begin the I state computations of cell (1,1) and theMstate computations of cell (2,1). TheMandDstate values fromcell (0,1)may then be used at that same time to start theDstate calculationsof cell (1,1).Oneclock cycle later, theM, I, andDstate values fromcell (0,1) areused to begin the I state computations of cell (0,2) and the M state computations of cell (1,2).

[0306] Now, the next cell to traverse may be cell (3,0). However, there is a ten-cycle latency after the start of I and D calculations before the results from cell (2,0) will be available. The hardware, therefore, inserts seven dead cycles into the compute pipeline. These are again shown as the cells with haplotype index less than zero in FIG. 16 along the same diagonal as cells (2,0), (1,1), and (0,2). After completing the dead cycle that has an effective cell position in the matrix of (‑7,‑9), theM, I, and D state values for cell (2,0) are available. These (e.g., theM and D state outputs of cell (2,0)) are now used straight away to start theD state computations of cell (3,0). And, so, computation for another ten cells in the diagonal begins.

[0307] Such processingmay continue until the end of the last full diagonal in the swath 35a, which, in this example (that has a read length of 35 and haplotype length of 14), will occur after the diagonal that begins with the cell at (hap, rd) coordinatesof (13,0) is completed.After thecell (4,9) inFigure16 is traversed, thenext cell to traverseshouldbecell (13,1). However, there is a ten-cycle latency after the start of the I and D calculations before the results from cell (12,1) will be available.

[0308] Thehardwaremaybeconfigured, therefore, to start operationsassociatedwith the first cell in thenext swath35b, such as at coordinates (0, 10). Following the processing of cell (0, 10), then cell (13, 1) can be traversed. The whole diagonal of cells beginning with cell (13, 1) is then traversed until cell (5, 9) is reached. Likewise, after the cell (5, 9) is traversed, thenext cell to traverse shouldbe cell (13, 2).However, as before theremaybea ten-cycle latencyafter the start of I andD calculations before the results from cell (12, 2) will be available. Hence, the hardwaremay be configured to start operations associated with the first cell in the second diagonal of the next swath 35b, such as at coordinates (1, 10), followed by cell (0, 11).

[0309] Following the processing of cell (0, 11), the cell (13, 2) can be traversed, in accordance with the methods disclosed above. The whole diagonal 35 of cells beginning with cell (13,2) is then traversed until cell (6, 9) is reached. Additionally, after the cell (6, 9) is traversed, the next cell to be traversed should be cell (13, 3). However, here again there may be a ten-cycle latency period after the start of the I and D calculations before the results from cell (12, 3) will be 49 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 available. The hardware, therefore,may be configured to start operations associatedwith the first cell in the third diagonal of the next swath 35c, such as at coordinates (2, 10), followed by cells (1, 11) and (0, 12), and likewise.

[0310] This continues as indicated, in accordance with the above until the last cell in the first swath 35a (the cell at (hap, rd) coordinates (13, 9)) is traversed, at which point the logic can be fully dedicated to traversing diagonals in the second swath35b, startingwith thecell at (9, 10). Thepatternoutlinedabove repeats for asmanyswathsof 10 readsasnecessary, until the bottom swath 35c (those cells in this example that are associated with read bases having index 30, or greater) is reached.

[0311] In the bottom swath 35, more dead cells may be inserted, as shown in FIG 16 as cells with read indices greater than 35 and with haplotype indices greater than 13. Additionally, in the final swath 35c, an additional row of cells may effectively beadded.These cells are indicatedat line35 inFIG. 16, and relate to adedicated clock cycle in eachdiagonal of the final swath where the final sum operations are occurring. In these cycles, the M and I states of the cell immediately above are added together, and that result is itself summedwith a running final sum (that is initialized to zero at the left edge of the HMM matrix 30).

[0312] Taking the discussion above as context, and in view of FIG. 16, it is possible to see that, for this example of read length of 35 and haplotype length of 14, there are 102 dead cycles, 14 cycles associatedwith final sumoperations, and 20 cycles of pipeline latency, for a total of 102+14+20 = 146 cycles of overhead. It can also be seen that, for any HMM job 20 with a read length greater than 10, the dead cycles in the upper left corner of FIG. 16 are independent of read length. It can also be seen that the dead cycles at the bottom and bottom right portion of FIG. 16 are dependent on read length, with fewest dead cycles for reads having mod(read length, 10) = 9 and most dead cycles for mod(read length, 10) = 0. It can further be seen that the overhead cycles become smaller as a total percentage of HMMmatrix 30 evaluation cycles as the haplotype lengths increase (biggermatrix, partially fixednumberof overheadcycles) or as the read lengths increase (note: this refers to thepercentageof overheadassociatedwith the final sum row in thematrix being reducedas read length -row- count-increases).Using suchhistogramdata from representativewhole humangenome runs, it has beendetermined that traversing the HMM matrix in the manner described above results in less than 10% overhead for the whole genome processing.

[0313] Furthermethodsmaybeemployed to reduce theamountof overheadcycles including:Havingdedicated logic for the final sumoperations rather than sharing adders with theMandD state calculation logic. This eliminates one row of the HMM matrix 30. Using dead cycles to begin HMM matrix operations for the next HMM job in the queue.

[0314] Each grouping of ten rows of the HMMmatrix 30 constitutes a "swath" 35 in the HMM accelerator function. It is noted that the length of the swathmay be increased or decreased so as tomeet the efficiency and / or throughput demands of thesystem.Hence, theswatch lengthmaybeabout five rowsor less toabout fifty rowsormore, suchasabout ten rows to about forty-five rows, for instance, about fifteen or about twenty rows to about forty rows or about thirty-five rows, including about twenty five rows to about thirty rows of cells in length.

[0315] With theexceptions noted in the section, above, related to harvesting cycles thatwould otherwise bedead cycles at the right edge of the matrix of FIG. 16, the HMM matrix may be processed one swath at a time. As can be seen with respect toFIG. 16, the states of the cells in the bottom rowof each swath 35a feed the state computation logic in the top row of the next swath 35b. Consequently, there may be a need to store (put) and retrieve (get) the state information for those cells in the bottom row, or edge, of each swath.

[0316] The logic to do thismay include one ormore of the following: when theM, I, andD state computations for a cell in the HMMmatrix 30 complete for a cell with mod(read index, 10) = 9, save the result to the M, I, D state storage memory. WhenMand I state computations (e.g.,whereDstate computationsdonot require information fromcells above them in the matrix) for a cell in theHMMmatrix 30begin for a cell withmod(read index, 10) = 0, retrieve thepreviously savedM, I, andD state information from the appropriate place in the M, I, D state storage memory. Note in these instances that M, I, and D state values that feed row 0 (the top row) M and I state calculations in the HMM matrix 30 are simply a predetermined constant value and do not need to be recalled frommemory, as is true for theM andD state values that feed column 0 (the left column) D state calculations.

[0317] As noted above, the HMM accelerator may or may not include a dedicated summing resource in the HMM hardware accelerator such that exist simply for the purpose of the final sum operations. However, in particular instances, as described herein, an additional rowmay be added to the bottom of theHMMmatrix 30, and the clock cycles associated with this extra row may be used for final summing operations. For instance, the sum itself may be achieved by borrowing (e.g., as per FIG. 13) an adder from the M state computation logic to do the M+I operation, and further by borrowing an adder from theD state computation logic to add the newly formedM+I sum to the running final sumaccumulation value. In such an instance, the control logic to activate the final sum operation may kick in whenever the read index that guides the HMM traversing operation is equal to the length of the inputted read sequence for the job. These operations canbe seen at line 34 toward the bottom of the sample HMM matrix 30 of FIG. 16.

[0318] Hence, as can be seen above, in one implementation, the variant caller may make use of the mapper and / or aligner engines to determine the likelihood as to where various reads originated, such as with respect to a given location, e.g., chromosomal location. In such instances, the variant caller may be configured to detect the underlying sequence at 50 EP 4 682 891 A2 5 10 15 20 25 30 35 40 45 50 55 that location, such as independently of other regions not immediately adjacent to it, such as by implementing the HMM operations set forth herein above. This is particularly useful and works well when the region of interest does not resemble any other region of the genome over the span of a single read (or a pair of reads for paired-end sequencing). However, a significant fraction of the human genome does not meet this criterion, which canmake variant calling, e.g., the process of reconstructing a subject’s genome from the reads that an NGS produces, challenging.

[0319] Particularly, though DNA sequencing has improved dramatically, variant calling remains a difficult problem, largely due to the genome’s redundant structure. As disclosed herein, however, the complexities presented by the genome’s redundancy may be overcome, at least in part, from a perspective driven by short read data. More particularly, the devices, systems, andmethods of employing the same as disclosed hereinmay be configured in such amanner so as to focus onHomologous orSimilar regions thatmay otherwise have been characterized by low variant calling accuracy. In certain instances, such low variant calling accuracy may stem from difficulties observed in read mapping and alignments with respect to homologous regions that typically may result in very low read MAPQs. Accordingly, presented herein are strategic implementations that accurately call variants (SNPs, Indels, and the like) in homologous regions, such as by jointly considering the information present in these homologous regions.

[0320] For instance, many regions of the genome are homologous, e.g., they have near-identical copies located elsewhere in the genome, e.g., in multiple locations, and as a result, the true source location of a read may be subject to considerable uncertainty. Specifically, if a groupof reads ismappedwith lowconfidence, e.g., due toapparent homology, a typical variant caller may ignore and not process the reads, even though they may contain useful information. In other instances, if a read is mis-mapped (e.g., the primary alignment is not the true source of the read), detection errors may result. More specifically, previously implemented short-read sequencing technologies have been susceptible to these problems, and conventional detection methods often leaves large regions of the genome in the dark.

[0321] In some instances, long-read sequencing can be employed to mitigate these problems, but it typically hasmuch higher cost and / or highererror rates, takes longer, and / or suffers fromother shortcomings.Therefore, in various instances, it may be beneficial to perform a multi-region joint detection operation as herein described. For instance, instead of considering each region in isolation and / or instead of performing and analyzing long read sequencing, multi-region joint detection (MRJD) methodologies may be employed, such as where the MRJD protocol considers multiple, e.g., all, locations fromwhichagroupof readsmayhaveoriginated, andattempts todetect theunderlyingsequences together, e.g., jointly, using all available information, which may be regardless of low or abnormal confidence and / or certainty scores.

[0322] For example, for a diploid organism with statistically uniform coverage, a brute force Bayesian calculation, as describedabove,maybeperformed inavariant call analysis.However, in abrute forceMLRDcomputation, thecomplexity of the calculation grows rapidly with the number of regionsN, and the number of candidate haplotypesK to be considered. Particularly, to consider all combinations of candidate haplotypes, the number of candidate solutions forwhich to calculate probabilities may often times be exponential. For instance, as described in greater detail below, in a brute force implementation, the number of candidate haplotypes includes the number of active positions, which if a graph-assembly technique is used to generate the list of candidate haplotypes in a variant call operation, such as in the building of a De Brujin graph as disclosed herein, then the number of active positions is the number of independent "bubbles" in the graph. Hence, such a brute-force calculation can be prohibitively expensive to implement, and as such brute force Bayesian calculations can be prohibitively complex.

[0323] Accordingly, in one aspect, as set forth in FIG. 17A, a method to reduce the complexity of such brute force calculations is herein provided. For instance, as disclosed above, though the speed and accuracy of DNA / RNA sequencing has improved dramatically, especially with respect to the methods disclosed herein, variant calling, e.g., the process of reconstructing a subject’s genome from the reads a sequencer produces, remains a difficult problem, largely due to the genome’s redundant structure. The devices, systems, and methods disclosed herein therefore are configured to reduce thecomplexitiespresentedby thegenome’s redundancy fromaperspectivedrivenbyshort readdata in contrast to long read sequencing. In particular, provided herein aremethods for performing very long read detection that accounts for homologous and / or similar regions of the genome that are usually characterized by low variant calling accuracy without necessarily having to perform long read sequencing.

[0324] For instance, in one embodiment, a system and method for performing multi region joint detection is provided. Specifically, in a first instance, a general variant calling operation may be performed such as employing the methods disclosed herein. Particularly, a general varia...

Claims

1. A method for constructing and using a chimeric reference standard specific to an individual subject for improving accuracy of analysis of a plurality of portions of a genome of a subject, the method comprising: determining a homology of a first portion of a genome to each reference standard of a plurality of reference standards stored in a database, wherein said genome is based on a biological sample obtained from said subject; selecting a portion of a first reference standard that is a best fit to the first portion of the subject's genome based on the determined homology of the first portion to the plurality of reference standards, wherein the first reference standard is a linear reference standard; determining a homology of a second portion of the genome to each reference standard of the plurality of reference standards; selecting a portion of a second reference standard that is a best fit to the second portion of the subject's genome based on the determined homology of the second portion to the plurality of reference standards, wherein the first portion is a portion of a first chromosome of the subject's DNA and the second portion is a portion of a second chromosome of the subject's DNA, or wherein the first portion is a portion of one strand of the subject's DNA and the second portion is a portion of the other strand of the subject's DNA; producing a chimeric reference standard specific to the individual subject by assembling the selected portions of the first reference standard and the second reference standard into one combined, chimeric reference; and matching the chimeric reference standard to said genome of the subject using a mapping, aligning, and / or variant calling procedure, thereby improving the accuracy of the analysis of the plurality of portions of the genome of the subject.

2. The method of claim 1, wherein the first portion is a portion of a first chromosome of the subject's DNA and the second portion is a portion of a second chromosome of the subject's DNA, and wherein the first chromosome and the second chromosome are homologous chromosomes.

3. The method of claim 1 or claim 2, wherein each reference standard of the plurality of reference standards corresponds to i) one or more genomes from one or more members of a family, or (ii) one or more genomes from one or more individuals having an ancestor in common with each other.

4. The method of any one of claim 1 through 3, wherein producing the chimeric reference comprises producing the chimeric reference genome in graph format.

5. The method of claim 4, wherein variations of the chimeric reference from one or both of the first reference standard and the second reference standard are represented by bubbles.

6. The method of claim 4, wherein variations of the chimeric reference from one or both of the first reference standard and the second reference standard are represented by branches.

7. The method of any one of claim 1 to 6, wherein producing the chimeric reference comprises performing a multi-region joint detection process.

8. The method of any one of claim 1 to 7, wherein the matching comprises a mapping procedure.

9. The method of claim 8, wherein the mapping procedure comprises an alt-aware type mapping procedure.

10. The method of any one of claim 1 to 7, wherein the matching comprises an aligning procedure.

11. The method of any one of claim 1 to 7, wherein the matching comprises a variant calling procedure.

12. The method of claim 11, wherein the matching comprises a Genome Analysis Tool Kit (GATK) haplotype caller procedure, a Hidden Markov Model (HHM) tool procedure, or a De Bruijn Graph function procedure.

13. The method of any one of claims 1 to 12, wherein the first portion and / or the second portion comprise(s) at least 1 or 2 million base pairs.

14. A non-transitory computer-readable medium storing software comprising instructions executable by one or more computers which, upon such execution, cause the one or more computers to perform operations for improving accuracy of analysis of a plurality of portions of a genome of a subject, the operations comprising the method of any one of claims 1 to 13.

15. A system for improving accuracy of nucleic acid sequence analysis of nucleic acid sequences from a subject, the system comprising: one or more computers and one or more storage devices storing instructions that are operable, when executed by the one or more computers, to cause the one or more computers to perform operations comprising the method of any one of claims 1 to 13.

Images