System and method for correlated error event mitigation for variant calling

Abstract

To disclose a method, a system and a device that include a computer program for improving the accuracy of variant calls by taking into account the index of correlation error events.SOLUTION: In one aspect, a method may include actions of accessing a pileup of sequence reads aligned to a first region of a reference genome, obtaining information describing one or more characteristics of each of the plurality of reads of the pileup, providing one or more inputs to a probability model describing the one or more characteristics of the plurality of reads of the pileup, wherein the probability model is configured to determine a score, for each hypothesis of one or more hypotheses selected based on the one or more inputs, that indicates whether each hypothesis is true, obtaining output information for each of the one or more hypotheses, and determining, based on the obtained output information, a likelihood that a true variant exists at the first position.SELECTED DRAWING: Figure 1

Description

[Technical field]

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS This application claims the benefit of U.S. Provisional Patent Application No. 62 / 710,348, filed February 16, 2018, entitled "Methods, Devices, and Systems for Performing Foreign Read Detection and Burst Error Detection," the entire disclosure of which is incorporated herein by reference. Sequence Listing This specification contains a Sequence Listing that has been submitted electronically as an ASCII text file named Sequence_Listing.txt. The ASCII text file, created on March 10, 2022, is 877 bytes in size. The material in the ASCII text file is incorporated herein by reference in its entirety. [Background technology]

[0002] A nucleic acid sequencer is an instrument configured to automate the process of sequencing nucleic acids, such as deoxyribonucleic acid (DNA) or ribonucleic acid (RNA). Nucleic acid sequencing is the process of determining the order of nucleotides within a genetic sequence.

[0003] A nucleic acid sequencing device is configured to receive a nucleic acid sample and generate one or more output data, called "reads," that represent the order of nucleotides in the nucleic acid sample. The nucleotides in a DNA sample can include one or more of guanine (G), cytosine (C), adenine (A), and thymine (T), in any combination. The nucleotides in an RNA sample can include one or more of G, C, A, and uracil (U), in any combination.

[0004] The DNA reads generated by the DNA sequencing device can be mapped and aligned to the known DNA sequence of a reference genome. Once mapped and aligned to the reference genome, the sequences of the mapped and aligned reads can be analyzed in terms of the reference genome to identify potential variations that exist between the sequences of the mapped and aligned reads and the reference genome. Summary of the Invention [Means for solving the problem]

[0005] Aspects of the present disclosure are directed to methods, systems, and apparatus, including computer programs, for taking into account correlated error events that are problematic for variant callers used to determine whether alternative (hereinafter "alt") alleles identified in pileups of mapped and aligned sequence reads are true variants.

[0006] According to one innovative aspect of the present disclosure, a method is disclosed for improving the accuracy of variant calling by considering a metric of correlated error events. In some implementations, the method includes accessing, by one or more computers and from one or more memory devices, a pileup of a plurality of sequence reads aligned to a first region of a reference genome; obtaining, by the one or more computers, information describing one or more characteristics of each of the plurality of reads of the pileup corresponding to a first location of the reference genome; and providing, by the one or more computers and based on the obtained information, one or more inputs to a probabilistic model describing one or more characteristics of the plurality of reads of the pileup, wherein the probabilistic model generates, for each hypothesis of the one or more hypotheses selected based on the one or more inputs, a sequence indicating whether the hypothesis is true. The method may include actions of a method configured to determine a core, the method including: providing, by one or more computers, obtaining output information for each of the one or more hypotheses, the output information for each of the one or more hypotheses (i) generated by a probabilistic model based on processing of one or more inputs to a probabilistic model describing one or more characteristics of a respective reading of the pileup, and (ii) indicating a score indicating whether the hypothesis is true; and determining, by the one or more computers and based on the obtained output information generated by the probabilistic model for each of the plurality of hypotheses, a likelihood that a true variant is present at the first location.

[0007] Other variations include corresponding systems, apparatus, and computer programs for performing the actions of the methods defined by instructions encoded in a computer readable storage device.

[0008] These and other variations may optionally include one or more of the following features: For example, in some implementations, determining, based on the obtained output information generated by the one or more computers and by the probabilistic model for each of the multiple hypotheses, the likelihood that a true variant is present at the first position may include: determining, by the one or more computers, based on the output information generated by the probabilistic model for each of the multiple hypotheses, a population score indicating the likelihood that a true variant is present; determining, by the one or more computers, whether the score generated by the population score meets a predetermined threshold; and adding, by the one or more computers, information indicating that a true variant is present at the first position based on determining that the population score meets the predetermined threshold to the VCF file.

[0009] In some implementations, the information indicating that a true variant is present at a first position can include information identifying (i) the first position, (ii) a candidate alternative allele at the first position, and (iii) a population score.

[0010] In some implementations, the information describing one or more characteristics of each read may include information describing (i) a mapping quality score for each sequence read of the pileup at the first position, and (ii) a read allele score for each sequence read of the pileup at the first position for each candidate allele at the first position.

[0011] In some implementations, the read allele score of each read in the pileup at the first position is a score for a particular candidate allele G m,φ Given a reading r i is based on the output generated by the P-HMM model for each of the readings at the first position, which indicates the probability of observing

[0012] In some implementations, the output information can include first output information regarding a first hypothesis of the one or more hypotheses that includes a probability that the sequence read at the first position indicates the occurrence of a homozygous reference with an exogenous allele that matches the alternative, and second output information regarding a second hypothesis of the one or more hypotheses that includes a probability that the sequence read at the first position indicates the occurrence of a homozygous alternative with an exogenous allele that matches the reference allele.

[0013] In some implementations, the information describing one or more characteristics of each sequence read can include information describing (i) a read orientation of each sequence read of the pile-up at the first position, (ii) a position of each base at the first position within each sequence read of the pile-up at the first position referenced to the 5' end of the sequence read, (iii) a read allele score of each sequence read of the multiple sequence reads for each candidate allele at the reference position, and (iv) a base quality score of each read of the base at the first position.

[0014] In some implementations, the read allele score of each sequence read of the pileup at the first position is a score for a particular candidate allele G m,φ Given an array read, r i is based on the output generated by the P-HMM model for each of the sequence reads at the first position, which indicates the probability of observing

[0015] In some implementations, the information describing one or more characteristics of each sequence read can include information describing (i) a read orientation of each sequence read of the pile-up at the first position, (ii) a position of each base at a first position within each sequence read of the pile-up at the first position referenced to the 5' end of the sequence read, and (iii) a read allele score of each sequence read of the multiple sequence reads for each candidate allele at the reference position.

[0016] In some implementations, the output information may include first output information regarding a first hypothesis of the one or more hypotheses, which may include a probability that the sequence read at the first position indicates the occurrence of a homozygous reference with a sequencing error matching the alternative allele, and second output information regarding a second hypothesis of the one or more hypotheses, which includes a probability that the sequence read at the first position indicates the occurrence of a homozygous alternative with a sequencing error matching the reference allele.

[0017] In some implementations, the information describing one or more characteristics of each sequence read can include information describing: (i) a read orientation of each sequence read of the pile-up at the first position; (ii) a position of each base at the first position, such as position "0" 142 within each sequence read, referenced to the 5' end of the sequence read; (iii) a mapping quality score of each sequence read of the pile-up at the first position; (iv) a read allele score of each sequence read of the multiple reads of each candidate allele at the reference position; and (v) a base quality score of each read of the aligned base at the first position.

[0018] In some implementations, the read allele score of each sequence read of the pileup at the first position is a score for a particular candidate allele G m,φ Given a reading r i is based on the output generated by the P-HMM model for each of the sequence reads at the first position, which indicates the probability of observing

[0019] In some implementations, the output information may include a first output information for a first hypothesis of the one or more hypotheses including a first probability that the sequence read at the first position indicates occurrence of a homozygous reference with an exogenous allele matching the alternative, a second probability that the sequence read at the first position indicates occurrence of a homozygous alternative with an exogenous allele matching the reference allele, a third output information for a first hypothesis of the one or more hypotheses including a probability that the sequence read at the first position indicates occurrence of a homozygous reference with a sequencing error matching the alternative allele, and a fourth output information for a second hypothesis of the one or more hypotheses including a probability that the sequence read at the first position indicates occurrence of a homozygous alternative with a sequencing error matching the reference allele.

[0020] In some implementations, the one or more memory devices receive the pile-up of aligned array reads from a field programmable gate array (FPGA) device, where the FPGA includes one or more configurable digital logic gates configured as a mapping and alignment unit for performing read mapping and alignment.

[0021] In some implementations, the computer is configured to access the one or more memory devices using one or more wired or wireless networks, the field programmable gate array (FPGA) device and the one or more memory devices are housed in an expansion card coupled to a circuit board of the sequencing device, the sequencing device is configured to generate sequence reads based on input samples and store the generated sequence reads in the one or more memory devices, and a mapping and alignment unit of the FPGA is configured to access the one or more memory devices to retrieve the generated sequence reads.

[0022] In some implementations, the computer and the sequencing device are each configured to access one or more memory devices using one or more wired or wireless networks, the field programmable gate array (FPGA) device and the one or more memory devices are housed in an expansion card coupled to a circuit board of a server located remotely from the computer and the sequencing device, the sequencing device is configured to generate sequence reads based on input samples and provide the generated sequence reads to the server using the one or more wired or wireless networks for storing the generated sequence reads in one memory device, and a mapping and alignment unit of the FPGA is configured to access the one or more memory devices to retrieve the generated sequence reads.

[0023] These and other aspects of the present disclosure are described in more detail in the following detailed description, which refers to the accompanying drawings. [Brief description of the drawings]

[0024] [Figure 1] FIG. 1 is a context diagram of an example system for correlated error mitigation in variant calling. [Diagram 2] 1 is a flow chart of an example of a process for correlated error mitigation in variant calling. [Diagram 3] 1 is a flow chart of an example process for mapping error mitigation in variant calling. [Figure 4] 11 is a line graph of an example function for converting example output values ​​from the mapping and alignment unit representing first mapping confidence scores to second mapping confidence scores (μ); [Diagram 5] 1 is a flow chart of an example of a process for sequencing error mitigation for variant calling. [Figure 6] 1 is another flow chart of an example process for correlated error mitigation in variant calling. [Figure 7]1 is another flow chart of an example process for correlated error mitigation in variant calling. [Figure 8] FIG. 1 is a block diagram of system components that can be used to implement a system for correlated error mitigation in variant calling. [Figure 9A] FIG. 1 is an example of a summary of experimental results from an implementation of a process for mitigating correlated errors in variant calling performed on a pile-up of sequencing reads, where the results show an example of a true positive result (SEQ ID NO: 1). [Figure 9B] FIG. 9B is an example of a modified set of pileup probability results of FIG. 9A. [Figure 10A] An example of a summary of experimental results from an implementation of a process for mitigating correlated errors in variant calling performed on a pile-up of sequencing reads, where the results show a low probability of occurrence of mapping errors (SEQ ID NO: 2). [Figure 10B] FIG. 10B is an example of a modified set of pileup probability results of FIG. 10A. [Figure 11A] An example of a summary of experimental results from an implementation of a process for mitigating correlated errors in variant calling performed on a pile-up of sequencing reads, where the results indicate a high probability that a candidate alternative allele is unlikely to be a true variant due to sequencing error (SEQ ID NO: 3). [Figure 11B] FIG. 11B is an example of a modified set of pileup probability results of FIG. 11A. [Figure 12A] An example of a summary of experimental results from an implementation of a process for correlated error mitigation in variant calling performed on a pile-up of sequencing reads where the results indicate a high probability that the candidate alternative allele is unlikely to be a true variant due to sequencing errors in both read orientations (SEQ ID NO: 4). [Figure 12B] 13 is an example of a modified set of pileup probability results from FIG. 12. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0025] Aspects of the present disclosure are directed to methods, systems, and apparatus, including computer programs, for taking into account errors that are problematic for variant callers.

[0026] Correlated error events cause traditional variant callers to generate genotyping errors because the internal probability calculations used by traditional variant callers are typically based on the assumption that errors are uncorrelated. There are two phenomena that tend to generate highly correlated errors in variant calling: (1) mapping error, and (2) sequence-specific error. Mapping error occurs when a read maps to a specific location in the reference genome other than the read's true origin. Sequence-specific error, called sequencing error or systematic error, occurs in this specification because certain base sequences tend to generate sequencing errors at a high rate. Both types of error can lead to high-confidence false positives and other genotyping errors in traditional variant callers.

[0027] Some variant callers attempt to mitigate these issues with ad-hoc rules for filtering reads and variants, but such rules do not approach the performance limits possible with more complex algorithms. Other variant callers use machine learning to recognize and suppress such errors, but machine learning has other drawbacks, such as "unstable" or "ambiguous" outputs that cannot be explained and therefore difficult to address scenarios that were not represented in the training data. This disclosure describes new methods for dealing with both types of errors. In particular, rather than relying on ad-hoc rules or a collection of machine learning to account for the correlated nature of these types of errors, this disclosure addresses the possible fact that certain errors are correlated in a pileup of readings into a probability calculation.

[0028] For the purposes of this disclosure, a "correlated error event" is a category of error that refers to two or more mapping errors or two or more sequencing errors. The process described herein may be applied to consider a single type of correlated error event, such as one or more mapping errors or one or more sequencing errors. Alternatively, the process described herein may be applied to consider multiple types of correlated errors, such as one or more mapping errors and one or more sequencing errors.

[0029] Correlated error event mitigation system

[0030] 1 is a context diagram of an example of a system 100 for correlated error mitigation in variant calling. System 100 includes a nucleic acid sequencer 110 and a secondary analysis unit 120.

[0031] The nucleic acid sequencer 110 is configured to perform a primary analysis of the biological sample 105 to convert the raw physical signals detected by the sequencer into an ordered series of nucleotide base calls with associated quality scores. The primary analysis is specific to the nature of the sequencing technology being employed. In some implementations, for example, nucleotides can be detected by sensing changes in fluorescence, charge, current, emitted light, or any combination thereof. In some embodiments, the biological sample includes DNA, RNA, PNA, LNA, chimeric or hybrid forms of nucleic acids.

[0032] The nucleic acid sample may be a purified sample or a raw DNA sample, including, for example, a lysate derived from a buccal swab, paper, fabric, or other substrate impregnated with saliva, blood, or other bodily fluid. In some implementations, the nucleic acid sample may include small or fragmented portions of DNA, such as genomic DNA. In some implementations, the target sequence may be present in one or more bodily fluids, including, but not limited to, blood, plasma, semen, urine, and serum. In some implementations, the target sequence may include nucleic acid obtained from non-human DNA, such as microbial, plant, or entomological DNA.

[0033] The primary analysis may include receiving the biological sample 105 and generating output data 112, referred to as one or more base calls, each having a quality score, assembled into a plurality of "reads," each representing an ordered set of nucleotides in a sequence fragment prepared from the received biological sample 105. In some implementations, the biological sample 105 may include a DNA sample, and the sequencing device 110 may perform the primary analysis to output a plurality of reads including an ordered sequence of nucleotides or bases from the DNA sample. In such implementations, the sequenced nucleotide order includes one or more of guanine (G), cytosine (C), adenine (A), and thymine (T) in any combination. In other implementations, the biological sample 105 may include an RNA sample. In such implementations, the sequenced nucleotide order includes one or more of G, C, A, and uracil (U) in any combination. Thus, although the example of FIG. 1 describes a DNA sequencing device 110 that generates output reads based on an input DNA sample, other implementations may include a sequencing device 110 that generates output reads based on an RNA sample. Depending on the sequencing method used, the reads, including an ordered sequence of contiguous base pairs, sequenced from one or more fragments of the biological sample 105 may vary in length from about 30 base pairs to 10,000 base pairs or more. For example, in some implementations, the read length of the sequenced fragments may be between about 150 base pairs to 500 base pairs, about 150 base pairs, about 250 base pairs, or about 300 base pairs. The reads may be single reads or paired end reads from fragments prepared from the biological sample 105.

[0034] In some implementations, the nucleic acid sequencer 110 includes a next-generation sequencer (NGS) configured to generate sequence reads 112 for a given sample 105 in a manner that achieves ultra-high throughput, scalability, and speed using massively parallel sequencing techniques. In various examples, NGS allows for rapid sequencing of entire genomes, the ability to expand into deeply sequenced target regions, the use of RNA sequencing (RNA-Seq) to discover novel RNA variants and splice sites, or quantification of mRNA for gene expression analysis, analysis of epigenetic factors such as genome-wide DNA methylation and DNA-protein interactions, sequencing of cancer samples to study rare somatic variants and tumor subclones, and the study of microbial diversity in humans or the environment.

[0035] The nucleic acid sequencing device 110 is configured to generate sequence reads 112 and provide the generated sequence reads 112 to a secondary analysis unit 120. The secondary analysis unit 120 may include one or more computers, such as one or more memory devices 122, a field programmable gate array 124, and a variant calling unit 130. The one or more computers may include one or more devices configured to perform one or more operations. The one or more computers may include only hardware, only software, or any combination thereof.

[0036] In some implementations, the secondary analysis unit 120 may be integrated with the nucleic acid sequencing device 110. In such implementations, for example, each of the one or more components of the secondary analysis unit 120 may be housed on an expansion card, such as a peripheral component interconnect (PCI) expansion card, and installed in the nucleic acid sequencing device 110. In other implementations, for example, each of the one or more components of the secondary analysis unit 120 may be part of a separate computer that is distinct from the nucleic acid sequencing device 110 and directly connected to the nucleic acid sequencing device 110 using an Ethernet cable, a USB cable, a USB-C cable, or the like. In still other implementations, for example, each of the components of the secondary analysis unit 120 is integrated into a cloud-based server that is remotely accessible by the nucleic acid sequencing device 110 using one or more wired or wireless networks, such as a local area network (LAN), a wide area network (WAN), a cellular network, the Internet, or a combination thereof. In yet other implementations, for example, one or more components of the secondary analysis unit 120 are integrated into the nucleic acid sequencer 110, and one or more components of the secondary analysis unit 120 are integrated into another computer, such as a cloud-based server. In such implementations, for example, the FPGA 124 used to implement the mapping and alignment unit 126 is integrated into the nucleic acid sequencer 110 and memory 122, and the variant calling unit 130 is integrated into another computer, such as a cloud-based server.

[0037] Each of these components described with reference to Figure 1, including the nucleic acid sequencer 110, the secondary analysis unit 120, and one or more other computers, such as one or more cloud-based servers, may alternatively or additionally be capable of communicating over one or more wired or wireless networks, including one or more of a LAN, a WAN, a cellular network, the Internet, or combinations thereof, if communication over a direct connection is not enabled. Similarly, each of the components of the secondary analysis unit 120 may be configured to communicate with each other or with components external to the secondary analysis unit 120 using one or more buses, one or more direct connections, or one or more networks to achieve interaction between the respective components as described herein.

[0038] The secondary analysis unit 120 is configured to receive the readings 112 and store the readings 112 in a first portion 122a of the memory 122. The field programmable gate array (FPGA) 124 is dynamically configurable to implement one or more modules of a genomic data analysis pipeline. For example, the FPGA 124 can be dynamically configured to implement a mapping and alignment unit 126, a pair of hidden Markov model (P-HMM) units 128, or both. In some implementations, the mapping and alignment unit 126 is a single functional module. In other implementations, the mapping and alignment unit 126 is separated into two separate functional modules including a dedicated mapping unit 126a and a dedicated alignment unit 126b. In some implementations, the FPGA 124 is configured to implement both the mapping and alignment unit 126 and the P-HMM unit 128 at a particular time.

[0039] However, in other implementations, the FPGA 124 may be dynamically reconfigured on demand to implement any of the specific genome analysis modules, or other computers described herein at any time. For example, the FPGA 124 may be initially configured to include a mapping and alignment unit 126, and then, once the mapping and alignment operations are performed by the FPGA 124 on the reads retrieved from the memory 122, the FPGA 124 may be dynamically reconfigured as a P-HMM unit 128. The FPGA 124 may be dynamically reconfigured on demand from one genome analysis module to another genome analysis module, or other computers, as dictated by a particular genome analysis workflow. For purposes of this disclosure, the terms unit and module are used interchangeably to mean one or more hardware components, one or more software components, or any combination thereof, configured to perform one or more specific operations.

[0040] With reference to the FPGA unit 124, implementation of the functional operations of each mapping and alignment unit 126 and P-HMM unit 128 can be accomplished in hardware by programming programmable digital logic gates using a hardware description programming language, such as Very High Speed ​​Integrated Circuit (VHSIC) Hardware Description Language (VHDL), to dynamically configure or reconfigure the programmable logic gates. Alternatively, the FPGA 124 can be used to implement the mapping and alignment unit 126 and P-HMM unit 128, although the disclosure need not be so limited. For example, in other implementations, one or more of the mapping and alignment unit 126 and P-HMM unit 128 may also be implemented on one or more computers, either locally to the DNA sequencing device or remotely from the DNA sequencing device, using software. In still other implementations, the FPGA 124 may also be configured to perform the functions of the variant calling unit 130 to perform an analysis of the modified probability results to determine whether such probability results should be included in the VCF file 170 generated by the variant calling unit 130.

[0041] However, the present disclosure is not limited to the use of the dynamically reconfigurable FPGA 124 to implement one or more of the mapping and alignment unit 126, the P-HMM module 128, or other computers of the secondary analysis unit 120, such as the variant calling unit 130 described herein. Alternatively, other types of programmable or non-programmable integrated circuits can be used. For example, one or more application specific integrated circuits (ASICs) can be programmed to perform the functions of one or more of the respective genomic analysis modules or other computers described herein. ASICs include integrated circuits that include one or more programmable logic circuits, similar to the FPGAs described herein, in that the digital logic gates of the ASICs are programmable using a hardware description language, such as VHDL. However, ASICs differ from FPGAs in that ASICs can only be programmed once and cannot be dynamically reconfigured once programmed. Furthermore, aspects of the present disclosure are not limited to using FPGAs or ASICs to implement the genomic analysis modules or other computers of the secondary analysis unit 120. Alternatively, any of the genomic analysis modules of the secondary analysis unit 120 or other computers may be implemented using one or more central processing units (CPUs), graphical processing units (GPUs), or any combination thereof that implement the genomic analysis module of the secondary analysis unit 120 or other computers through the execution of software instructions.

[0042] In some implementations, the mapping and alignment unit 126 can be implemented using an FPGA 124 configured to map and align the generated reads 112 stored in a first portion 122a of the memory 122 to a reference genome stored in another portion 122b of the memory 122. However, the present disclosure is not limited to storing reads in the memory 122, accessing reads from the memory 122, storing a reference genome in the memory 122, or accessing a reference genome in the memory 122. Instead, in some implementations, the generated reads 112, the reference genome, or both, can be stored in a memory device in a cloud-based server accessible via one or more networks.

[0043] The mapped and aligned reads can be output by the mapping and alignment unit 126 for storage in the memory 122 in a third portion 122c of the memory 122. In some implementations, the write to the memory 122 including the mapped and aligned reads from the FPGA 124, referred to as stored in the third portion 122c, may actually be stored in the memory 122 in a manner that overwrites the original generated reads 112 output by the nucleic acid sequencer 110 and stored in the first portion 122a. Thus, although multiple stages of information are shown as being stored in the first portion 122a, the second portion 122b, and the third portion 122c of the memory 122, respectively, there is no requirement of the present disclosure that all of the data described by the present disclosure as being stored in one of these respective portions of the memory 122 be present in the memory 122 at any particular time during the execution of the processes disclosed by the present disclosure, however, all of the data described herein may be present when stored in the memory 122 at the same time.

[0044] In some implementations, memory 122 may include a single memory device or multiple memory devices. The use of additional memory devices can reduce data access latency and increase throughput by allowing reads and writes to high-speed memory devices, such as flash memory, as opposed to read or write requests to one or more disk storage devices that are accessed using multiple levels of a memory hierarchy used to create the illusion of high-speed memory.

[0045] Similarly, in some implementations, implementing each genome analysis module or other computer of the secondary analysis unit 120 using integrated circuits such as FPGA 124, ASIC, CPU, GPU, or combinations thereof can include a single FPGA 124, a single ASIC, a single CPU, a single GPU, or any combination thereof. Alternatively, or in addition, implementing each genome analysis module or other computer of the secondary analysis unit 120 using integrated circuits such as FPGA 124, ASIC, CPU, GPU, or combinations thereof can include multiple FPGAs 124, multiple ASICs, multiple CPUs, or multiple GPUs, or any combination thereof. Using additional integrated circuits such as multiple FPGAs to implement genome analysis units or other computers of the secondary analysis unit 120 can reduce the time it takes to perform secondary analysis operations such as mapping, alignment, P-HMM probability calculation, and variant calling. In some implementations, implementing these secondary analysis operations using FPGAs can reduce the time it takes to complete these secondary analysis operations from 24 hours or more to 30 minutes or less. In some implementations, using multiple FPGAs to perform these secondary analysis operations can result in these secondary analysis operations being completed in as little as five minutes.

[0046] The output of the mapping and alignment unit 126 includes a pile-up of the reads that are mapped and aligned to the reference genome. The pile-up may include a text-based format for summarizing the base calls of the aligned reads from the DNA sample to the reference genome or a portion of the reference genome. This output may be stored in one or more portions 122b, 122c of the memory 122 in a computer-readable binary format that can be accessed and analyzed by one or more of the FPGA 124, the P-HMM unit 128, and the variant calling unit 130. Alternatively, the output of the mapping and alignment unit 126 may be stored in the memory using a computer-readable binary format on one or more remote computers, such as a remote cloud server accessed using one or more networks. A human-friendly depiction of the output of the mapping and alignment unit 126 of the FPGA 124 may be depicted on a graphical user interface of a user device. An example of such a graphical user interface is shown with reference to the interface 140.

[0047] The interface 140 can be provided for display on a user interface of a user device that can access the memory 122 of the secondary analysis unit 120. For example, in some implementations, the secondary analysis unit 120 has an attached display device. Alternatively, or in addition, other devices with a display device, such as a smartphone or tablet, can connect to the same network as the secondary analysis unit 120, access the memory 122, and then display the pileup, such as the pileup 141 of the interface 140. In such implementations, the variant calling unit 130 can (i) access the output of the FPGA 124 that mapped and aligned the acquired reads 112 to the reference genome stored in the memory, and (ii) generate rendering data that, when rendered on the display device, outputs data representing the reads mapped and aligned to the reference genome by the FPGA 124 and stored in the memory 122 for display on the user device in a human-friendly manner that can be read by a person using the interface 140.

[0048] Interface 140 shows that the output of FPGA 124 includes pileup 141 of mapped and aligned reads. In this example, interface 140 represents 14 reads, each with 14 horizontal lines. The reads are grouped based on the DNA strand from which the reads were generated. For example, pileup 141 of mapped and aligned reads in the example of FIG. 1 includes a first direction extending leftward from the 5'1 end of the read to the 3'1 end of the read. 1 includes a first set of rearward aligned readings and a second set of forward aligned readings extending in a second direction to the right from the 5'2 end of the readings toward the 3'2 end of the readings. Thus, in this example of interface 140 of 14 mapped and aligned readings output from FPGA 124, the bottom seven readings represent the first set of mapped and aligned readings and the top seven readings represent the second set of mapped and aligned readings. The 5' or 3' end of each of the two middle readings, not shown in interface 140, occurs outside window 140. Although this example used to explain the concepts of the present disclosure illustrates a pileup of only 14 readings, the present disclosure is not so limited. Also, although pileup 141 displays rearward aligned readings at the bottom of the pileup and forward aligned readings at the top of the pileup, other alternatives may exist. For example, as shown with reference to the embodiments of Figures 7-10, forward aligned readings can be presented at the bottom of pile-up 141 and rear aligned readings can be presented at the top of the pile-up.

[0049] Although an example of an interface 140 capable of displaying information describing characteristics of a sequence read is provided, there is no requirement that any implementation of the present disclosure output information describing characteristics of a sequencer read on a display device, or that the variant calling unit 130, or other components described herein, access information from interface 140. Instead, interface 140 is provided merely to illustrate an example of the type of information that constitutes characteristics of a sequence read.

[0050] The present disclosure can use the FPGA 124 to generate pileups of mapped and aligned reads including tens of reads, hundreds of reads, thousands of reads, or more, required for a particular genome sequence analysis workflow. As an example, high-throughput next-generation sequencing of a nucleic acid sample can result in hundreds of thousands of short reads that need to be mapped and aligned to one or more regions of a reference genome sequence, or a portion thereof. Mapping and aligning such a large amount of reads can result in a large number of overlapping or overlapping short sequence nucleic acid reads. In some implementations, for example, the number of overlapping or overlapping short sequence reads can include 1x, 5x, 10x, 30x, 100x, or more coverage of one or more respective reference locations of a reference genome sequence. "30x coverage" refers to mapping and alignment of short reads including, for example, a pileup of 30 or more overlapping reads to one or more reference locations of a reference genome sequence. As another example, "5x coverage" refers to the mapping and alignment of short reads, for example, including a pile-up of 5 or more overlapping reads to one or more reference locations of a reference genome sequence.

[0051] Methods for reducing correlated errors

[0052] Therefore, new read processing methods need to be designed to accurately and efficiently map and align large numbers of reads to a reference genome sequence, or a portion thereof. For example, data resulting from sequencing the human genome can result in hundreds of millions of short reads that typically need to be mapped and aligned to locations in the complete reference genome before the short reads can be further analyzed for potential variants to determine their biological, diagnostic, and / or therapeutic relevance.

[0053] The pileup of overlapping reads allows the comparison of each of the different reads at a specific reference location of the reference genome sequence. The analysis of multiple overlapping reads at a specific reference location allows accurate determination of whether there is a true variation, variant, or deviation in the reads that are mapped and aligned to a specific location of the reference genome, or if there is an error in any one of the reads at the problematic location in the pileup. For example, if only one or two reads out of 30 reads detect a specific nucleotide at position "X" of the reference genome sequence, and each of the other 28 or 29 reads support the determination that another nucleotide is present at position "X", then the two outlying reads can be excluded as an error at position "X".

[0054] Analysis of pileups of overlapping reads allows for more accurate analysis of reads compared to methods that do not evaluate the similarity or difference of overlapping reads, and can determine how the subject's genome differs from a reference genome, e.g., a model genome. For example, analysis of pileups of overlapping reads can more accurately identify errors, such as chemical errors, mechanical errors, and read errors, and distinguish such errors from true variants. More specifically, if a subject has a true variant at position "X" in the reference genome, the majority of the reads in the pileup should support the existence of a true variant, e.g., by the majority of the reads that contain the true variant. A statistical model, such as that described herein, can then be implemented to determine the true genetic sequence of the subject with all its true variants from the reference genome.

[0055] Thus, in various examples, once the reads of the nucleic acid sample are generated, their sequence order is aligned, and the generated reads are mapped to a reference genome or a portion thereof, the true genetic sequence of the subject's genome can be determined. Once the true sample genome is determined, one or more true variations can be determined based on a comparison between the true sample genome and the reference genome or a portion thereof. Once one or more variations between the true sample genome and the reference genome or a portion thereof are determined, a list of all true variants or deviations between the sample genome and the reference genome is determined and called. Such variations may be due to various reasons and may have biological, diagnostic, and / or therapeutic relevance.

[0056] The exemplary pileup illustrated by interface 140 shows that the output of the mapping and alignment unit 126 of the FPGA 124 includes a base quality score 143 and a mapping confidence score 144 for each of the pileup reads. The base quality score 143 includes a value indicating a confidence level that the base called for a read at a particular position of interest, such as position "0" 142, is accurate. In the example of FIG. 1, the base quality scores are represented by a range of values ​​defined by a high base quality score of "41", indicating a high level of confidence that the base call for the read at position "0" 142 is accurate, and a low base quality score of "2", indicating a low level of confidence that the base call for the read at position "0" 142 is accurate. In some implementations, the base quality score is an output of the nucleic acid sequencer 110 received by the secondary analysis unit 120 and is expressed as a Phred scaled probability of a base call error Q phred-base can be determined using Q phred-base =-10*log10(P e-base ) In this example, P e-base is the probability of a base calling error for a particular read. In some implementations, a low base quality score may be indicative of a sequencing error.

[0057] The mapping confidence score 144 indicates a confidence level that the obtained read 112 has been correctly mapped to the reference genome 145 by the mapping and alignment unit 126 at a particular position of interest, such as position "0" (indicated by reference number 142). In the example of FIG. 1 , the mapping confidence score is represented by a range of values ​​defined by a high mapping confidence score of "250", indicating a high level of confidence that the read has been correctly mapped to the reference genome 145 at position "0" 142, and a low mapping confidence score of "0", indicating a low level of confidence that the read has been correctly mapped to the reference genome 145 at position "0" 142. In some implementations, the mapping confidence score is an output of the mapping and alignment unit 126 and is a Phred scale Q of the mapping error. phred-mapping can be determined using Q phred-mapping =-10*log10(P e-mapping ) In this example, P e-mapping is the probability of mapping error for a particular read. The value of the mapping confidence score 144 may be proportional to the difference between the best alignment score from an alignment algorithm, such as a Smith-Waterman aligner, and the second best score of the aligner. In some implementations, adjustments to the method may be made to account for the number of secondary alignments. In some implementations, a low mapping score may be a factor indicative of a mapping error.

[0058] Interface 140 also shows that the output of the FPGA includes the base nucleotide called at position "0" 142. In the example of FIG. 1, the top 12 reads of pileup 141 were determined to have the same base call at position "0" 142 as the reference genome. The exemplary interface of 140 represents this determination by not depicting the letters A, C, G, or T, which represent nucleotides, for each of the top 12 reads at position "0" 142. Thus, based on review of the information depicted in interface 140, the top 12 reads can be determined to have a base call of "A" (adenine) at position "0" 142. Interface 140 also shows that an alternative allele of G (guanine) was determined as the base call for the last two reads of the pileup. G (guanine) is an allele because it is different from the nucleotide base of the reference genome at position "0".

[0059] The output of the FPGA also includes additional information that can be determined from an analysis of the pileup 141 shown on the interface 140. First, information describing the strand direction of each read can be determined, also referred to as the read orientation of each read. For example, the interface 140 shows that each of the alternative alleles occurs in the same strand direction or the same read orientation. In this example, such information is evidenced by the fact that an alternative allele (e.g., "G") occurs in the first set of reads in a backward aligned orientation. As another example, the information describing the strand can also include the strand direction of each read of the pileup, which can be determined with reference to their respective 3' and 5' ends. Second, information describing the location of the alternative allele within each read can also be determined. For example, the proximity between the alternative alleles of each read at position "0" 142 from 5' to the end of the read can be determined. In this example, the determined alternative allele "G" occurs far from the 5' end of the respective read because each alternative allele is close to the opposite 3' end. The further from the 5' end the alternative allele occurs, the more likely the alternative allele is associated with a sequencing error. Third, a base quality for each read at position "0" 142 can be determined. As an example, reads with alternative allele "G" have base qualities of 6 and 2, respectively. Fourth, for each read at reference position "0", a mapping confidence score for each read can be determined. As an example, reads with alternative allele "G" at position "0" 142 have mapping confidence scores of 45 and 3, respectively. For purposes of this example, the information shown in interface 140 is merely an example to illustrate features of the present disclosure. However, real world examples of such interfaces are shown with reference to Figures 6-9.

[0060] Each type of information displayed by interface 140 can be referred to as one or more properties of the DNA read. Properties include read-specific properties such as base quality score 143 or mapping confidence score 144. Additional examples of read-specific properties are described with reference to interface 140 below.

[0061] Information displayed by interface 140 is also stored in memory 122 . As an example, an interface used to generate interface 140 may use information describing characteristics of the DNA reads, such as information indicative of a pileup of mapped and aligned reads 141, information indicative of a reference genome 145 (or a portion thereof), information indicative of a base quality score for each read 143, information indicative of a mapping confidence score 144 for each read, information indicative of a base call for each read at a reference position (e.g., position "0" 142), information indicative of whether each read contains an alternative allele, information indicative of the identity of the reads that have alternative alleles, information indicative of the location of each alternative allele with reference to the 5' end of the read that contains the alternative allele, information indicative of the direction (or orientation) of each read (e.g., aligned forward or aligned reverse), information indicative of the direction (or orientation) of each read that has an alternative allele (e.g., aligned forward or aligned reverse), and information indicative of whether it has been determined that the alternative allele occurs in a first set of reads in a forward aligned orientation or in a second set of reads in a reverse aligned orientation to generate interface 140. Information describing, indicating, or otherwise representing each of these characteristics may be stored in memory 122, for example at locations 122b, 122c. By way of example, these characteristics may be stored in memory 122 in a machine-readable binary format.

[0062] The variant calling unit 130 can obtain from the memory 122 information describing the characteristics of the mapped and aligned reads and shown on the interface 140. For some of the inputs, the variant calling unit 130 can use the information describing the characteristics of the mapped and aligned reads from the memory 122 to generate inputs to one or more probability models 131 of the variant calling unit 130. The variant calling unit 130 can provide the generated inputs as inputs to the one or more probability models 131. In some implementations, the variant calling unit 130 can also request that the P-HMM unit 128 generate one or more probabilities for the inputs to the one or more probability models 131 used by the variant calling unit 130. As an example, the variant calling unit 130 can request that the P-HMM unit 128 generate one or more probabilities for the inputs to the one or more probability models 131 used by the variant calling unit 130. m,φ For each read, the variant calling unit 130 may be requested to determine the probability of observing the read ri given r. In such implementations, for example, the variant calling unit 130 may use the probability values ​​returned by the P-HMM unit 128 as read alleles. The variant calling unit 130 may then provide information describing the properties of the mapped and aligned reads, including (i) information describing the properties retrieved from the memory 112 and / or a remote memory and (ii) information probabilities calculated by the P-HMM unit 128 describing one or more properties of the pileup described by the interface 140, or any combination thereof. In some implementations, the variant calling unit 130 may calculate or receive calculation results from another computer of the secondary analysis unit 120, which represent alternative forms of read allele scores. These different forms of read allele scores are described in more detail below.

[0063] It should be noted that memory 122 stores information describing, supporting, or both, one or more characteristics of the read described by interface 140. In some cases, the type of information described with reference to interface 140 may need to be derived from information actually stored in memory 122. For example, in some implementations, memory 122 may store the location of a candidate alternative allele, such as "G" in interface 140, and the location of the 5' end of the read that includes the candidate alternative allele, such as "G" in interface 140. Then, based on that information stored in memory 122, variant calling unit 130, or other components of secondary analysis unit 120, can determine the distance of the candidate alternative allele "G" from the 5' end of the read. In such a case, memory 122 does not need to store the distance of the candidate alternative allele "G" from the 5' end, since such information can be derived from information stored in memory 122. Other types of information describing characteristics of sequence reads can be similarly derived from actually stored information describing the characteristics.

[0064] The probability models 131 described herein are used by the variant calling unit 130 as described herein to determine probability scores of various candidate genotypes being true. The improvements presented by the present disclosure improve the accuracy of the variant caller in ways that traditional probability models cannot account for the occurrence of correlated error events, such as two or more mapping errors and / or two or more sequencing errors. These new probability models 131 include a mapping error probability model 132 and a sequencing error probability model 134. These probability models provide technical advantages because they are not limited by rule-based decision-making or the characteristics of a given training dataset.

[0065] Mapping error probability model

[0066] The mapping error probability model 132 is designed to take into account mapping errors, for example errors that occur when multiple regions of a reference genome, such as a first region and a second region, contain similar and possibly nearly identical base sequences. In such a case, the mapping and alignment unit may erroneously map a set of reads to the first region rather than the second region due to the similarity of the base sequences of the respective regions. The possibility of mapping errors may be exacerbated when there are naturally occurring variants at one or more locations in the first region, which may result in a sequence identical to that of the second region. To take into account such errors, the mapping error model receives as input to the probability model a set of information obtained by the variant calling unit 130 from the memory 122, from the P-HMM unit 128, or a combination thereof, which describes the characteristics of the mapped and aligned reads.

[0067] To determine the probability that a mapping error occurred, the mapping error probability model 132 receives input obtained by the variant calling unit 130, including (i) a read allele score for each read of the pileup stored in memory 122 for each candidate allele at a reference position, such as reference position "0" 142, and (ii) a mapping quality score for each read of the pileup stored in memory 122. In some implementations, the read allele scores include a probability of allele G m,φ Given a DNA molecule containing i The value P(r i |G m,φ ) can be included. This value P(r i |G m,φ There are various means to calculate or estimate .

[0068] In implementations using haplotype-based callers such as the Genome Analysis Toolkit (GATK) or Dragen® platforms, a de Brun graph is used to identify haplotype H kA list including haplotypes H can be generated, where the haplotypes represent sequences of bases that extend in one or both directions beyond a reference position, such as reference position "0" 142. A hidden Markov model (HMM), such as P-HMM unit 128, can then be used to generate a list including the haplotypes H k Given a DNA molecule containing i The read haplotype score P(r i |H k ) can be calculated.

[0069] In the HMM computation, rather than assuming that the alignment returned by the mapper / aligner is correct, the probabilities of multiple possible alignments can be summed, taking into account possible uncertainties in the alignment. Then, in some implementations, read allele scores can be assigned the highest score possible over the haplotype containing the allele.

[0070]

number

[0071] A detailed description of the use of the variant calling unit 130 to calculate such probabilities using an HMM unit is described in more detail in U.S. Patent Publication No. 2016 / 0306922, the entire contents of which are incorporated herein by reference.

[0072] Variant callers other than haplotype-based callers can use simpler estimates to reduce computational complexity. For example, a variant caller can assume that the alignment from the mapper / aligner is correct and estimate the score accordingly. To detect SNPs, such a variant caller would make the base call b aligned to a reference position, such as reference position "0" 142, as follows: i and base quality q i Based on this, a read allele score can be estimated for each read in the pileup stored in memory 122.

[0073]

number

[0074] In the case of indels, such variant callers can assign a score based on the length of the indel, whether the indel is an insertion or deletion, and the surrounding sequence context (e.g., the duration and length of short tandem repeats).

[0075] Regardless of whether a column-wise implementation associated with SNP detection or a more general implementation associated with SNP and indel detection is used, the output of the mapping error probability model 132 includes one or more probabilities, including a score, indicating the likelihood that one or more mapping errors occurred at a reference position, such as position "0" 142. In some implementations, the mapping error probability model 132 is configured to output two probability scores for two different hypotheses, including (i) the likelihood that the read at the reference position 142 indicates the occurrence of a homozygous reference with an exogenous allele matching the substitution, and (ii) the likelihood that the read at the reference position 142 indicates the occurrence of a homozygous substitution with an exogenous allele matching the reference allele (165).

[0076] Sequencing error probability model

[0077] The sequencing error probability model 134 is a probability model that accounts for sequencing errors that may occur because certain combinations of nucleotides may confuse the sequencing algorithm and generate erroneous sequences. Similar to the mapping error model 132 described above, in some implementations, there may be variation in the input that can be provided to the sequencing error probability model 134 based on the complexity of the variant calling unit used. More complex haplotype variant callers may use read allele scores calculated using equation (1) above, while other variant callers other than haplotype variant callers may use read alleles calculated using equation (2). In some implementations, the type of read allele score used may be determined based on whether the variant calling unit 130 is used to detect only SNPs or to detect SNPs and indels.

[0078] Regardless of the type of read allele score used, the sequencing error probability model 134 receives as input from the variant calling unit 130 a set of information retrieved by the variant calling unit 130 that describes characteristics of the mapped and aligned reads. To determine the probability that a sequencing error occurred, the sequencing error probability model 134 receives input generated or otherwise obtained by the variant calling unit 130 including: (i) the read orientation of each read of the pileup stored in memory 122, (ii) the position of each base at a reference position, such as position "0" 142 within each read, that references the 5' end of the read, (iii) the read allele score of each read of the pileup stored in memory 122 for each candidate allele at a reference position, such as reference position "0" 142, and (iv) a base quality score of each read of the base at a reference position, such as position "0".

[0079] Other variations of the sequencing error probability model can be configured to receive other inputs. For example, in some implementations, the base quality score of each read of a base at a reference position, such as position "0" 142, is not required as an input. The base quality score of each read of a base at reference position "0" 142 is an example of a scenario in which the read allele score is determined using Equation (2). In such a case, the base quality score of each read of a base at a reference position, such as reference position "0" 142, can be derived from the read allele score determined using Equation (2). However, there may be other implementations in which the base quality score of each read of a base at a reference position, such as position "0" 142, is not required as a dedicated input, as it can instead be derived from another received input.

[0080] In yet other implementations, another fourth input may be provided to the sequencing error probability model. For example, the sequencing error probability model may be configured to receive an input describing a number of homopolymers in the reference genome 145 preceding a reference location, such as position "0" 142, in the same direction (or read orientation) of a read containing a candidate alternative allele, such as the candidate alternative allele "G" at position "0" of the interface 140 of FIG. 1. Note, for example, that three reference alleles "G" occur in the reference genome 145 prior to the reference location 142 that is the same as the candidate alternative allele "G". This number of homopolymers may be input as another input to the sequencing error probability model. This input describing the number of homopolymers may be added to improve the accuracy of the model.

[0081] The output of the sequencing error probability model 134 includes one or more probabilities, including a score, indicating the likelihood that one or more sequencing errors occurred at a reference position, such as position "0" 142. In some implementations, the sequencing error probability model 134 is configured to output two probability scores for two different hypotheses, including (i) the likelihood that the read at the reference position 142 indicates the occurrence of a homozygous reference with a sequencing error matching the alternative allele (166), and (ii) the likelihood that the read at the reference position 142 indicates the occurrence of a homozygous alternative with a sequencing error matching the reference allele (167).

[0082] The variant calling unit 130 uses the conventional probability model and one or more thereof to generate a set of modified probability outcomes 150 for each of the multiple hypotheses. The multiple hypotheses may be determined based on inputs provided to the one or more probability models 131, the particular one or more probability models 131 used, or a combination thereof. The generated set of modified probability outcomes 150 may include a probability value for each of the multiple hypotheses indicating the likelihood that the respective hypothesis is true.

[0083] As an example, in some implementations, if the variant calling unit 130 generates or otherwise provides only input for the mapping error model 132, then the modified set of probability results 150 includes conventional probabilities for hypotheses 161, 162, 163, and non-conventional probabilities for hypotheses 164, 165, each of which considers the possibility that one or more mapping errors occur at the reference location 142, each of which is described in more detail below. As another example, if the variant calling unit 130 generates or otherwise provides only input for the sequencing error model 134, then the modified set of probability results 150 includes conventional probabilities for hypotheses 161, 162, 163, and non-conventional probabilities for hypotheses 166, 167, each of which considers the possibility that one or more sequencing errors occur at the reference location 142, each of which is described in more detail below. However, if the variant calling unit 130 generates or otherwise provides inputs for both the mapping error model 132 and the sequencing error model 134, then the modified set of probability results 150 generated by the variant calling unit 130 includes conventional probabilities 161, 162, 163, and non-conventional probabilities 164, 164, 166, 167, each of which is described in more detail below. The modified set of probability results 150 may be generated and provided in a computer-readable binary format. The modified set of probability results 150 improves upon the results of the conventional probability calculations typically performed by the variant calling unit 129 in that the modified probability results 150 may also include additional probability scores for one or more additional hypotheses 164, 165, 166, 167 that can be used by the variant calling unit 130 to account for the occurrence of one or more mapping errors, one or more sequencing errors, or a combination of both, using the modified probability calculations.

[0084] The human readable version of the modified set of probability results 150 may be shown using a graphical user interface on a display of a user device that accesses the set of probability results 150 generated by the variant calling unit 130. An example of such a graphical user interface is shown in the example of FIG. 1 with reference to interface 160. In some implementations, the display may include, for example, a display device that is coupled to the secondary analysis unit 120. Each of the probabilities in the modified set of probabilities 150 is described below with reference to the display of the probabilities in interface 160. However, these probabilities may also be obtained and analyzed in a machine readable form by the variant calling unit 130.

[0085] The conventional probabilities calculated by the variant caller 130 may include a set of probabilities determined using conventional probability models. These conventional probability models are configured to determine a probability score for each of three hypotheses including (i) the likelihood that a read at the reference position 142 indicates the occurrence of a homozygous reference 161, (ii) the likelihood that a read at the reference position 142 indicates the occurrence of a heterozygous alternative 162, and (iii) the likelihood that a read at the reference position 142 indicates the occurrence of a homozygous alternative. A homozygous reference occurs when both alleles at the reference position 142 are the same. In such a case, the alternative allele at the reference position 142 does not occur. A heterozygous alternative occurs when one of the alleles at the reference position 142 is the alternative allele and the other allele at the reference position 142 is the reference allele. A homozygous alternative occurs when both alleles at the reference position 142 are the alternative allele. The conventional probability calculations that generate these three hypotheses assume that all reads in the pileup are correctly mapped and that sequencing errors are uncorrelated across the reads.

[0086] However, mapping errors do occur, and mapping errors and sequencing errors tend to be highly correlated between reads. To account for the occurrence of these errors, the present disclosure employs a variant calling unit 130 configured to generate a modified set of probability results 150 using one or more modified probability models 131, including probability scores for four or more additional non-conventional hypotheses. In some implementations, such as a single alternative allele (e.g., a reference allele and a first alternative allele), the modified set of probability results 150 can include probability scores for four additional non-conventional hypotheses 164, 165, 166, 167 described herein. However, in other implementations where there is more than a single alternative allele (e.g., a reference allele, a first alternative allele, and a second alternative allele), then the modified set of probability results 150 can include four or more non-conventional hypotheses described herein. In such a scenario, each combination of the first allele, the second allele, and the reference allele has a generated probability score for a set of non-conventional hypotheses corresponding to the four non-conventional hypotheses 164, 165, 166, 167 described herein. The modified set of probability results 150 provides an improvement to the variant calling unit 130 when compared to a conventional variant caller that does not utilize the probability scores output by the one or more probability models 131 by allowing the variant calling unit 130 to consider the possibility that one or more correlated error events occurred at a reference location, such as reference location 142 of pileup 141.

[0087] The modified set of probability results 150 includes respective non-conventional probability scores for different hypotheses that account for the potential occurrence of one or more mapping errors. These additional probabilities include (i) the likelihood that a read at the reference position 142 indicates the occurrence of a homozygous reference with an exogenous allele matching the alternative (164), and (ii) the likelihood that a read at the reference position 142 indicates the occurrence of a homozygous alternative with an exogenous allele matching the reference allele (165). The exogenous allele can include an allele mapped to a base nucleotide of the reference genome 145 at the reference position 142 that was the result of a mapping error. The exogenous allele may be erroneously mapped to a first region of the reference genome that has a substantially similar or identical sequence of nucleotide bases to one or more second regions of the reference genome 145.

[0088] The modified set of probability results 150 includes respective non-conventional probabilities that account for the potential occurrence of one or more sequencing errors. These additional probabilities include (i) the likelihood that a read at reference position 142 indicates the occurrence of a homozygous reference with a sequencing error matching the alternate allele (166), and (ii) the likelihood that a read at reference position 142 indicates the occurrence of a homozygous alternate with a sequencing error matching the reference allele (167).

[0089] The set of modified probability results 150 collectively represent the probability that a true variant of the base nucleotide at the reference position 142 exists. Additionally, each probability of the set of modified probability results 150, including probabilities 161, 162, 163, 164, 165, 166, 167, provides a particular probability score that a particular genotype represented by each hypothesis exists, where the particular genotype can include a homozygous reference, a heterozygous alternative, a homozygous alternative, a homozygous reference with an exogenous allele matching the alternative, a homozygous alternative with an exogenous allele matching the reference allele, a homozygous reference with a sequencing error matching the alternative allele, or a homozygous alternative with a sequencing error matching the reference allele.

[0090] The set of modified probability results 150 can be used by the variant calling unit 130 to determine whether a candidate alternative allele of one or more sample reads is a true variant of a reference allele at a reference location of interest. Referring to the example of FIG. 1, the variant calling unit 130 can use the set of modified probability results 150 to determine whether a candidate alternative allele "G" shown in the interface 140 is a true variant of a reference allele "A" at a reference location. For example, the variant calling unit 130 processes the set of modified probability results 150 to determine at 136 whether data identifying a true variant at a reference location should be included in a variant call format (VCF) file 170 generated by the variant calling unit 130.

[0091] The variant calling unit 130 can determine whether a candidate alternative allele of one or more reads of a sample, such as the data describing a read having a candidate alternative allele "G" shown on the interface 140, should represent a true variant by determining a population score based on the modified probability result 150 and evaluating the population score using one or more pre-defined thresholds. The population score can indicate the likelihood that a true variant is present. In some implementations, the population score can be determined using equation (14). However, other methods of determining the population score are within the scope of this disclosure. If the computer determines at 136 that the population score meets the pre-defined threshold, then the computer can add information to the VCF file indicating that a true variant is present at the first position. The information added to the VCF file can include, for example, data indicating the position of the candidate allele, an identifier of the alternative allele, a genotype of the alternative allele, and a population score. Alternatively, if the computer determines at 136 that the population score does not meet the pre-defined threshold, then the computer can discard the output data and determine that the output data indicates the occurrence of a false positive at the reference position.

[0092] In the example shown in FIG. 1 , the variant calling unit 130 may determine 136 that the population score based on the probabilities of the set of altered probabilities 150 does not meet a predefined threshold and does not include information identifying the candidate alternative allele "G" in the VCF file 170 as a true variant. This is because the population score indicates a high likelihood that one or more sequencing errors are present because the candidate alternative allele "G" occurs only on a single strand at a low base quality forward aligned position away from the 5' end of the respective read, based on the set of altered probabilities 140. Thus, the variant calling unit 130 may determine, based on an evaluation of the set of altered probabilities 150, that the candidate alternative allele "G" is not a true variant at the reference location 142, but is instead a false positive, and no true variant is present at that reference location 142.

[0093] The probabilities shown in interface 160 are merely examples and are provided for purposes of illustrating one example of the present disclosure. The probabilities shown in 160 are not the result of the actual information described in FIG. 1 being put into the actual probability model described by this specification.

[0094] Figure 2 is a flow chart of an example process 200 for correlated error event mitigation for variant calling. The process 200 is described below as being implemented by a computer, such as the variant calling unit of Figure 1. In some implementations, the variant calling unit may be implemented using hardware, such as one or more FPGAs, ASICs, CPUs, GPUs, or combinations thereof, using software, or a combination thereof.

[0095] A computer can begin performing process 200 by accessing (step 210) a pileup of aligned array reads stored in one or more memory devices. The aligned array reads may have been generated using a mapping and alignment unit implemented using configurable logic gates of an FGPA device. In some implementations, the accessed reads may be stored in one or more memory devices and include a first set of aligned array reads in a forward direction and a second set of aligned array reads in a reverse direction. Each set of array reads corresponds to a particular read orientation or direction.

[0096] The computer may continue performing process 200 by obtaining information describing one or more characteristics of each reading in the pileup at the first location of the pileup of array reads (step 220). The one or more characteristics may include attributes of the readings in the pileup at the first location that may be used to account for the probability of occurrence of one or more correlated error events.

[0097] The computer may continue performing process 200 by providing (step 230) one or more inputs to a probabilistic model describing one or more characteristics of the pile-up reading at the first location. The one or more characteristics associated with the pile-up reading at the first location may include (i) information describing the one or more characteristics obtained from the one or more memory devices, (ii) information generated by one or more models, such as a P-HMM model, based on processing of the one or more models of information describing the one or more characteristics obtained from the one or more memory devices, or a combination thereof. In some implementations, the probabilistic model is configured to determine, for each hypothesis of the one or more hypotheses selected based on the one or more inputs, a score for each of the hypothesis indicating that the hypothesis is true.

[0098] The computer may continue performing process 200 by obtaining output information for each of the one or more hypotheses based on the one or more inputs. The output information for each hypothesis may (i) be generated by the probabilistic model based on processing of one or more inputs to the probabilistic model that describe one or more characteristics of a respective reading of the pile-up, and (ii) may indicate a probability that the hypothesis is true (step 240). In some implementations, the computer may obtain such output information for each of the one or more hypotheses, or a subset of the one or more hypotheses. Whether a particular hypothesis is included in the output information may be determined based on the inputs provided to the probabilistic model.

[0099] In some implementations, the one or more hypotheses may include, based on one or more inputs to a probabilistic model as described above, (i) the likelihood that a read at the reference position indicates the occurrence of a homozygous reference, (ii) the likelihood that a read at the reference position indicates the occurrence of a heterozygous substitution, (iii) the likelihood that a read at the reference position indicates the occurrence of a homozygous substitution, (iv) the likelihood that a read at the reference position indicates the occurrence of a homozygous reference with an exogenous allele matching the substitution, (v) the likelihood that a read at the reference position indicates the occurrence of a homozygous substitution with an exogenous allele matching the reference allele, (vi) the likelihood that a read at the reference position indicates the occurrence of a homozygous reference with a sequencing error matching the substitution allele, and (vii) the likelihood that a read at the reference position indicates the occurrence of a homozygous substitution with a sequencing error matching the reference allele, or any combination thereof.

[0100] The computer may continue performing process 200 by determining 250 the likelihood that a true variant is present at the first position based on the obtained output data generated by the probabilistic model for each of the multiple hypotheses. In some implementations, this determination may be made based on an individual or collective evaluation of the determined probability scores for each of the one or more hypotheses against one or more predefined thresholds, for example.

[0101] For example, the computer can determine a population score based on the output data generated by the probabilistic model for each of the multiple hypotheses. The population score can indicate the likelihood that a true variant exists. In some implementations, the population score can be determined using equation (14). However, other methods of determining the population score are within the scope of the present disclosure. If the computer determines that the population score meets a predetermined threshold, then the computer can add information to the VCF file indicating that a true variant exists at the first position. Alternatively, if the computer determines that the population score does not meet the predetermined threshold, then the computer can discard the output data and determine that the output data indicates the occurrence of a false positive at the reference position.

[0102] Figure 3 is a flow chart of an example process 300 of mapping error mitigation for variant calling. The process 300 is described below as being implemented by a computer, such as the variant calling unit of Figure 1. In some implementations, the variant calling unit may be implemented using hardware, such as one or more FPGAs, ASICs, CPUs, GPUs, or combinations thereof, using software, or a combination thereof.

[0103] A computer can begin performing process 300 by accessing (step 310) a pileup of aligned array reads stored in one or more memory devices. The aligned array reads may have been generated using a mapping and alignment unit implemented using configurable logic gates of an FGPA device. In some implementations, the accessed reads may be stored in one or more memory devices and include a first set of aligned array reads in a forward direction and a second set of aligned array reads in a reverse direction. Each set of array reads corresponds to a particular read orientation or direction.

[0104] The computer can continue performing process 300 by obtaining (320) information describing (i) a mapping confidence score for each of the reads in the pile-up of aligned sequence reads stored in one or more memories, and (ii) a read allele score for each of the reads in the pile-up of aligned sequence reads stored in one or more memories for each candidate allele at the reference position.

[0105] In some implementations, the mapping confidence score may include the output of the mapping and alignment unit 126, and may be calculated using the Phred scale of the mapping error, Q phred-mapping can be determined using Q phred-mapping =-10*log10(P e-mapping ) In this example, P e-mapping is the probability of mapping error for a particular read. The value of the mapping confidence score 144 can be proportional to the difference between the best alignment score from an alignment algorithm, such as the Smith-Waterman aligner, and the second best score of the aligner.

[0106] The read allele score can be determined in a number of different ways. For example, a more complex haplotype variant caller can use a read allele score calculated using Equation (1) above. Alternatively, other variant callers other than haplotype variant callers can use a read allele score calculated using Equation (2) above. In some implementations, the type of read allele score used can be determined based on whether the variant calling unit 130 is used to detect only SNPs or to detect SNPs and indels. For example, some implementations that use a variant calling unit to detect SNPs and indels can then use a read allele score calculated using Equation (1) above. As another example, other implementations that use a variant calling unit to detect only SNPs can then use a read allele score calculated using Equation (2).

[0107] The computer can continue performing process 300 by providing one or more inputs to a probability model describing the obtained information describing (i) a mapping confidence score for each of the reads in the pile-up of aligned sequence reads stored in one or more memories, and (ii) a read allele score for each of the reads in the pile-up of aligned sequence reads stored in one or more memories for each candidate allele at a reference position (step 330).

[0108] The computer can continue performing process 300 by obtaining (340) output information for each of the one or more hypotheses based on the one or more inputs obtained in 320. In the example of process 300, the obtained inputs include inputs of a mapping error probability model including (i) a mapping confidence score for each of the reads and (ii) a read allele score for each of the reads for each candidate allele at the reference position. Thus, based on receiving these inputs of the mapping error probability model, the computer generates output information for one or more hypotheses that consider the occurrence of one or more mapping errors. Such hypotheses include (i) the likelihood that the read at the reference position indicates the occurrence of a homozygous reference with an exogenous allele matching the substitution and (ii) the likelihood that the read at the reference position indicates the occurrence of a homozygous substitution with an exogenous allele matching the reference allele.

[0109] The output information includes information generated by the mapping error probability model based on processing of the probabilistic model of inputs to the probabilistic model describing (i) a mapping confidence score for each of the reads of the pile-up of aligned sequence reads stored in the one or more memories and (ii) a read allele score for each of the reads of the pile-up of aligned sequence reads stored in the one or more memories for each candidate allele at the reference position for each of these hypotheses. Further, the obtained output information includes, for each particular hypothesis that takes into account the occurrence of one or more mapping errors, a score indicating the likelihood that the hypothesis is true.

[0110] The computer may continue performing process 300 by determining (step 350) the likelihood that a true variant is present at the first position based on the obtained output data generated by the probabilistic model for each of the multiple hypotheses. In some implementations, this determination may be made based on an individual or collective evaluation of the determined probability scores for each of the one or more hypotheses against one or more predefined thresholds, for example.

[0111] For example, the computer can determine a population score based on the output data generated by the probabilistic model for each of the multiple hypotheses. The population score can indicate the likelihood that a true variant exists. In some implementations, the population score can be determined using equation (14). However, other methods of determining the population score are within the scope of the present disclosure. If the computer determines that the population score meets a predetermined threshold, then the computer can add information to the VCF file indicating that a true variant exists at the first position. Alternatively, if the computer determines that the population score does not meet the predetermined threshold, then the computer can discard the output data and determine that the output data indicates the occurrence of a false positive at the reference position.

[0112] The process 300 above describes a method of using a probabilistic model that can be used to account for the likelihood of mapping errors. Examples of probabilistic models that can be used to account for the likelihood of one or more mapping errors are described in more detail below.

[0113] In one implementation, the probability model is modified to incorporate the possibility that a pile-up of array reads stored in a memory device may contain multiple incorrectly mapped reads resulting in the occurrence of one or more mapping errors. The probability model can be applied to the following scenarios: (1) For each read r i The mapping quality μ indicates the confidence on the Fred scale that the mapping was done correctly. i Therefore, R={r i ,μ i :i=1...N R}, (2) the secondary sequences (i.e., other loci in the genome to which the reads align well) may be unknown and / or too numerous to tabulate, and (3) as input, i and candidate allele G m,φ The base quality P(ri |G m,φ ) is given.

[0114] In one implementation, the present disclosure modifies traditional variant calling probability models that complement the list of candidate genotypes with extended candidate genotypes, representing the hypothesis that the pile-up of sequence reads stored in the memory device contains a mixture of local and exogenous alleles.

[0115] In the case of a diploid genome with one foreign allele, the disclosed methodology involves the addition of the foreign allele F m The expanded candidate genotype G' m =[G m,1 G m,2 F m ]. Define the local allele G m,1 and G m,2 has an allele frequency of (1-β) / 2, while the foreign allele F m are each assumed to have an allele frequency β, where β is unknown.

[0116] Each expanded candidate genotype G m In contrast, the model is

[0117]

number

[0118] where U(t) is the Heavyside Unit function.

[0119]

number

[0120] P0(F m ) is the genotype occurring at the locus of interest [ρ F m ] and ρ is the reference allele at the locus of interest.

[0121]

number

[0122] is the joint probability P(G' m ,R) is an estimate.

[0123] A probabilistic model for determining the likelihood of occurrence of one or more mapping errors is the mapping quality μ of the incorrectly mapped reads. i and the prior probability of a variant (or cluster of variants) that mismaps read i. A probability model that determines the likelihood of one or more mapping errors occurring is based on the relationship between p F is that a sufficient number of variants occur at different locations to map the reads with a mapping quality of μ = -10log 10 (p F / P0(F m )) represents the prior probability of mismapping. The term 1 above is the mapping quality indicator μ i is nonzero if and only if p does not exceed this threshold. F It increases as p decreases. In general, the number of variants that have arisen at distant sites is unknown, and F Sweep the

[0124]

number

[0125] We find the value that maximizes , and thus test all hypotheses about the number of distant variants that an exogenous read may end up with.

[0126] In some implementations, the complexity of the probabilistic model for determining the likelihood of one or more mapping errors occurring is optimized. β and p F are swept independently and the results are evaluated over a continuous range of values,

[0127]

number

[0128] Note that evaluating p may appear to have high computational complexity. Depending on the resolution, the computational complexity may be very high. However, F is μ i needs to be swept only over a discrete set of values ​​corresponding to the values ​​of

[0129]

number

[0130] Typically, this is a small number. In particular, there is no need to sweep the β term at all; instead, the optimal value can be estimated as the fraction of readings where term 1 exceeds term 2, and this estimate usually has little effect on the results. Finally, the computational cost of this operation is usually insignificant compared to other parts of the system (e.g., computing a Hidden Markov Model (HMM)).

[0131] In some implementations, the mapping confidence score may need to be adjusted. In some implementations, the mapping confidence score reported by the mapping and alignment unit may represent an estimate of the Phred-scale confidence that the reading was correctly mapped, but in reality, this estimate may be inaccurate. Therefore, depending on the mapper, it may be beneficial to adjust the mapping confidence score to better match the true likelihood of incorrect mapping. In some implementations, the output value of the mapping and alignment unit, which represents a first mapping confidence score, such as a MAPQ score, is adjusted to the mapping confidence score μ using a function 400 shown in FIG. i can be converted to

[0132] In some implementations, the term β may be limited to the range [0, 0.5]. A value β=0.5 corresponds to a scenario where all readings for an alternative reference location are mapped to the reference location of interest. Although there may be multiple sources of extraneous readings, implying higher possible values ​​of β, limiting β to 0.5 may improve overall accuracy by recovering some true positives that would otherwise be suppressed.

[0133] In some implementations, the number of candidates is G m,1 =G m,2 =G m In this case, equation (3) can be simplified to:

[0134]

number

[0135] In some implementations, the above formula may assume an input that includes only reads that overlap with the genotyping event, i.e., reads that favor one allele over another. However, in some implementations, there may be ambiguity about which reads overlap with the event. In such scenarios, the following formula avoids the complications associated with including non-overlapping reads in the calculation:

[0136]

number

[0137] In some implementations, another variation of the probabilistic model for determining the likelihood that one or more mapping errors have occurred is θ ... i Let denote the strand orientation of read i, then values ​​0 and 1 indicate forward and reverse strand orientations, respectively. Using a modified probabilistic mapping model with strand awareness to determine the likelihood that there is one or more occurrences of mapping errors, three hypotheses can be evaluated. These three hypotheses include that the extraneous read occurs only on the forward strand, only on the reverse strand, or on both strands. Solutions include:

[0138]

number

[0139] Starting from equation (5), this variant becomes:

[0140]

number

[0141] Where Λ0={i:θ i =0}, Λ1={i:θ i =1}, Λ2={i:θ i =0 or 1}.

[0142] 5 is a flow chart of an example of a process for variant calling sequencing error mitigation. Process 500 is described below as being implemented by a computer, such as the variant calling unit of FIG. 1. In some implementations, the variant calling unit can be implemented using hardware, such as one or more FPGAs, ASICs, CPUs, GPUs, or combinations thereof, using software, or a combination thereof.

[0143] A computer may begin performing process 500 by accessing a pile-up of aligned array reads stored in one or more memory devices (step 510). The aligned array reads may have been generated using a mapping and alignment unit implemented using configurable logic gates of an FGPA device. In some implementations, the accessed reads may be stored in one or more memory devices and include a first set of aligned array reads in a forward direction and a second set of aligned array reads in a reverse direction. Each set of array reads corresponds to a particular read orientation or direction.

[0144] The computer can continue performing process 500 by obtaining (520) information describing: (i) a read orientation for each of the reads in the pile-up of aligned sequence reads stored in one or more memories; (ii) a position of each base at a reference location that references the 5' end of the read for each of the reads in the pile-up of aligned sequence reads stored in one or more memories; (iii) a read allele score for each of the reads in the pile-up of aligned sequence reads stored in one or more memories for each candidate allele at the reference location; and (iv) a base quality score for each read of the base at a reference location, such as position "0."

[0145] The read allele score can be determined in a number of different ways. For example, a more complex haplotype variant caller can use a read allele score calculated using Equation (1) above. Alternatively, other variant callers other than haplotype variant callers can use a read allele score calculated using Equation (2) above. In some implementations, the type of read allele score used can be determined based on whether the variant calling unit 130 is used to detect only SNPs or to detect SNPs and indels. For example, some implementations that use a variant calling unit to detect SNPs and indels can then use a read allele score calculated using Equation (1) above. As another example, other implementations that use a variant calling unit to detect only SNPs can then use a read allele score calculated using Equation (2).

[0146] The computer can continue performing process 500 by providing one or more inputs to a probabilistic model (step 530) describing the obtained information describing: (i) a read orientation for each of the reads in the pile-up of aligned sequence reads stored in one or more memories; (ii) a position of each base at a reference location referencing the 5' end of the read for each of the reads in the pile-up of aligned sequence reads stored in one or more memories; (iii) a read allele score for each of the reads in the pile-up of aligned sequence reads stored in one or more memories for each candidate allele at the reference location; and (iv) a base quality score for each read of the base at a reference location, such as position "0".

[0147] The computer can continue performing process 500 by obtaining (540) output information for each of the one or more hypotheses based on the one or more inputs obtained in 520. In an example of process 500, the obtained inputs include inputs for a sequencing error probability model including: (i) a read orientation for each of the reads of the pile-up of aligned sequence reads stored in one or more memories, (ii) a position of each base at a reference location that references the 5' end of the read for each of the reads of the pile-up of aligned sequence reads stored in one or more memories, (iii) a read allele score for each of the reads of the pile-up of aligned sequence reads stored in one or more memories for each candidate allele at the reference location, and (iv) a base quality score for each read of the base at the reference location, such as position "0". Thus, based on receiving these inputs of the sequencing error probability model, the computer generates output information for one or more hypotheses that consider the occurrence of one or more sequencing errors. Such hypotheses include (i) the possibility that the read at the reference position indicates the occurrence of a homozygous reference with a sequencing error matching the alternative allele and (ii) the possibility that the read at the reference position indicates the occurrence of a homozygous alternative with a sequencing error matching the reference allele.

[0148] The obtained output information includes information generated by the sequencing error probability model based on processing the probability model into a probability model that describes (i) for each of these hypotheses: (i) a read orientation for each of the reads of the pile-up of aligned sequence reads stored in the one or more memories, (ii) a position of each base at a reference location that references the 5' end of the read for each of the reads of the pile-up of aligned sequence reads stored in the one or more memories, (iii) a read allele score for each of the reads of the pile-up of aligned sequence reads stored in the one or more memories for each candidate allele at the reference location, and (iv) a base quality score for each read of the base at a reference location, such as position "0". Additionally, the obtained output information includes, for each particular hypothesis that takes into account the occurrence of one or more sequencing errors, a score indicating the likelihood that the hypothesis is true.

[0149] The computer may continue performing process 500 by determining the likelihood that a true variant is present at the first position based on the obtained output data generated by the probabilistic model for each of the multiple hypotheses (step 550). In some implementations, this determination may be made based on an individual or collective evaluation of the determined probability scores for each of the one or more hypotheses against one or more predefined thresholds, for example.

[0150] For example, the computer can determine a population score based on the output data generated by the probabilistic model for each of the multiple hypotheses. The population score can indicate the likelihood that a true variant exists. In some implementations, the population score can be determined using equation (14). However, other methods of determining the population score are within the scope of the present disclosure. If the computer determines that the population score meets a predetermined threshold, then the computer can add information to the VCF file indicating that a true variant exists at the first position. Alternatively, if the computer determines that the population score does not meet the predetermined threshold, then the computer can discard the output data and determine that the output data indicates the occurrence of a false positive at the reference position.

[0151] The process 500 above describes a method of using a probability model that can be used to account for the likelihood of sequencing errors. Examples of probability models that can be used to account for the likelihood of one or more sequencing errors, also referred to as systematic errors, are described in more detail below.

[0152] In one implementation, the probability model considers a probability P(r i |G m,φ ) is modified to take into account the observation that the base quality is not represented by the base quality used to calculate

[0153] A probabilistic model to determine the likelihood of occurrence of a sequencing error can be applied to the following scenarios:

[0154] (1) Errors appear to occur independently in each chain direction, since it is common for one chain direction to be damaged while the other chain is error-free.

[0155] (2) Errors are more likely to occur further from the 5' end of a read, and error rates often drop off sharply at a certain distance from the 5' end. Thus, when reads for a given strand are listed in order of decreasing distance from the 5' end, all errors are contained within the first subset of the list.

[0156] (3) Errors are often accompanied by a decrease in average base quality across the subset of erroneous reads, but not all erroneous reads have poor base quality, and in many cases the average base quality is not low enough to reflect the true error rate associated with these error events.

[0157] (4) Errors are often preceded by a homopolymer that matches the error; for example, a T==>G error is often preceded by a sequence of G.

[0158] Thus, the present disclosure provides a probabilistic model for determining the likelihood of occurrence of sequencing errors that takes into account the above four properties.

[0159] In one implementation, a probabilistic model for determining the likelihood of occurrence of a sequencing error that takes into account the above four properties can be realized by starting with the following term definitions:

[0160] Let θ denote the chain direction θ=0,1.

[0161] If the reading r θ,i If we order i=i by strand orientation and distance from the locus of interest to the 5' end, then i=1 is the furthest from the 5' end and is most likely to be affected by an error event.

[0162] Suppose q θ,i , the reading r θ,i Let denote the Fred scale base quality of the bases aligned with the locus of interest.

[0163] Suppose

[0164]

number

[0165] Let i = 1...n be the ordered readings for the strand direction θ. θ Let the average base quality of a subset of

[0166]

number

[0167] It is.

[0168] If we assume that the extended candidate genotype G' m =[G m,1 G m,2 E m,0 E m,1 ] is defined, then E m,θ is the error allele in strand direction θ.

[0169] Suppose L E,θ The base E just before the error in the chain direction θ m,θ is the length of the homopolymer that matches

[0170] Expanded candidate genotype G' m For each, calculate the following:

[0171]

number

[0172] Where:

[0173]

number

[0174] denotes the prior probability of an error event affecting a subset of reads as a function of the average base quality of the subset and the length of the homopolymer matching the error.

[0175] quantity

[0176]

number

[0177] is the joint probability P(G' m represents an estimate of the probability of error (R, R) of a sequence of reads that are ordered by decreasing distance from the 5' end, starting from the first, and that affect a subset of consecutive reads, but do not affect reads outside that subset, (2) error events occur independently on each strand, and (3) the prior probability of such an error event depends on the average base quality across the subset of reads affected by the error event and the strand orientation E m,θ is a function of the length of the homopolymer that corresponds to the base θ immediately preceding the error.

[0178] In some implementations, the number of candidates is reduced to G m,1 =G m,2 =G m We can only test the evaluation case where . In this case, the expression simplifies to:

[0179]

number

[0180] In some implementations, α θ The value of α θ = 0.5. Therefore, we can rewrite (7) as follows:

[0181]

number

[0182] In some implementations, the prior probability function of (6)

[0183]

number

[0184] is expressed as a general function that can have a wide range of shapes. In theory, this function can be trained on real data. However, some implementations of current probabilistic models for determining sequencing error have limited use in training this function on real data.

[0185]

number

[0186] can be used.

[0187] In some implementations, the above equations can be adapted to accommodate the likelihood of different error alleles for each strand direction. However, in some implementations, the probability model for determining the likelihood of one or more sequencing errors is E m,0 =E m,1 In such an implementation, the subscript θ is removed and then the error allele is denoted by E m It can be expressed as:

[0188] Figure 6 is another flow chart of an example process 600 of correlated error mitigation for variant calling. The process 600 is described below as being implemented by a computer, such as the variant calling unit of Figure 1. In some implementations, the variant calling unit can be implemented using hardware, such as one or more FPGAs, ASICs, CPUs, GPUs, or combinations thereof, using software, or a combination thereof.

[0189] A computer can begin performing process 600 by accessing (step 610) a pileup of aligned array reads stored in one or more memory devices. The aligned array reads may have been generated using a mapping and alignment unit implemented using configurable logic gates of an FGPA device. In some implementations, the accessed reads may be stored in one or more memory devices and include a first set of aligned array reads in a forward direction and a second set of aligned array reads in a reverse direction. Each set of array reads corresponds to a particular read orientation or direction.

[0190] The computer can continue performing process 600 by obtaining (620) information describing: (i) a read orientation for each of the reads in the pile-up of aligned sequence reads stored in one or more memories; (ii) a position of each base at a reference location referencing the 5' end of the read for each of the reads in the pile-up of aligned sequence reads stored in one or more memories; (iii) a mapping confidence score for each of the reads; (iv) a read allele score for each of the reads in the pile-up of aligned sequence reads stored in one or more memories for each candidate allele at the reference location; and (v) a base quality score for each read of the base at a reference location, such as position "0."

[0191] In some implementations, the mapping confidence score may include the output of the mapping and alignment unit 126, and may be calculated using the Phred scale of the mapping error, Q phred-mapping can be determined using Q phred-mapping =-10*log10(P e-mapping ) In this example, P e-mapping is the probability of mapping error for a particular read. The value of the mapping confidence score 144 can be proportional to the difference between the best alignment score from an alignment algorithm, such as the Smith-Waterman aligner, and the second best score of the aligner.

[0192] The read allele score can be determined in a number of different ways. For example, a more complex haplotype variant caller can use a read allele score calculated using Equation (1) above. Alternatively, other variant callers other than haplotype variant callers can use a read allele score calculated using Equation (2) above. In some implementations, the type of read allele score used can be determined based on whether the variant calling unit 130 is used to detect only SNPs or to detect SNPs and indels. For example, some implementations that use a variant calling unit to detect SNPs and indels can then use a read allele score calculated using Equation (1) above. As another example, other implementations that use a variant calling unit to detect only SNPs can then use a read allele score calculated using Equation (2).

[0193] The computer can continue performing process 600 by providing (630) one or more inputs to a probabilistic model describing the obtained information describing: (i) a read orientation for each of the reads in the pile-up of aligned sequence reads stored in one or more memories; (ii) a position of each base at a reference location referencing the 5' end of the read for each of the reads in the pile-up of aligned sequence reads stored in one or more memories; (iii) a mapping confidence score for each of the reads; (iv) a read allele score for each of the reads in the pile-up of aligned sequence reads stored in one or more memories (step 630) for each candidate allele at the reference location; and (v) a base quality score for each read of the base at a reference location, such as position "0".

[0194] The computer can continue performing process 600 by obtaining (640) output information for each of the one or more hypotheses based on the one or more inputs obtained in 620. In an example of process 600, the obtained inputs include inputs for the mapping error probability model and the sequencing error probability model including: (i) a read orientation for each of the reads of the pileup of aligned sequence reads stored in one or more memories, (ii) a position of each base at a reference location that references the 5' end of the read for each of the reads of the pileup of aligned sequence reads stored in one or more memories, (iii) a mapping confidence score for each of the reads, (iv) a read allele score for each of the reads of the pileup of aligned sequence reads stored in one or more memories for each candidate allele at the reference location, and (v) a base quality score for each read of the base at the reference location, such as position "0". Thus, based on receiving these inputs of the mapping error probability model and the sequencing error probability model, the computer generates output information for one or more hypotheses that consider the occurrence of one or more mapping errors and one or more sequencing errors. Such hypotheses include: (i) the possibility that the read at the reference location indicates the occurrence of a homozygous reference with an exogenous allele matching the substitution; (ii) the possibility that the read at the reference location indicates the occurrence of a homozygous substitution with an exogenous allele matching the reference allele; (iii) the possibility that the read at the reference location indicates the occurrence of a homozygous reference with a sequencing error matching the substitution allele; and (iv) the possibility that the read at the reference location indicates the occurrence of a homozygous substitution with a sequencing error matching the reference allele.

[0195] The obtained output information includes information generated by the mapping error probability model and the sequencing error probability model based on processing of each probability model into a probability model that describes, for each of these hypotheses, (i) a read orientation for each of the reads of the pileup of aligned sequence reads stored in one or more memories, (ii) a position of each base at a reference location that references the 5' end of the read for each of the reads of the pileup of aligned sequence reads stored in one or more memories, (iii) a mapping confidence score for each of the reads, (iv) a read allele score for each of the reads of the pileup of aligned sequence reads stored in one or more memories for each candidate allele at the reference location, and (v) a base quality score for each read of the base at a reference location, such as position "0" (620). Additionally, the obtained output information includes, for each particular hypothesis that considers the occurrence of one or more mapping errors or one or more sequencing errors, a score indicating the likelihood that the hypothesis is true.

[0196] The computer may continue performing process 600 by determining (step 650) the likelihood that a true variant is present at the first position based on the obtained output data generated by the probabilistic model for each of the multiple hypotheses. In some implementations, this determination may be made based on an individual or collective evaluation of the determined probability scores for each of the one or more hypotheses against one or more predefined thresholds, for example.

[0197] For example, the computer can determine a population score based on the output data generated by the probabilistic model for each of the multiple hypotheses. The population score can indicate the likelihood that a true variant exists. In some implementations, the population score can be determined using equation (14). However, other methods of determining the population score are within the scope of the present disclosure. If the computer determines that the population score meets a predetermined threshold, then the computer can add information to the VCF file indicating that a true variant exists at the first position. Alternatively, if the computer determines that the population score does not meet the predetermined threshold, then the computer can discard the output data and determine that the output data indicates the occurrence of a false positive at the reference position.

[0198] The processes described with reference to the flowcharts of Figures 2, 3, and 5 generally describe a process for correlated error event mitigation of variant calling. These respective processes describe a comprehensive process for correlated error event mitigation with reference to Figure 2, a process for mapping error mitigation with reference to Figure 3, and a process for sequencing error mitigation with reference to Figure 5. However, when sufficient inputs are obtained by the computer to take into account one or more mapping errors and one or more sequencing errors, the present disclosure does not require that separate disjoint probability calculations for mapping errors and sequencing errors, respectively, be performed. Instead, in some implementations, a complete probability model can be used in a process such as process 600 to take into account mapping errors and sequencing errors. Although a combined probability model for taking into account mapping errors and sequencing errors is described below, it is not necessary to use a combined probability model. Instead, as described with reference to other implementations, separate disjoint probability models may be relied upon.

[0199] In the following discussion, we derive the complete probability calculation for the case of a diploid genome with at most one foreign allele and one systematic error allele. However, this can be directly extended to more alleles, taking into account the discussion below. The extended candidate genotype is G' m =[G m,1 G m,2 F m E m In the formula below, the notation 0 denotes the reference allele, 1 denotes the first alternative allele, 2 denotes the second alternative allele, etc. F m or E m A dash indicates the absence of an exogenous or systematic error allele.

[0200] Candidate G' m When does not contain any foreign or error alleles, the following formula is generated:

[0201]

number

[0202] G' m When G includes an exogenous allele, equation (4) or one of the variants above can be used. m When contains an error allele, equation (10) or one of its variants can be used. In general, it is not necessary to test both the cases of foreign and error alleles, since it is very rare to find pileups that are affected by both types of errors simultaneously.

[0203] As an example of a candidate list for the common case of a single ALT allele (REF=0, ALT=1), the following extended candidate genotypes can be tested:

[0204]

number

[0205] For each (unexpanded) candidate G m For the joint probability P(G m ,R) is G m is the maximum number of candidates that match

[0206]

number

[0207] and the posterior probability is simply

[0208]

number

[0209] where

[0210]

number

[0211] It is.

[0212] The examples described below with reference to Figures 9A to 12B show examples of the calculations described with reference to Figures 3, 5 and 6 for various pile-ups.

[0213] 7 is another flow chart of an example process 700 of correlated error mitigation for variant calling. The process 700 is described below as being implemented by a computer, such as the variant calling unit of FIG 1. In some implementations, the variant calling unit may be implemented using hardware, such as one or more FPGAs, ASICs, CPUs, GPUs, or combinations thereof, using software, or a combination thereof.

[0214] A computer storing one or more probability models can begin performing process 700 by receiving (710) input data including information describing one or more characteristics of a read from a pileup of aligned sequence reads. In some implementations, the one or more characteristics can include (i) information describing a read orientation for each of the reads in the pileup of aligned sequence reads, (ii) information describing a position of each base at a reference location that references the 5' end of the read for each of the reads in the pileup of aligned sequence reads, (iii) a mapping quality score for each of the reads in the pileup of aligned sequence reads, (iv) a read allele score for each of the reads in the pileup of aligned sequence reads for each candidate allele at the reference location, and (v) a base quality score for each read of the base at a reference location, such as position "0". In some implementations, all or a subset of these inputs can be provided as inputs to the probability model, as well as references to other probability models herein. In some implementations, the one or more probability models stored by the computer can include a mapping error probability model and a sequencing error probability model.

[0215] In some implementations, the mapping confidence score may include the output of the mapping and alignment unit 126, and may be calculated using the Phred scale of the mapping error, Q phred-mapping can be determined using Q phred-mapping =-10*log10(P e-mapping ) In this example, P e-mapping is the probability of mapping error for a particular read. The value of the mapping confidence score 144 can be proportional to the difference between the best alignment score from an alignment algorithm, such as the Smith-Waterman aligner, and the second best score of the aligner.

[0216] The read allele score can be determined in a number of different ways. For example, a more complex haplotype variant caller can use a read allele score calculated using Equation (1) above. Alternatively, other variant callers other than haplotype variant callers can use a read allele score calculated using Equation (2) above. In some implementations, the type of read allele score used can be determined based on whether the variant calling unit 130 is used to detect only SNPs or to detect SNPs and indels. For example, some implementations that use a variant calling unit to detect SNPs and indels can then use a read allele score calculated using Equation (1) above. As another example, other implementations that use a variant calling unit to detect only SNPs can then use a read allele score calculated using Equation (2).

[0217] The computer can continue performing process 700 by determining (720) a set of one or more hypotheses based on the received input. For example, if the one or more inputs include information describing characteristics of the reads associated with a probabilistic model that takes into account one or more mapping errors, then the computer can determine a set of one or more hypotheses that take into account one or more mapping errors. Such hypotheses can include, for example, (i) the likelihood that the read at the reference position indicates the occurrence of a homozygous reference with an exogenous allele that matches the substitution, and (ii) the likelihood that the read at the reference position indicates the occurrence of a homozygous substitution with an exogenous allele that matches the reference allele.

[0218] Alternatively, or in addition, for example, if the one or more inputs include information describing characteristics of a read associated with a probabilistic model that takes into account one or more sequencing errors, then the computer can determine a set of one or more hypotheses that take into account one or more sequencing errors. Such hypotheses can include (i) the likelihood that the read at the reference position indicates the occurrence of a homozygous reference with a sequencing error matching the alternative allele, and (ii) the likelihood that the read at the reference position indicates the occurrence of a homozygous alternative with a sequencing error matching the reference allele.

[0219] In some implementations, the one or more received inputs may cause the computer to determine a set of hypotheses that consider both mapping errors and sequencing errors. Such a determination may be performed by the computer, for example, when the one or more received inputs include information describing the characteristics of the reads associated with a probabilistic model that considers both one or more mapping errors and one or more sequencing errors. The set of hypotheses that consider both one or more mapping errors and one or more sequencing errors may include, for example, (i) the likelihood that the read at the reference position indicates the occurrence of a homozygous reference with an exogenous allele that matches the substitution, (ii) the likelihood that the read at the reference position indicates the occurrence of a homozygous substitution with an exogenous allele that matches the reference allele, (iii) the likelihood that the read at the reference position indicates the occurrence of a homozygous reference with a sequencing error that matches the substitution allele, and (iv) the likelihood that the read at the reference position indicates the occurrence of a homozygous substitution with a sequencing error that matches the reference allele.

[0220] The computer may continue performing process 700 by determining a respective hypothesis score for each of the one or more hypotheses at 730, which indicates the probability that the respective hypothesis determined at step 720 is true. Determining the score for each hypothesis by the computer may include, for example, using one or more probability models to determine a probability score for each of the one or more hypotheses. An example of a probability model that may be used to determine a score for each mapping error related hypothesis is described above with reference to equation (3). However, other variations of other probability models are also described above to calculate a score for each mapping error hypothesis. An example of a probability model that may be used to determine a score for each sequencing error related hypothesis is described above with reference to equation (9). However, other variations of other probability models are also described above to calculate a score for each sequencing error hypothesis.

[0221] The computer can continue performing process 700 by providing (740) output data generated by the probabilistic model, including a score for each of the one or more hypotheses determined in step 720. In some implementations, the provided output data can be provided to a second computer configured to determine the likelihood of a true variant being present based on the output data. In some implementations, the second computer can be the computer that provided the input in step 710. In some implementations, the second computer can be another genomic analysis model of a secondary analysis unit, or other computer module, of which the computer is a part. In some implementations, the second computer can be a remote computer with a secondary analysis unit, including a mapper and alignment module, but not including a variant calling module, that can communicate with the computer using one or more networks.

[0222] However, in other implementations, the computer need not provide input to a second computer. Instead, the output information provided at 740, including a set of scores for each of the one or more hypotheses, can be used by a computer, such as a variant caller, to determine whether a candidate allele at a reference position is a true variant or a false positive, as described with reference to FIG.

[0223] System Components

[0224] FIG. 8 is a block diagram of system components that can be used to implement a system for correlated error mitigation in variant calling.

[0225] Computing device 800 is intended to represent various forms of digital computers, such as laptops, desktops, workstations, personal digital assistants, servers, blade servers, mainframes, and other suitable computers. Computing device 850 is intended to represent various forms of mobile devices, such as personal digital assistants, mobile phones, smart phones, and other similar computing devices. In addition, computing device 800 or 850 may include a Universal Serial Bus (USB) flash drive. The USB flash drive may store an operating system and other applications. The USB flash drive may include input / output components, such as a wireless transmitter or a USB connector that may be inserted into a USB port of another computing device. The components shown herein, their connections and relationships, and their functions are for illustrative purposes only and are not intended to limit the implementation of the invention described and / or claimed in this document.

[0226] Computing device 800 includes a processor 802, a memory 804, a storage device 808, a high-speed interface 808 connecting to memory 804 and a high-speed expansion port 810, and a low-speed interface 812 connecting to a low-speed bus 814 and storage device 808. Each of components 802, 804, 808, 808, 810, and 812 are interconnected using various buses and may be mounted on a common motherboard or in other manners as needed. Processor 802 may process instructions for execution within computing device 800, including instructions stored in memory 804 or storage device 808, to display graphical information of a GUI on an external input / output device, such as a display device 816 coupled to high-speed interface 808. In other implementations, multiple processors and / or multiple buses may be used, along with multiple memories and multiple types of memories, as needed. Multiple computing devices 800 may also be connected with each device providing a portion of the operations required, for example, as a server bank, a group of blade servers, or a multiprocessor system.

[0227] The memory 804 stores information within the computing device 800. In one implementation, the memory 804 is a volatile memory unit(s). In another implementation, the memory 804 is a non-volatile memory unit(s). The memory 804 may also be another form of computer-readable medium, such as a magnetic or optical disk.

[0228] The storage device 808 can provide mass storage for the computing device 800. In one implementation, the storage device 808 can be or include a computer-readable medium, such as a floppy disk device, a hard disk device, an optical disk device, a tape device, a flash memory, or other similar solid-state memory device, or an array of devices, including a storage area network or other configuration of devices. The computer program product can be tangibly embodied in an information carrier. The computer program product can also include instructions that, when executed, perform one or more methods, such as those described above. The information carrier is a computer or machine-readable medium, such as the memory 804, the storage device 808, or a memory on the processor 802.

[0229] The high-speed controller 808 manages bandwidth-intensive operations while the low-speed controller 812 manages low-bandwidth-intensive operations. Such allocation of functions is merely exemplary. In one implementation, the high-speed controller 808 is coupled to the memory 804, the display device 816, e.g., via a graphics processor or accelerator, and to a high-speed expansion port 810 that can accept various expansion cards (not shown). In this implementation, the low-speed controller 812 is coupled to the storage device 808 and the low-speed expansion port 814. The low-speed expansion port, which can include various communication ports, e.g., USB, Bluetooth, Ethernet, wireless Ethernet, can be connected to one or more input / output devices, such as a keyboard, a pointing device, a microphone / speaker set, a scanner, or a network device, such as a switch or router, e.g., via a network adapter. The computing device 800 can be implemented in several different forms, as shown in the figure. For example, it can be implemented as a standard server 820 or multiple times in a group of such servers. It can also be implemented as part of a rack server system 824. It may further be implemented in a personal computer, such as a laptop computer 822. Alternatively, components from computing device 800 may be combined with other components in a mobile device (not shown), such as device 850. Each such device may include one or more of computing devices 800, 850, and the overall system may be made up of multiple computing devices 800, 850 in communication with each other.

[0230] The computing device 800 can be implemented in a number of different forms, as shown in the figure. For example, it can be implemented as a standard server 820, or multiple times in a group of such servers. It can also be implemented as part of a rack server system 824. It can also be implemented in a personal computer, such as a laptop computer 822. Alternatively, components from the computing device 800 can be combined with other components in a mobile device (not shown), such as device 850. Each such device can include one or more of the computing devices 800, 850, and the entire system can be made up of multiple computing devices 800, 850 in communication with each other.

[0231] Computing device 850 includes, among other components, a processor 852, memory 864, and input / output devices such as a display device 854, a communication interface 866, and a transceiver 868. Device 850 may also include a storage device, such as a microdrive or other device, to provide additional storage. Each of components 850, 852, 864, 854, 866, and 868 are interconnected using various buses, and some of the components may be mounted on a common motherboard or in other manners as desired.

[0232] The processor 852 can execute instructions within the computing device 850, including instructions stored in the memory 864. The processor can be implemented as a chipset of chips including separate analog and digital processors. Additionally, the processor can be implemented using any of a number of architectures. For example, the processor 810 can be a CISC (Complex Instruction Set Computer) processor, a RISC (Reduced Instruction Set Computer) processor, or a MISC (Minimum Instruction Set Computer) processor. The processor can provide coordination of other components of the device 850, such as, for example, control of a user interface, applications launched by the device 850, and wireless communication by the device 850.

[0233] The processor 852 can communicate with a user through a control interface 858 and a display interface 856 coupled to a display device 854. The display device 854 can be, for example, a TFT (thin film transistor liquid crystal) display or an OLED (organic light emitting diode) display, or other suitable display technology. The display interface 856 can have appropriate circuitry for driving the display device 854 to present graphical and other information to the user. The control interface 858 can receive commands from a user and translate them for submission to the processor 852. Additionally, an external interface 862 can be provided in communication with the processor 852 to enable short-range communication of the device 850 with other devices. The external interface 862 can be provided, for example, for wired communication in some implementations or for wireless communication in other implementations, and multiple interfaces can also be used.

[0234] The memory 864 stores information within the computing device 850. The memory 864 can be implemented as one or more of a computer-readable medium, a volatile memory unit(s), or a non-volatile memory unit(s). An expansion memory 874 can also be provided and connected to the device 850 through an expansion interface 872, which can include, for example, a SIMM (single in-line memory module) card interface. Such expansion memory 874 can provide additional storage space for the device 850 or can store applications or other information for the device 850. In particular, the expansion memory 874 can include instructions that perform or supplement the above processes and can also include secure information. Thus, for example, the expansion memory 874 can be provided as a security module for the device 850 and can be programmed with instructions that enable secure use of the device 850. Furthermore, secure applications can be provided via a SIMM card along with additional information, such as disposing identification information on the SIMM card in an unhackable manner.

[0235] The memory may include, for example, flash memory and / or NVRAM memory, as described below. In one implementation, the computer program product is tangibly embodied in an information carrier. The computer program product includes instructions that, when executed, perform one or more methods, such as those described above. The information carrier is a computer or machine-readable medium, such as memory 864, expansion memory 874, or memory on processor 852, which may be received, for example, via transceiver 868 or external interface 862.

[0236] Device 850 can communicate wirelessly through communication interface 866, which can include digital signal processing circuitry as needed. Communication interface 866 can provide for communication under various modes or protocols, such as GSM voice calls, SMS, EMS, or MMS messaging, CDMA, TDMA, PDC, WCDMA, CDMA2000, or GPRS, among others. Such communication can occur, for example, through radio frequency transceiver 868. Additionally, short-range communication can occur, such as using Bluetooth, Wi-Fi, or other such transceivers (not shown). In addition, GPS (Global Positioning System) receiver module 870 can provide additional navigation and location-related wireless data to device 850, which can be used as needed by applications launched on device 850.

[0237] Device 850 can also communicate audibly using voice codec 860, which can receive voice information from a user and convert it into usable digital information. Similarly, voice codec 860 can generate sounds audible to a user, such as through a speaker in a handset of device 850. Such sounds can include sounds from a voice call, can include recorded sounds, such as voice messages, music files, and the like, and can also include sounds generated by applications running on device 850.

[0238] The computing device 850 may be implemented in a number of different forms, as shown in the figure, such as as a mobile phone 880, or as part of a smartphone 882, personal digital assistant, or other similar mobile device.

[0239] Various implementations of the systems and methods described herein can be realized in digital electronic circuitry, integrated circuits, specially designed ASICs (application-specific integrated circuits), computer hardware, firmware, software, and / or combinations of such implementations. These various implementations can include implementations in one or more computer programs executable and / or interpretable on a programmable system including at least one programmable processor, which can be special purpose or general purpose, coupled to receive and transmit data and instructions between a storage system, at least one input device, and at least one output device.

[0240] These computer programs (also known as programs, software, software applications, or code) contain machine instructions for a programmable processor and may be implemented in high level procedural and / or object-oriented programming languages, and / or assembly / machine languages. As used herein, the terms "machine-readable medium" and "computer-readable medium" refer to any computer program product, apparatus, and / or device used to provide machine instructions and / or data to a programmable processor, including a machine-readable medium that receives the machine instructions as a machine-readable signal, e.g., magnetic disks, optical disks, memory, programmable logic devices (PLDs). The term "machine-readable signal" refers to any signal used to provide machine instructions and / or data to a programmable processor.

[0241] To provide for user interaction, the systems and techniques described herein can be implemented in a computer having a display device, e.g., a CRT (cathode ray tube) or LCD (liquid crystal display) monitor, for displaying information to the user, and a keyboard and pointing device, e.g., a mouse or trackball, by which the user can provide input to the computer. Other types of devices can also be used to provide for user interaction, e.g., feedback provided to the user can be any form of sensory feedback, e.g., visual feedback, auditory feedback, or haptic feedback, and input from the user can be received in any form, including acoustic, speech, or tactile input.

[0242] The systems and techniques described herein can be implemented in a computing system that includes back-end components, e.g., data servers, or includes middleware components, e.g., application servers, or includes front-end components, e.g., client computers, having a graphical user interface or web browser through which a user can interact with an implementation of the systems and techniques described herein, or any combination of such back-end, middleware, or front-end components. The components of the system can be interconnected by any form or medium of digital data communication, e.g., a communications network. Examples of communications networks include a local area network ("LAN"), a wide area network ("WAN"), and the Internet.

[0243] A computing system may include clients and servers. Generally, clients and servers are remote from each other and typically interact through a communication network. The relationship of client and server arises by virtue of computer programs running on the respective computers and having a client-server relationship to each other.

[0244] Several embodiments have been described. Nevertheless, it should be understood that various modifications may be made without departing from the spirit and scope of the present invention. Additionally, the logic flows depicted in the figures do not require the particular order shown, or sequence, to achieve desirable results. Additionally, other steps may be provided or steps may be omitted from the described flows, and other components may be added or removed from the described systems. Accordingly, other embodiments are within the scope of the following claims.

[0245] [Example]

[0246] The subject matter of the present disclosure is further illustrated with reference to the following examples, which in no way limit the scope of the disclosure.

[0247] The examples provided in this section provide actual examples illustrating the computation of the probabilistic models described herein on actual pileups. Each pileup plot includes the MAPQ mapping confidence score for each read, the base quality at a reference position, such as position "0", and the average base quality per strand (blue=forward, red=backward).

[0248] Example 1 - Exemplary true positive variants

[0249] FIG. 9A is an example of a summary of experimental results from an implementation of a process for correlated error mitigation in variant calling performed on a pile-up of sequencing reads, with the results showing an example of a true positive result.

[0250] More specifically, Figure 9A displays an image 940 of a pileup 941 of aligned sequence reads showing an example of a typical true positive. A true positive is an example of a pileup 941 of reads that has a true variant at a reference location 942. In this example, it can be seen that the characteristics of the reads indicate that the pileup 941 of reads has a candidate alternative allele "C" at the reference location 942 that differs from the reference allele "T" in the reference genome to which it is mapped and aligned.

[0251] In this example, and with reference to image 940 of pileup 941 in FIG. 9A , it can be seen that the alternative allele frequencies for "C" are balanced between the first forward read direction (or orientation) and the second reverse read direction (or orientation), the MAPQ mapping confidence scores are high for each of the reads at the reference location, most of the mapping confidence scores 944 are up to "250", the average base quality for both strand directions (or orientations) is high, and most of the base quality scores 943 are "35" or higher. As a result, the probability score for candidate [10|-|-] 962 is high, as shown in the probability score results column 980 and the normalized probability score results column 980. Thus, information identifying or otherwise associated with the true alternative "C" can be included in the VCF file.

[0252] A complete set of modified probability results 960 is shown in FIG. 9B. Candidate [10|-|-] 962 corresponds to hypothesis 162 and includes the probability that a read at reference position 142 indicates the occurrence of a heterozygous alternative allele. Probability score results column 980 and normalized probability score results column 990 show the probability score results and normalized score results, respectively, for each of the other conventional hypotheses 961, 962, 963 and non-conventional hypotheses 964, 965, 966, 967. For clarity, conventional hypothesis 961 corresponds to conventional hypothesis 161 above, conventional hypothesis 962 corresponds to conventional hypothesis 162 above, and conventional hypothesis 963 corresponds to conventional hypothesis 163 above. Furthermore, non-conventional hypothesis 964 corresponds to non-conventional hypothesis 164 above, non-conventional hypothesis 965 corresponds to non-conventional hypothesis 165 above, non-conventional hypothesis 966 corresponds to non-conventional hypothesis 166 above, and non-conventional hypothesis 967 corresponds to non-conventional hypothesis 167 above.

[0253] Example 2 - Low probability of mapping errors

[0254] FIG. 10A is an example of a summary of experimental results from an implementation of a process for correlated error mitigation in variant calling performed on a pile-up of sequencing reads, where the results indicate a low probability of occurrence of mapping errors.

[0255] More specifically, FIG. 10A displays an image 1040 of a pileup 1041 of aligned sequence reads, showing, for example, possible foreign reads, or possible mapping errors. In this example, the pileup has a low alternative allele frequency, and a rather low MAPQ mapping confidence score 1044 for reads with alternative alleles. Both of these factors will be exploited by the mapping error probability model provided by this disclosure. Thus, the result is that it is unlikely to represent a true variant. A review of the probability scores for each hypothesis indicates that one cannot draw a conclusion with high confidence that the "C" allele is a foreign read, nor can one make a heterozygous call with high confidence. Thus, the output information may be discarded without adding any information identifying the candidate alternative alleles to the VCF file.

[0256] A complete set of modified probability results 1060 is shown in FIG. 10B. Probability score results column 1080 and normalized probability score results column 1090 show the probability score results and normalized score results, respectively, for each of the other conventional hypotheses 1061, 1062, 1063 and non-conventional hypotheses 1064, 1065, 1066, 1067. For clarity, conventional hypothesis 1061 corresponds to conventional hypothesis 161 above, conventional hypothesis 1062 corresponds to conventional hypothesis 162 above, and conventional hypothesis 1063 corresponds to conventional hypothesis 163 above. Additionally, non-conventional hypothesis 1064 corresponds to non-conventional hypothesis 164 above, non-conventional hypothesis 1065 corresponds to non-conventional hypothesis 165 above, non-conventional hypothesis 1066 corresponds to non-conventional hypothesis 166 above, and non-conventional hypothesis 1067 corresponds to non-conventional hypothesis 167 above.

[0257] Example 3 - Unlikely true variants due to sequencing errors

[0258] FIG. 11A is an example of a summary of experimental results from an implementation of a process for correlated error mitigation of variant calling performed on pile-ups of sequencing reads, where the results indicate a high probability that a candidate alternative allele is unlikely to be a true variant due to sequencing error.

[0259] More specifically, FIG. 11A displays an image 1140 of a pile-up 1141 of aligned sequence reads exhibiting systematic, or sequencing, errors.

[0260] In this example, all of the "G" alternative alleles occur in a single read orientation in the forward aligned direction. Furthermore, all of the "G" alternative alleles occur in the subset of forward oriented reads furthest from the 5' end of the reads. The subset of reads with the "G" alternative allele has very low base quality, as evident from the base quality score 1143. Furthermore, the "G" alternative allele matches two base alleles in the reference genome that present the base allele at the current reference position 1145. The sequencing error probability model takes into account the aforementioned characteristics of the reads, and a probability score is output for each of the seven hypotheses 1161, 1162, 1163, 1164, 1165, 1166, 1167, shown in the probability score result 1180 and normalized probability score result column 1190, supporting, with high confidence, that the "G" alternative allele is unlikely to support the true variant. Thus, the output information may be discarded without adding any information identifying the candidate alternative allele to the VCF file.

[0261] A complete set of modified probability results 1160 is shown in FIG. 11B. Probability score results column 1180 and normalized probability score results column 1190 show the probability score results and normalized score results, respectively, for each of the other conventional hypotheses 1161, 1162, 1163 and non-conventional hypotheses 1164, 1165, 1166, 1167. For clarity, conventional hypothesis 1161 corresponds to conventional hypothesis 161 above, conventional hypothesis 1162 corresponds to conventional hypothesis 162 above, and conventional hypothesis 1163 corresponds to conventional hypothesis 163 above. Additionally, non-conventional hypothesis 1164 corresponds to non-conventional hypothesis 164 above, non-conventional hypothesis 1165 corresponds to non-conventional hypothesis 165 above, non-conventional hypothesis 1166 corresponds to non-conventional hypothesis 166 above, and non-conventional hypothesis 1167 corresponds to non-conventional hypothesis 167 above.

[0262] Example 4 - Unlikely true variants due to sequencing errors in both reading orientations

[0263] FIG. 12A is an example of a summary of experimental results from an implementation of the process for correlated error mitigation of variant calling performed on a pile-up of sequencing reads, where the results indicate a high probability that the candidate alternative allele is unlikely to be a true variant due to sequencing errors in both read orientations.

[0264] More specifically, FIG. 12A displays an image 1240 of a pile-up 1241 of aligned sequence reads that exhibit systematic, or sequencing, errors in both read orientations.

[0265] In this example, the base quality is extremely degraded, as is evident from the base quality score 1243. For reads sequenced in the forward direction, the "G" alternative allele is restricted to a subset of reads where the locus of interest is far from the 5' end. Furthermore, the "G" alternative allele matches three preceding reference alleles in the reference genome that precede the reference allele at reference position 1242. The sequencing error probability model takes into account the aforementioned characteristics of the reads, and a probability score is output for each of the seven hypotheses 1261, 1262, 1263, 1264, 1265, 1266, 1267, shown in the probability score result 1280 and normalized probability score result column 1290, supporting with high confidence that the "G" alternative allele is unlikely to support the true variant. Thus, the output information may be discarded without adding any information identifying the candidate alternative allele to the VCF file.

[0266] A complete set of modified probability results 1260 is shown in Figure 12B. Probability score results column 1280 and normalized probability score results column 1290 show the probability score results and normalized score results, respectively, for each of the other conventional hypotheses 1161, 1262, 1263, and non-conventional hypotheses 1264, 1265, 1266, 1267. For clarity, conventional hypothesis 1261 corresponds to conventional hypothesis 161 above, conventional hypothesis 1262 corresponds to conventional hypothesis 162 above, and conventional hypothesis 1263 corresponds to conventional hypothesis 163 above. Furthermore, non-conventional hypothesis 1264 corresponds to non-conventional hypothesis 164 described above, non-conventional hypothesis 1265 corresponds to non-conventional hypothesis 165 described above, non-conventional hypothesis 1266 corresponds to non-conventional hypothesis 166 described above, and non-conventional hypothesis 1267 corresponds to non-conventional hypothesis 167 described above.

[0267] Other embodiments

[0268] While the present invention has been described in conjunction with the drawings and its detailed description, it is to be understood that the foregoing description is intended to be illustrative, and not limiting, of the scope of the invention as defined by the appended claims. Other aspects, advantages, and modifications are within the scope of the appended claims. [Explanation of symbols]

[0269] 100 A System for Reducing Correlated Errors in Variant Calling 105 Biological Samples 110 Nucleic acid sequencer 112 readings 120 Secondary Analysis Unit 122 Memory 124 Field Programmable Gate Array (FPGA) 126 Mapping and Alignment Units 128 Paired Hidden Markov Model (P-HMM) units 130 Valiant Calling Unit 131 Probability Models 132 Mapping Error Probability Model 134 Sequencing Error Probability Model 140 Interface 141 Pileup 142 Reference position 143 Base Quality Score 144 Mapping Confidence Score 145 Reference Genomes 150 Probability Results 160 Interface 170 files 200 processes 210 Steps 220 Steps 230 Steps 240 steps 250 steps 300 processes 310 Steps 330 Steps 350 steps 500 processes 510 Steps 530 steps 550 steps 600 processes 610 Steps 630 steps 650 steps 700 processes 800 computing devices 802 Processor 804 Memory 808 Components 810 Components 812 Low Speed ​​Controller 822 Laptop Computer 824 Rack Server System 850 Computing Devices 852 processor 854 Display device 856 Display Interface 858 Control Interface 860 Audio Codec 862 External Interface 864 Memory 866 Communication Interface 868 Transceiver 870 (Global Positioning System) Receiver Module 874 Expandable Memory 880 Mobile Phones 882 Smartphones 940 images 941 Reading pileup 942 Reference position 943 Base Quality Score 944 Mapping Confidence Score 961 Traditional Hypothesis 962 Traditional Hypothesis 963 Traditional Hypothesis 964 Non-conventional Hypothesis 965 Non-conventional Hypothesis 966 Non-conventional Hypothesis 967 Non-conventional Hypothesis 980 Probability score result column 1040 images 1041 Pileup of aligned sequence reads 1044 Mapping Confidence Score 1061 Conventional Hypothesis 1062 Conventional Hypothesis 1063 Conventional Hypothesis 1064 Non-conventional Hypothesis 1065 Non-conventional Hypothesis 1066 Non-conventional Hypothesis 1067 Non-conventional Hypothesis 1080 Probability score result column 1140 images 1141 Pileup of aligned sequence reads 1143 Base Quality Score 1145 Current reference position 1160 Complete set of modified probability results 1161 Conventional Hypothesis 1162 Conventional Hypothesis 1163 Conventional Hypothesis 1164 Non-conventional Hypothesis 1165 Non-conventional Hypothesis 1166 Non-conventional Hypothesis 1167 Non-conventional Hypothesis 1180 Probability score result column 1240 images 1241 Pileup of aligned sequence reads 1242 Reference position 1243 Base Quality Score 1260 Complete set of modified probability results 1261 Conventional Hypothesis 1262 Conventional Hypothesis 1263 Traditional Hypothesis 1264 Non-conventional Hypothesis 1265 Non-conventional Hypothesis 1266 Non-conventional Hypothesis 1267 Non-conventional Hypothesis 1280 Probability Score Results

Claims

1. 1. A method for improving the accuracy of variant calling by taking into account a measure of correlated error events, comprising: accessing, by one or more computers and from one or more memory devices, a pileup of a plurality of sequence reads aligned to a first region of a reference genome; obtaining, by the one or more computers, information describing one or more characteristics of each of the plurality of reads of the pileup that correspond to a first location of the reference genome; providing, by the one or more computers and based on the obtained information, one or more inputs to a probabilistic model describing the one or more characteristics of the plurality of readings of the pile-up, the probabilistic model being configured to determine, for each hypothesis of one or more hypotheses selected based on the one or more inputs, a score indicating whether the hypothesis is true; obtaining, by the one or more computers, output information for each of the one or more hypotheses, the output information for each of the one or more hypotheses (i) generated by the probabilistic model based on processing of the one or more inputs to the probabilistic model describing the one or more characteristics of each reading of the pile-up, and (ii) indicating a score indicating whether the hypothesis is true; determining a likelihood that a true variant is present at the first location based on the obtained output information generated by the one or more computers and by the probabilistic model for each of the plurality of hypotheses; A method comprising:

2. determining a likelihood that a true variant is present at the first location based on the obtained output information generated by the one or more computers and by the probabilistic model for each of the plurality of hypotheses; determining, by the one or more computers, a population score based on the output information generated by the probabilistic model for each of the plurality of hypotheses, the population score being indicative of a likelihood that the true variant is present; determining, by the one or more computers, whether the score generated by the population score meets a predetermined threshold; adding, by the one or more computers, information to a VCF file indicating that a true variant is present at the first position based on a determination that the population score meets the predetermined threshold; The method of claim 1 , comprising:

3. 3. The method of claim 2, wherein the information indicating that a true variant is present at the first position comprises information identifying (i) the first position, (ii) a candidate alternative allele at the first position, and (iii) the population score.

4. 2. The method of claim 1, wherein the information describing the one or more characteristics of each read comprises information describing (i) a mapping quality score for each sequence read of the pileup at the first position, and (ii) a read allele score for each sequence read of the pileup at the first position for each candidate allele at the first position.

5. The read allele score of each read of the pileup at the first position is a score for a particular candidate allele G ｍ,φ Given a reading r ｉ 5. The method of claim 4, wherein the probability of observing a first position is based on an output generated by a P-HMM model for each of the readings at the first position, the output being indicative of a probability of observing a first position.

6. The output information is a first output information regarding a first hypothesis of the one or more hypotheses, the first output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous reference with an exogenous allele matching the alternative; and and second output information regarding a second hypothesis of the one or more hypotheses, the second output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous substitution with an exogenous allele that matches a reference allele. The method of claim 4 , comprising:

7. 2. The method of claim 1, wherein the information describing one or more characteristics of each of the sequence reads comprises information describing: (i) a read orientation of each sequence read of the pileup at the first position, (ii) a position of each base at the first position within each sequence read of the pileup at the first position referenced to a 5' end of the sequence read, (iii) a read allele score of each sequence read of the plurality of sequence reads for each candidate allele at a reference position, and (iv) a base quality score of each read of the base at the first position.

8. The read allele score of each sequence read of the pileup at the first position is determined to be a particular candidate allele G ｍ,φ Given an array read r ｉ 8. The method of claim 7, wherein the probability of observing a sequence read at each of the first positions is based on an output generated by a P-HMM model, the P-HMM model indicating a probability of observing a sequence read at each of the first positions.

9. 2. The method of claim 1, wherein the information describing the one or more characteristics of each of the sequence reads comprises information describing: (i) a read orientation of each sequence read of the pileup at the first position, (ii) a position of each base at the first position within each sequence read of the pileup at the first position, referenced to a 5' end of the sequence read, and (iii) a read allele score of each sequence read of the plurality of sequence reads for each candidate allele at the first position.

10. The output information is a first output information regarding a first hypothesis of the one or more hypotheses, the first output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous reference with a sequencing error corresponding to the alternative allele; and and second output information regarding a second hypothesis of the one or more hypotheses, the second output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous substitution with a sequencing error matching the reference allele. The method of claim 7, comprising:

11. 2. The method of claim 1, wherein the information describing the one or more characteristics of each of the sequence reads comprises information describing: (i) a read orientation of each sequence read of the pileup at the first position; (ii) a position of each base at the first position, such as position "0" 142 within each sequence read, referenced to the 5' end of the sequence read; (iii) a mapping quality score of each sequence read of the pileup at the first position; (iv) a read allele score of each sequence read of the plurality of reads for each candidate allele at a reference position; and (v) a base quality score of each read of the aligned base at the first position.

12. The read allele score of each sequence read of the pileup at the first position is determined to be a particular candidate allele G ｍ,φ Given a reading r ｉ 12. The method of claim 11, wherein the probability of observing a sequence read at each of the first positions is based on an output generated by a P-HMM model, the P-HMM model indicative of a probability of observing a sequence read at each of the first positions.

13. The output information is a first likelihood that the sequence read at the first location indicates the occurrence of a homozygous reference with an exogenous allele that matches the alternative; a second likelihood that the sequence read at the first location indicates the occurrence of a homozygous substitution with an exogenous allele that matches a reference allele; a third output information regarding a first hypothesis of the one or more hypotheses, the third output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous reference with a sequencing error corresponding to the alternative allele; and the likelihood that the sequence read at the first position indicates the occurrence of a homozygous substitution with a sequencing error that matches the reference allele; The method of claim 11 , further comprising: fourth output information relating to a second hypothesis of the one or more hypotheses comprising:

14. 2. The method of claim 1, wherein the one or more memory devices receive the pile-up of aligned array reads from a field programmable gate array (FPGA) device, the FPGA including one or more configurable digital logic gates configured as a mapping and alignment unit for performing read mapping and alignment.

15. the computer is configured to access the one or more memory devices using one or more wired or wireless networks; a field programmable gate array (FPGA) device and the one or more memory devices are housed on an expansion card coupled to a circuit board of the sequencing device; the sequencing apparatus is configured to generate sequence reads based on input samples and store the generated sequence reads in the one or more memory devices; 15. The method of claim 14, wherein the mapping and alignment unit of the FPGA is configured to access the one or more memory devices to obtain the generated sequence reads.

16. the computer and the sequencing apparatus are each configured to access the one or more memory devices using one or more wired or wireless networks; the field programmable gate array (FPGA) device and the one or more memory devices are contained on an expansion card coupled to a circuit board of a server located remotely from the computer and the sequencing apparatus; the sequencing device is configured to generate sequence reads based on input samples and provide the generated sequence reads to the server using the one or more wired or wireless networks for storing the generated sequence reads in the one or more memory devices; 15. The method of claim 14, wherein the mapping and alignment unit of the FPGA is configured to access the one or more memory devices to obtain the generated sequence reads.

17. 1. A system comprising: one or more computers; and one or more storage devices storing instructions that, when executed by the one or more computers, cause the one or more computers to: accessing, by one or more computers and from one or more memory devices, a pileup of a plurality of sequence reads aligned to a first region of a reference genome; obtaining, by the one or more computers, information describing one or more characteristics of each of the plurality of reads of the pileup that correspond to a first location of the reference genome; providing, by the one or more computers and based on the obtained information, one or more inputs to a probabilistic model describing the one or more characteristics of the plurality of readings of the pile-up, the probabilistic model being configured to determine, for each hypothesis of one or more hypotheses selected based on the one or more inputs, a score indicating whether the hypothesis is true; obtaining, by the one or more computers, output information for each of the one or more hypotheses, the output information for each of the one or more hypotheses (i) generated by the probabilistic model based on processing of the one or more inputs to the probabilistic model describing the one or more characteristics of the readings of each of the pile-ups, and (ii) indicating a score indicating whether the hypothesis is true; determining a likelihood that a true variant is present at the first position based on the obtained output information generated by the one or more computers and by the probabilistic model for each of the plurality of hypotheses; and one or more storage devices operable to cause the system to perform operations including:

18. determining a likelihood that a true variant is present at the first location based on the obtained output information generated by the one or more computers and by the probabilistic model for each of the plurality of hypotheses; determining, by the one or more computers, a population score based on output information generated by the probabilistic model for each of the plurality of hypotheses, the population score being indicative of a likelihood that the true variant is present; and determining, by the one or more computers, whether the score generated by the population score meets a predetermined threshold; adding, by the one or more computers, information to a VCF file indicating that a true variant is present at the first position based on a determination that the population score meets the predetermined threshold; 20. The system of claim 19, comprising:

19. 20. The system of Claim 18, wherein the information indicating that a true variant is present at the first position comprises information identifying: (i) the first position; (ii) a candidate alternative allele at the first position; and (iii) the population score.

20. 20. The system of claim 17, wherein the information describing the one or more characteristics of each read comprises information describing (i) a mapping quality score for each sequence read of the pileup at the first position, and (ii) a read allele score for each sequence read of the pileup at the first position for each candidate allele at the first position.

21. The read allele score of each read of the pileup at the first position is a score for a particular candidate allele G ｍ,φ Given a reading r ｉ 21. The system of claim 20, wherein the probability of observing a first position is based on an output generated by a P-HMM model for each of the readings at the first position, the output being indicative of a probability of observing a first position.

22. The output information is a first output information regarding a first hypothesis of the one or more hypotheses, the first output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous reference with an exogenous allele matching the alternative; and 21. The method of claim 20, further comprising: second output information regarding a second hypothesis of the one or more hypotheses comprising a probability that the sequence read at the first position indicates the occurrence of a homozygous substitution with an exogenous allele that matches a reference allele.

23. 18. The system of claim 17, wherein the information describing the one or more characteristics of each sequence read comprises information describing: (i) a read orientation of each sequence read of the pileup at the first position, (ii) a position of each base at the first position within each sequence read of the pileup at the first position referenced to a 5' end of the sequence read, (iii) a read allele score of each sequence read of the plurality of sequence reads for each candidate allele at a reference position, and (iv) a base quality score of each read of the base at the first position.

24. The read allele score of each sequence read of the pileup at the first position is determined to be a particular candidate allele G ｍ,φ Given an array read r ｉ 24. The system of claim 23, wherein the probability of observing a sequence read at each of the first positions is based on an output generated by a P-HMM model, the P-HMM model indicative of a probability of observing a sequence read at each of the first positions.

25. 18. The system of claim 17, wherein the information describing the one or more characteristics of each sequence read comprises information describing: (i) a read orientation of each sequence read of the pileup at the first position, (ii) a position of each base at the first position within each sequence read of the pileup at the first position referenced to the 5' end of the sequence read, and (iii) a read allele score of each sequence read of the plurality of sequence reads for each candidate allele at the first position.

26. The output information is a first output information regarding a first hypothesis of the one or more hypotheses, the first output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous reference with a sequencing error corresponding to the alternative allele; and and second output information regarding a second hypothesis of the one or more hypotheses, the second output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous substitution with a sequencing error matching the reference allele.

26. The system of claim 25, comprising:

27. 18. The system of claim 17, wherein the information describing the one or more characteristics of each sequence read comprises information describing: (i) a read orientation of each sequence read of the pileup at the first position; (ii) a position of each base at the first position, such as position "0" 142 within each sequence read, referenced to the 5' end of the sequence read; (iii) a mapping quality score of each sequence read of the pileup at the first position; (iv) a read allele score of each sequence read of the plurality of reads for each candidate allele at the reference position; and (v) a base quality score of each read of the aligned base at the first position.

28. The read allele score of each sequence read of the pileup at the first position is determined to be a particular candidate allele G ｍ,φ Given a reading r ｉ 28. The system of claim 27, wherein the probability of observing a sequence read at each of the first positions is based on an output generated by a P-HMM model, the P-HMM model indicative of a probability of observing a sequence read at each of the first positions.

29. The output information is a first likelihood that the sequence read at the first location indicates the occurrence of a homozygous reference with an exogenous allele that matches the alternative; a second likelihood that the sequence read at the first location indicates the occurrence of a homozygous substitution with an exogenous allele that matches a reference allele; a third output information regarding a first hypothesis of the one or more hypotheses, the third output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous reference with a sequencing error corresponding to the alternative allele; and and fourth output information regarding a second hypothesis of the one or more hypotheses, the fourth output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous substitution with a sequencing error matching the reference allele.

28. The system of claim 27, comprising:

30. 20. The system of claim 17, wherein the one or more memory devices receive the pile-up of aligned array reads from a field programmable gate array (FPGA) device, the FPGA including one or more configurable digital logic gates configured as a mapping and alignment unit for performing read mapping and alignment.

31. the computer is configured to access the one or more memory devices using one or more wired or wireless networks; a field programmable gate array (FPGA) device and the one or more memory devices are housed on an expansion card coupled to a circuit board of the sequencing device; the sequencing apparatus is configured to generate sequence reads based on input samples and store the generated sequence reads in the one or more memory devices; 31. The system of claim 30, wherein the mapping and alignment unit of the FPGA is configured to access the one or more memory devices to obtain the generated sequence reads.

32. the computer and the sequencing apparatus are each configured to access the one or more memory devices using one or more wired or wireless networks; the field programmable gate array (FPGA) device and the one or more memory devices are contained on an expansion card coupled to a circuit board of a server located remotely from the computer and the sequencing apparatus; the sequencing device is configured to generate sequence reads based on input samples and provide the generated sequence reads to the server using the one or more wired or wireless networks for storing the generated sequence reads in the one or more memory devices; 31. The system of claim 30, wherein the mapping and alignment unit of the FPGA is configured to access the one or more memory devices to obtain the generated sequence reads.

33. A non-transitory computer-readable medium storing software including instructions executable by one or more computers, which, upon such execution, cause said one or more computers to: accessing, by the one or more computers and from one or more memory devices, a pileup of a plurality of sequence reads aligned to a first region of a reference genome; obtaining, by the one or more computers, information describing one or more characteristics of each of the plurality of reads of the pileup that correspond to a first location of the reference genome; providing, by the one or more computers and based on the obtained information, one or more inputs to a probabilistic model describing the one or more characteristics of the plurality of readings of the pile-up, the probabilistic model being configured to determine, for each hypothesis of one or more hypotheses selected based on the one or more inputs, a score indicating whether the hypothesis is true; obtaining, by the one or more computers, output information for each of the one or more hypotheses, the output information for each of the one or more hypotheses (i) generated by the probabilistic model based on processing of the one or more inputs to the probabilistic model describing the one or more characteristics of the respective readings of the pile-up, and (ii) indicating a score indicating whether the hypothesis is true; determining a likelihood that a true variant is present at the first position based on the obtained output information generated by the one or more computers and by the probabilistic model for each of the plurality of hypotheses; A non-transitory computer-readable medium for causing operations to be performed, including:

34. determining a likelihood that a true variant is present at the first location based on the obtained output information generated by the one or more computers and by the probabilistic model for each of the plurality of hypotheses; determining, by the one or more computers, a population score based on output information generated by the probabilistic model for each of the plurality of hypotheses, the population score being indicative of a likelihood that the true variant is present; and determining, by the one or more computers, whether the score generated by the population score meets a predetermined threshold; adding, by the one or more computers, information to a VCF file indicating that a true variant is present at the first position based on a determination that the population score meets the predetermined threshold; 34. The computer readable medium of claim 33, comprising:

35. 35. The computer-readable medium of Claim 34, wherein the information indicating that a true variant is present at the first position comprises information identifying (i) the first position, (ii) a candidate alternative allele at the first position, and (iii) the population score.

36. 34. The computer-readable medium of claim 33, wherein the information describing the one or more characteristics of the respective reads comprises information describing (i) a mapping quality score for each sequence read of the pileup at the first position, and (ii) a read allele score for each sequence read of the pileup at the first position for each candidate allele at the first position.

37. The read allele score of each read of the pileup at the first position is a score for a particular candidate allele G ｍ,φ Given a reading r ｉ 37. The computer-readable medium of claim 36, wherein the probability of observing a first position is based on an output generated by a P-HMM model for each of the readings at the first position, the output being indicative of a probability of observing a first position.

38. The output information is a first output information regarding a first hypothesis of the one or more hypotheses, the first output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous reference with an exogenous allele matching the alternative; and and second output information regarding a second hypothesis of the one or more hypotheses, the second output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous substitution with an exogenous allele that matches a reference allele.

37. The computer readable medium of claim 36, comprising:

39. 34. The computer-readable medium of claim 33, wherein the information describing the one or more characteristics of each sequence read comprises information describing: (i) a read orientation of each sequence read of the pileup at the first position, (ii) a position of each base at the first position within each sequence read of the pileup at the first position referenced to a 5' end of the sequence read, (iii) a read allele score of each sequence read of the plurality of sequence reads for each candidate allele at a reference position, and (iv) a base quality score of each read of the base at the first position.

40. The read allele score of each sequence read of the pileup at the first position is determined to be a particular candidate allele G ｍ,φ Given an array read r ｉ 40. The computer-readable medium of claim 39, wherein the probability of observing a sequence read at each of the first positions is based on an output generated by a P-HMM model, the P-HMM model indicating a probability of observing a sequence read at each of the first positions.

41. 40. The computer-readable medium of claim 39, wherein the information describing the one or more characteristics of each sequence read comprises information describing: (i) a read orientation of each sequence read of the pileup at the first position, (ii) a position of each base at the first position within each sequence read of the pileup at the first position referenced to the 5' end of the sequence read, and (iii) a read allele score of each sequence read of the plurality of sequence reads for each candidate allele at the first position.

42. The output information is a first output information regarding a first hypothesis of the one or more hypotheses, the first output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous reference with a sequencing error corresponding to the alternative allele; and 40. The computer-readable medium of claim 39, further comprising second output information regarding a second hypothesis of the one or more hypotheses comprising a probability that the sequence read at the first position indicates the occurrence of a homozygous substitution with a sequencing error matching the reference allele.

43. 34. The computer-readable medium of claim 33, wherein the information describing the one or more characteristics of the respective sequence reads comprises information describing: (i) a read orientation of each sequence read of the pileup at the first position; (ii) a position of each base at the first position, such as position "0" 142 within each sequence read, referenced to the 5' end of the sequence read; (iii) a mapping quality score of each sequence read of the pileup at the first position; (iv) a read allele score of each sequence read of the plurality of reads for each candidate allele at the reference position; and (v) a base quality score of each read of the aligned base at the first position.

44. The read allele score of each sequence read of the pileup at the first position is determined to be a particular candidate allele G ｍ,φ Given a reading r ｉ 44. The computer-readable medium of claim 43, wherein the probability of observing a sequence read at each of the first positions is based on an output generated by a P-HMM model, the P-HMM model indicating a probability of observing a sequence read at each of the first positions.

45. The output information is a first likelihood that the sequence read at the first location indicates the occurrence of a homozygous reference with an exogenous allele that matches the alternative; a second likelihood that the sequence read at the first location indicates the occurrence of a homozygous substitution with an exogenous allele that matches a reference allele; a third output information regarding a first hypothesis of the one or more hypotheses, the third output information including a probability that the sequence read at the first position indicates the occurrence of a homozygous reference with a sequencing error corresponding to the alternative allele; and 44. The computer-readable medium of claim 43, further comprising fourth output information regarding a second hypothesis of the one or more hypotheses comprising a probability that the sequence read at the first position indicates the occurrence of a homozygous substitution with a sequencing error matching the reference allele.

46. 34. The computer-readable medium of claim 33, wherein the one or more memory devices receive the pile-up of aligned sequence reads from a field programmable gate array (FPGA) device, the FPGA including one or more configurable digital logic gates configured as a mapping and alignment unit for performing read mapping and alignment.

47. the computer is configured to access the one or more memory devices using one or more wired or wireless networks; a field programmable gate array (FPGA) device and the one or more memory devices are housed on an expansion card coupled to a circuit board of the sequencing device; the sequencing apparatus is configured to generate sequence reads based on input samples and store the generated sequence reads in the one or more memory devices; 37. The computer-readable medium of claim 36, wherein the mapping and alignment unit of the FPGA is configured to access the one or more memory devices to obtain the generated sequence reads.

48. the computer and the sequencing apparatus are each configured to access the one or more memory devices using one or more wired or wireless networks; the field programmable gate array (FPGA) device and the one or more memory devices are contained on an expansion card coupled to a circuit board of a server located remotely from the computer and the sequencing apparatus; the sequencing device is configured to generate sequence reads based on input samples and provide the generated sequence reads to the server using the one or more wired or wireless networks for storing the generated sequence reads in the one or more memory devices; 37. The computer-readable medium of claim 36, wherein the mapping and alignment unit of the FPGA is configured to access the one or more memory devices to obtain the generated sequence reads.

Images