Systems and methods for germline SNV and indel variant calling

Abstract

System and methods are described for performing germline variant calling on non-Illumina generated sequencing data, such as nanopore sequencing data. As a staring point, existing variant callers developed and optimized for use with Illumina sequencing data can be used to generate an initial set of variant calls. As these variant calls are not optimized for other types of sequencing data, such as sequencing data generated by a sequencing-by-expansion (SBX) chemistry, further steps are taken to filter / regenotype the variant calls to eliminate false positives. A machine learning algorithm is trained based on features extracted from the variant call information to remove the false positives, thereby improving the germline variant caller as used with SBX-based sequencing data.

Description

PATENTClient Reference No.: P39561-WO-ISYSTEMS AND METHODS FOR GERMLINE SNV AND INDEL VARIANT CALLINGCROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This Application claims priority to United States Provisional Patent Application No. 63 / 683,191, filed on August 14, 2024, and United States Provisional Patent Application No. 63 / 729,228, filed on December 6, 2024, each of which is herein incorporated by reference in its entirety.BACKGROUND

[0002] The development of affordable and rapid molecular sequencing technologies has enabled the development of targeted therapeutics that rely on the use of deoxyribonucleic acid (DNA) biomarkers to identify patients that are suitable for receiving targeted therapies. For example, mutations in certain genes, such as genes involved in cell proliferation, are known to lead to certain types of cancers that can be treated very effectively with specific types of drugs. Other mutations are known to confer resistance to certain therapies. Some genetic variants are associated with hereditary diseases such as cystic fibrosis and certain cancers that may be passed on to future generations. Therefore, there is a need for improved systems and methods to identify germline variants from sequencing data.SUMMARY

[0003] The present disclosure provides systems and methods for performing variant calling on sequencing data, and more specifically to calling germline single nucleotide variants, insertions, and deletions in the sequencing data.

[0004] In accordance with a first aspect of the present disclosure, a method for germline variant calling is provided. The method can include receiving sequencing data generated by a sequencer instrument performing a sequencing assay; generating a plurality of consensus reads from the sequencing data; aligning the plurality of consensus reads to a reference sequence to generate alignment data; identifying at least one callable region of the sequencing data; identify ing, for each callable region in the at least one callable region using a germline variant caller, variant information for a plurality of observed variants in the callable region based on the alignment data; extracting a plurality of features from the variant information; and training a machine learning (ML) model to predict a corrected genotype forPATENTClient Reference No.: P39561-WO-I at least one variant of the plurality of observed variants based on the extracted plurality of features.

[0005] In accordance with an embodiment of the first aspect, the germline variant caller comprises a GATK Haplotype Caller. Optimized settings for the GATK Haplotype Caller comprise at least one of: mapping-quality -threshold-for-genotyping 1; minimum-mapping- quality 1 ; enable-dynamic-read-disqualification-for-genotyping; min-base-quali ty-score 0; adaptive-pruning true; and smith- waterman FASTEST AV AIL ABLE.

[0006] In accordance with an embodiment of the first aspect, the plurality of features includes at least one feature of each of the following: one or more read alignment features from the alignment data; one or more variant call features from the variant information; one or more reference context features based on the reference sequence; and one or more features extracted from a public dataset. In some embodiments, the public dataset comprises the gnomAD public database.

[0007] In accordance with an embodiment of the first aspect, the plurality of features for an observed variant at a specific variant site comprise at least one of: distribution of quality scores among reads that support a reference allele at the variant site; distribution of quality scores among reads that support an alternative allele at the variant site; a type of the observed variant; a count of reads that support at least one of the reference allele or the alternative allele; a distance of the variant site from the end of each read, or a mean, minimum, maximum, or sum of distances of the variant site from a plurality of reads that support a particular allele; a number of other variant calls within a given distance of the variant site; an allele frequency of the alternative allele based on the public dataset; a number of bases immediately preceding or immediately following the variant site; and / or a depth of coverage for the alternative allele.

[0008] In accordance with an embodiment of the first aspect, the ML model is a gradient boosting machine (GBM) model. In some embodiments, the ML model is lightGBM, and the ML model is configured with the following parameters: (1) 'objective': 'multiclass'; (2) 'boosting_type': 'gbdf; (3) 'metric': 'multi_logloss'; (4) 'leaming_rate': 0.01; (5) 'nrounds': 2000; (6) 'bagging_fraction': 0.9; (7) Teature_fraction': 0.5; (8) 'bagging_freq': 1; (9) 'num leaves': 64; and (10) 'min data in leaf : 20.PATENTClient Reference No.: P39561-WO-1

[0009] In accordance with an embodiment of the first aspect, the method further includes: applying the trained machine learning model to a second set of variant information generated by the germline variant caller to re-genotype observed variants in the second set of variant information that were incorrectly genotyped by the germline variant caller; generating a set of filtered germline / regenotyped variant calls from the second set of variant information; and generating an updated genotype for the sample based on the set of filtered germline / regenotyped variant calls. In some embodiments, the set of filtered germline / regentotyped variant calls comprises at least one single nucleotide variant (SNV). In some embodiments, the observed variants comprise at least one InDei.

[0010] In accordance with an embodiment of the first aspect, the sequencing data comprises sequencing-by-expansion (SBX) sequencing data measured by a nanopore-based sequencer instrument. In some embodiments, the germline variant caller is optimized for non- SBX sequencing data, e.g., Illumina-based sequencing data.

[0011] In accordance with a second aspect of the present disclosure, a system for germline variant calling is provided. The system includes: a memory storing sequencing data generated by a sequencer instrument performing a sequencing assay; and at least one processor. The at least one processor is configured to: generate a plurality of consensus reads from the sequencing data, align the plurality of consensus reads to a reference sequence to generate alignment data, identify at least one callable region of the sequencing data, identify, for each callable region in the at least one callable region using a germline variant caller, variant information for a plurality of observed variants in the callable region based on the alignment data, extract a plurality7of features from the variant information, and train a machine learning (ML) model to predict a corrected genotype for at least one variant of the plurality of observed variants based on the extracted plurality of features.

[0012] In accordance with an embodiment of the second aspect, the system further includes the sequencer instrument communicatively coupled to the at least one processor. In some embodiments, the sequencer instrument is a nanopore-based sequencer instrument configured to perform a sequencing-by-expansion (SBX) sequencing assay.

[0013] In accordance with an embodiment of the second aspect, the plurality of features includes at least one feature of each of the following: one or more read alignment features from the alignment data; one or more variant call features from the variant information; onePATENTClient Reference No.: P39561-WO-1 or more reference context features based on the reference sequence; and one or more features extracted from a public dataset.

[0014] In accordance with an embodiment of the second aspect, the one or more features for an observed variant at a specific variant site comprise at least one of: distribution of quality scores among reads that support a reference allele at the variant site; distribution of quality scores among reads that support an alternative allele at the variant site; a type of the observed variant; a count of reads that support at least one of the reference allele or the alternative allele; a distance of the variant site from the end of each read, or a mean, minimum, maximum, or sum of distances of the variant site from a plurality of reads that support a particular allele; a number of other variant calls within a given distance of the variant site; an allele frequency of the alternative allele based on the public dataset; a number of bases immediately preceding or immediately following the variant site; and / or a depth of coverage for the alternative allele.

[0015] In accordance with an embodiment of the second aspect, the ML model is a gradient boosting machine (GBM) model.

[0016] In accordance with a third aspect of the present disclosure, a non-transitory computer-readable medium is provided. The non-transitory computer-readable medium stores instructions that, executed by at least one processor, cause a computer device to: receive sequencing data generated by a sequencer instrument performing a sequencing assay; generate a plurality of consensus reads from the sequencing data; align the pl urality of consensus reads to a reference sequence to generate alignment data; identify at least one callable region of the sequencing data; identify, for each callable region in the at least one callable region using a germline variant caller, variant information for a plurality of observed variants in the callable region based on the alignment data; extract a plurality of features from the variant information; and train a machine learning (ML) model to predict a corrected genotype for at least one variant of the plurality of observed variants based on the extracted plurality of features.BRIEF DESCRIPTION OF THE DRAWINGSPATENTClient Reference No.: P39561-WO-I

[0017] The novel features of the invention are set forth with particularity in the claims that follow. A better understanding of the features and advantages described in the present disclosure will be obtained by reference to the accompanying drawings, included as follows.

[0018] FIG. 1 sets forth an illustrative system including a sequencing device communicatively coupled to a computing system, in accordance with at least some embodiments of the present disclosure.

[0019] FIG. 2 is an illustration of the workflow for a method of calling germline SNVs from sequencing data using a variant caller, in accordance with at least some embodiments of the present disclosure.

[0020] FIG. 3 is an illustration of the workflow for a method of calling germline InDeis from sequencing data using a variant caller, in accordance with at least some embodiments of the present disclosure.

[0021] FIG. 4 is an illustration of various criteria for excluding certain reads for use in the method of calling germline InDeis illustrated in FIG. 2, in accordance with at least some embodiments of the present disclosure..

[0022] FIG. 5 illustrates incorporating lossless encoding to assign partial support from R1 and R2 to improve InDei calling performance, in accordance with at least some embodiments of the present disclosure..

[0023] FIGS. 6 A & 6B show flow charts of a method for performing germline variant calling using a machine learning model to correct genotype information generated by a germline variant caller, in accordance with at least some embodiments.

[0024] FIG. 7 illustrates an example computer system that may be utilized to implement techniques disclosed herein.DETAILED DESCRIPTION

[0025] System and methods are described for performing germline variant calling on sequencing data, such as sequencing data called using a nanopore-based sequencer instrument. Existing germline variant callers, developed and optimized for use with sequencing data generated by industry leading Illumina, Inc. sequencers, can be used toPATENTClient Reference No.: P39561-WO-1 generate an initial set of variant calls. These variant calls may not have the same accuracy when the germline variant caller is applied to sequencing data from other types of sequencer instruments, which can exhibit different error rates compared to sequencers from Illumina, Inc. The germline variant caller is not optimized for other types of sequencing data, such as sequencing data generated by a sequencing-by-expansion (SBX) chemistry. Therefore, further steps are taken to filter the variant calls generated by the existing germline variant caller to eliminate false positives.

[0026] In an embodiment, a machine learning algorithm, such as a gradient boosting machine (GBM), is trained based on features extracted from the variant call information to remove the false positives, thereby improving the germline variant caller as used with SBX- based sequencing data. The features can also include alignment features, reference features, and features from public datasets.Sequencing

[0027] Prepared nucleic acid molecules of interest (e.g., a sequencing library) can be sequenced using a sequencing assay as part of the procedure for determining sequencing reads for a plurality of genome positions (e.g., loci) or as part of a whole genome sequencing application. Any of a number of sequencing technologies or sequencing assays can be utilized. The term "Next Generation Sequencing (NGS)" as used herein refers to sequencing methods that allow for massively parallel sequencing of clonally amplified molecules and of single nucleic acid molecules (or of nucleic acid analogues).

[0028] Non-limiting examples of sequence assays that are suitable for use with the methods disclosed herein include nanopore sequencing (US Pat. Appl. Publ. Nos.2013 / 0244340, 2013 / 0264207, 2014 / 0134616, 2015 / 0119259 and 2015 / 0337366), Sanger sequencing, capillary array sequencing, thermal cycle sequencing (Sears et al., Biotechniques. 13:626-633 (1992)), solid-phase sequencing (Zimmerman et al., Methods Mol. Cell Biol., 3:39-42 (1992)), sequencing with mass spectrometry such as matrix-assisted laser desorption / ionization time-of-flight mass spectrometry (MALDI-TOF / MS; Fu et al., Nature Biotech., 16:381-384 (1998)), sequencing by hybridization (Drmanac et al., Nature Biotech., 16:54-58 (1998), and NGS methods, including but not limited to sequencing bysynthesis (e.g.. HiSeq™. MiSeq™, or Genome Analyzer, each available from Illumina),PATENTClient Reference No.: P39561-WO-I sequencing by ligation (e.g., SOLiD™, Life Technologies), ion semiconductor sequencing (e.g., Ion Torrent™, Life Technologies), and SMRT® sequencing (e.g., Pacific Biosciences).

[0029] Commercially available sequencer instruments include: sequencing-by- hybridization platforms from Affymetrix, Inc., now ThermoFisher Scientific, (Sunnyvale, Calif), sequencing-by-synthesis platforms from Illumina, Inc. (San Diego, Calif.) and Helicos Biosciences (Cambridge, Mass.), and sequencing-by-ligation platform from Applied Biosystems (Foster City. Calif). Other sequencing technologies include, but are not limited to, the Ion Torrent technology from ThermoFisher Scientific (Waltham, MA); as well as nanopore sequencing from Roche Sequencing Solutions, Inc. (Santa Clara, Cal.) and Oxford Nanopore Technologies (Oxford, United Kingdom).Sequencing Systems

[0030] FIG. 1 sets forth an illustrative system 100 including a sequencing device 110 (which may be alternately referred to as a sequencer instrument) communicatively coupled to a computing system 102. Sequencing device 110 can be coupled to computing system 102 either directly (e.g., through one or more communication cables) or through network 130, which may be the Internet or any other combination of wide-area, local area, wired, and / or wireless netw orks. In some embodiments, computing system 102 may be included in or integrated with the sequencing device 110. In some embodiments, sequencing device 110 may sequence (e g., perform a biochemical assay) a sample containing genetic material and produce resulting sequencing data. The sequencing data can be sent to computing system 102 (e.g., through network 130) or stored on a storage device and at a later stage transferred to computing system 102 (e.g., through network 130). In some embodiments, computing system 102 may or may not include a display 108 and one or more input devices (not illustrated) for receiving commands from a user or operator (e.g. a technician or a geneticist). In some embodiments, computing system 102 and / or sequencing device 110 can be accessed by users or other devices remotely through network 130. Thus, in some embodiments various methods discussed herein may be run remotely on computing system 102.

[0031] Computing system 102 may include one computing device or a combination of a number of computing devices of any type, such as personal computers, laptops, network servers (e.g.. local servers or servers included on a public / private / hybrid cloud), mobile devices, etc., where some or all of the devices can be interconnected. Computing system 102PATENTClient Reference No.: P39561-WO-I may include one or more processors (not illustrated), each of which can have one or more logic cores. In some embodiments, computing system 102 can include one or more general- purpose processors (e.g., CPUs), special-purpose processors such as graphics processors (GPUs), digital signal processors, or any combination of these and other types of processors. In some embodiments, some or all processors in computing system can be implemented using customized or customizable circuitry, such as application specific integrated circuits (ASICs) or field programmable gate arrays (FPGAs). Computing system 102 can also in some embodiments retrieve and execute non-transitory computer-readable instructions stored in one or more memories or storage devices (not illustrated) integrated into or otherwise communicatively coupled to computing system 102. The memory / storage devices can include any combination of non-transiloiy computer readable storage media including semiconductor memory chips of various ty pes (e.g.. dynamic random access memory (DRAM), static random access memory (SRAM), synchronous dynamic random access memory (SDRAM), flash memory, programmable read-only memory, etc.) and so on. Magnetic and / or optical disks can also be used. The memories / storage devices can also include removable storage media that can be readable and / or writeable; examples of such media include compact disc (CD), read-only digital versatile disc (e.g., DVD-ROM. dual-layer DVD-ROM), read-only and recordable Blu-ray® disks, ultra-density optical disks, flash memory cards (e.g., SD cards, mini-SD cards, micro-SD cards, etc.), and so on. In some embodiments, data and other information (e.g. sequencing data) can be stored in one or more remote locations, e.g., cloud storage, and synchronized with other the components of system 100.

[0032] In some embodiments, the sequencing device 110 can generate sequencing data by a sequencing by expansion (SBX) process. Examples of the SBX process include those described in U.S. Patent Application No. 17 / 456,342 (U.S. Publication No.US20220411458A1), entitled “Translocation control elements, reporter codes, and further means for translocation control for use in nanopore sequencing,” fded November 23, 2021 , which is herein incorporated by reference in its entirety. During library7preparation in the SBX process, a number of surrogate molecules are derived from and characterize nucleic acid material provided in a sample.

[0033] More particularly, the SBX process may translate a sequence of DNA into a measurable surrogate molecule called an Xpandomer. Xpandomer synthesis based on the natural function of DNA replication uses expandable nucleotide triphosphates (X-NTPs) thatPATENTClient Reference No.: P39561-WO-1 act as substrates for template-dependent, polymerase-based replication. These Xpandomer molecules are then processed by a sequencer instrument (e.g., the sequencing device 110) to measure the sequence in the original DNA template. As the Xpandomer molecule transits through nanometer-sized openings in an electrode-resistant membrane (a “nanopore”), each nanopore corresponding to a selective channel, a distinct electrical signal is generated for each base reporter and identifiable to enable highly accurate and high throughput nanoporebased nucleic acid sequencing (also referred to generally as nanopore sequencing).

[0034] Sequencing device 110 can generate a plurality of sequence reads corresponding to a genetic sample (e.g., a sample comprising a patient’s DNA or RNA material). For example, a sequence may be identified by processing a sample, which may include (for example) a blood, saliva, or tissue biopsy collected from a subject. Sequence reads can be obtained either directly from sequencing device 110, or from one or more local or remote volatile or non-volatile memories, storage devices, or databases communicatively coupled to computing system 102. Sequence reads can be pre-processed (e.g., pre-aligned) or they can be “raw,” in which case a downstream method may include a preprocessing (e.g., prealigning) step. Also, while in some embodiments, entire sequence reads (as generated by sequencing device 110) can be obtained, in other embodiments only sections of the sequence reads can be obtained. Thus, “obtaining a sequence read,” as used herein, refers generally to obtaining one or more sections of one or more (e.g., adjacent) sequence reads.Bioinformatics Workflow Overview

[0035] The output of an NGS sequencer is generally processed by a bioinformatics pipeline that processes the raw signal from the NGS sequencer and translates the raw signal into base calls, often referred to as raw reads, which are typically stored in a FASTQ file that combines the raw reads with associated quality data. This portion of the bioinformatics pipeline is often referred to as primary analysis.

[0036] The next section of the bioinformatics pipeline is called secondary analysis, and it takes the raw reads generated by the primary analysis, and performs several tasks, including alignment and variant calling.PATENTClient Reference No.: P39561-WO-1

[0037] Tertiary analysis is the final portion of the bioinformatics pipeline and uses the variant calling information to general medical insights that health care practitioners can use to improve treatments for their patients.Secondary Analysis

[0038] New sequencing technologies, such as nanopore-based sequencers, generate sequencing data with different characteristics than sequencing data generated by the current market leading sequencers, such as Illumina sequencers. For example, these differences can include differences in raw read accuracy and differences in the error profiles. Because Illumina sequencers currently dominate the market, the vast majority of the secondaryanalysis software tools that have been developed are custom tailored to process the type of data that is generated by the Illumina sequencers. These software tools, which typically work very well with data from Illumina sequencers, may not work well with data generated by new next generation sequencing technologies, such as nanopore sequencers. Consequently, there is a need to develop new secondary analysis tools that work well with the new sequencing technologies that are currently being developed. In addition, although the variant calling methods described herein may be particularly effective with nanopore sequencing data, the methods can also be used wi th other types of sequencing data.Germline SNV and Multi-Allelic Variant Calling

[0039] FIG. 2 illustrates an embodiment of a secondary analysis bioinformatics pipeline for calling germline single nucleotide variants (SNVs) and multi-allelic variants from sequencing data, and particularly from nanopore based sequencing data. The method includes five main steps:(1) Alignment of consensus reads to a reference genome;(2) Determination of active regions (sites where the alignments show signs of variants);(3) Assembly of possible haplotypes within each active region and haplotype- informed read re-alignment;(4) Identification of candidate variants based on assembled haplotype alignments; and(5) Genotyping variants and removing false positives.PATENTClient Reference No.: P39561-WO-I

[0040] The method starts with consensus reads 200 from, for example, a nanopore sequencer instrument 110 or another type of sequencer device. A ty pical format for receiving these reads is in a file using the FASTQ format. The consensus reads 200 stored in the file can be accessed and processed by the pipeline to perform an alignment 202 of the consensus reads 200 against a reference sequence. The aligned sequences can be stored in a SAM or BAM file that can be sent to a germline variant caller 204, such as GATK HaplotypeCaller, for example.

[0041] The germline variant caller 204, such as the GATK HaplotypeCaller, performs a variety of steps in order to call germline variants. First, a callable region can be identified from the aligned sequences as a region with a sufficient depth of read coverage to allow for variant calling. For example, the callable region can have a minimum depth of read coverage of at least 5, 10, 15, 20, 25, 30, 35, 40. 45. or 50. In some embodiments, the required depth of coverage can vary depending on a variety of factors, such as the types of bases in the region (i.e., GC rich regions or other type of motif which may have a higher error rate during sequencing would benefit from higher depth of coverage), and the quality scores of the bases in the region.

[0042] Next, within the callable regions, mismatches or gaps between the base calls of the reads and the reference sequence are used to identity' active regions that contain potential SNVs (single nucleotide variants) and InDeis (insertions and deletions).

[0043] For each active region, an assembly graph, such as a Debruijn graph for example, is generated by starting with the reference sequence, which is decomposed into a series of k- mers (short sequences k bases long), with each successive k-mer overlapping the previous k- mer by k-1 bases. The k-mers can be represented as nodes that can be joined by lines called edges. The edges can be weighted to keep track of the number of k-mers found in the sample, with the weights initially set to zero. This forms a reference graph.

[0044] Next, each sequence read can similarly be decomposed into a series of k-mers and can be matched to the reference graph. Each time two successive k-mers are matched to the graph, the weight of the edge joining the two k-mer nodes is incremented by one. If a k-mer cannot be matched to the graph, a new node and edge is added. This is repeated for all the k- mers in the sequence read. Under this method, the weights indicate the number of times aPATENTClient Reference No.: P39561-WO-1 particular k-mer was found in the sequence reads. In addition, the weights can also be used to determine the most likely path(s) through the graph.

[0045] In some embodiments, the assembly graph can then be optionally pruned by removing sections of the graph that are supported by an edge weight that is fewer than a threshold value, such as 2. A weight of 2, for example, means that 2 reads in the sample support that segment of the graph. The threshold can be increased, for more aggressive pruning, which will result in faster processing and higher specificity, but with lower sensitivity. Decreasing the threshold, means less pruning, which lowers the specificity, but increases the sensitivity. In some embodiments, the threshold can be set to zero, which means that the graph is not pruned for maximum sensitivity'.

[0046] After the graph has been assembled and optionally pruned, haplotype sequences can be generated from the graph by traversing all paths in the graph, with a likelihood score calculated for each haploty pe sequence as the product of the transition probabilities of the path edges. The probability of an edge can be calculated as the weight of the edge divided by the sum of the weights of all the edges that share the same source node. The haplotypes with the highest likelihood scores can be used for candidate variant detection.

[0047] These haplotypes are then aligned to the reference sequence in order to identify candidate variants. The alignment can be done using a Smith-Waterman alignment (SWA), for example, and can be stored in a SAM or BAM fde. The candidate variants can be determined by comparing the aligned sequence to the reference sequence and can be stored in a VCF (variant call format) file.

[0048] The steps described above can be performed, for example, using optimized settings for GATK HaplotypeCaller, an open source germline variant caller 204. The optimized settings include: — mapping-quality-threshold-for-genotyping 1 -RF MappingQualityReadFilter -minimum-mapping-quality 1 -enable-dynamic-read- disqualification-for-genotyping -min-base-quality-score 0 -adaptive-pruning true -smith- waterman FASTEST AVAILABLE. These parameters are changed from the defaults to take into account the differences in quality7scores and error modes in sequencing data generated by a nanopore sequencer as compared to industry standard Illumina sequencers, while balancing high sensitivity with fast computation times. The optimal parameters can bePATENTClient Reference No.: P39561-WO-1 identified by maximizing sensitivity based on a grid search, -min-base-quality -score 1 can be used instead to give a faster run time.

[0049] The candidate variants generated by the method described above has a high sensitivity (i.e., at least 90, 91, 92, 93, 94, 95, 96, 97, 98, or 99%), meaning most of the actual variants in the sample have been called, but with a relatively high number of false positive variant calls, meaning some of the candidate variants are in fact not actual variants. Selecting parameters so that the GATK HaplotypeCaller (or another variant caller) has a high sensitivity increases or maximizes the number of variant calls generated by the variant caller, which is desirable since the method described herein includes an additional filtering step to remove false positives. At least part of the reason for the high number of false positive variant calls is a result of the differences between the nanopore sequencing data and the Illumina sequencing data, as mentioned above with respect to the differences in uality scores and error modes. Therefore, to reduce the number of false positive variant calls, an additional filtering and re-genotyping step 208 can be added to the pipeline.

[0050] The additional filtering and re-genotyping step 208 can be based on a machine learning algorithm that filters and re-genotypes the germline small variant calls and initial genotype that result from running GATK HaplotypeCaller with the settings described above.

[0051] The additional machine learning filtering step 208 can use input features 206 extracted from the output of the germline variant caller 204, such as GATK HaplotypeCaller. These machine learning features 206 can include: (1) consensus read alignment features; (2) variant call features; (3) reference context surrounding the variant calls features; and (4) public dataset features, such as features derived from the gnomad small variant callset.

[0052] Read alignment features can be extracted by processing the read pileup at each position across the genome, or at one or more user-specified genomic regions, and for each site or position where at least one read diverges from the reference genome sequence, the extracted feature is the counts of reads which support the reference allele compared to each potential alternate allele. This method of feature extraction is highly flexible and supports excluding reads from these counts if they fail to meet thresholds of mapping quality, base quality at the divergent site, distance from variant site to the closer end of the read, and family size (where family size can be determined by clustering the reads into families based on the read’s unique molecular identifier (UMI), start position, and end position). In addition,PATENT Client Reference No.: P39561-WO-1 the method can filter out reads where the supported allele is ambiguous (e g., if there is a zero-quality base in very' close proximity, i.e. adjacent or separated by a single homopolymer, to the divergent site, or if there is a low quality' base at the site, such as a quality 18 base or worse).

[0053] The full list of features at each variant site is as follows:(1) Features related to the distribution of mapping qualities among reads which support the reference allele at the variant site include mean quality, minimum quality, maximum quality', sum of qualities, proportion of qualities below 20, proportion of qualities below 30, proportion of qualities below 40, and proportion of qualities below 60;(2) The same features as in (1), but for reads which support the ALT allele rather than the reference allele;(3) The type of variant being considered (SNV, Multi-allelic SNV, Multi-allelic SNV + InDei);(4) Counts of alt-supporting reads and ref-supporting reads;(5) Features related to how far the variant site is from the end of the read, for reads that support the reference allele at a given variant site, including mean, minimum, maximum, and sum of distances;(6) The same features as in (5), but for alt-supporting reads rather than ref-supporting reads;(7) The number of other variant calls within 100 base pair of the variant site from the initial variant calls (i.e., from the HaplotypeCaller VCF);(8) The allele frequency of the ALT allele being considered within a small variant public dataset (i.e., the gnomad small variant public dataset);(9) The two bases immediately before (when reading from the 5 ’ to 3 ’ end) the variant site in the reference genome;PATENTClient Reference No.: P39561-WO-1(10) The two bases immediately after (when reading from the 5’ to 3’ end) the variant site in the reference genome;(11) The genotype reported by the variant caller (i.e. in the GATK HaplotypeCaller VCF); and(12) The allele depths reported by the variant caller (i.e. in the GATK HaplotypeCaller VCF).

[0054] Using these features, machine learning models, such as a gradient boosting machine model such as LightGBM, are trained using a different subset of features for SNVs and multi-allelic variants. For SNVs, all of the features described above are used. For multi- allelic variants, information for both ALT alleles present in the HaplotypeCaller VCF file is recorded, but due to the typically smaller number of training examples for multi-allelics, a smaller number of features is used. The full list of features used for multi-allelic is:(1) REF allele;(2) ALT allele 1;(3) ALT allele2;(4) # of reads supporting REF allele;(5) # of reads supporting ALT1 ;(6) # of reads supporting ALT2;(7) REF allele depth as recorded by HaplotypeCaller;(8) ALT1 allele depth as recorded by HaplotypCaller; and(9) ALT2 allele depth as recorded by HaplotypeCaller.

[0055] For features 4-6, the same quality filtering as for the SNV features is used.

[0056] Other machine learning models can be used, such as random forest models, XGBoost, etc. Once these features 206 are aggregated for each variant type, the model computes for each variant the likelihood of the variant having different genotypes, includingPATENT Client Reference No.: P39561-WO-1 a “0 / 0” genotype indicating that the call was a false positive from HaplotypeCaller (i.e., since both alleles match the reference sequence, then there should be no variant called at this location). Then, the variant calls are processed using these likelihoods, and an updated VCF file 210 is produced with genotypes changed to reflect the maximum-likelihood genotypes predicted by the machine learning model, with false positives being filtered out (marked as FILTER=FAIL).

[0057] In an embodiment, the following non-default parameters are used for LightGBM model training:(1) 'objective': 'multiclass';(2) 'boosting_type': 'gbdf;(3) 'metric': 'multi_logloss';(4) 'leaming_rate': 0.01;(5) 'nrounds': 2000;(6) 'bagging_fraction': 0.9;(7) 'feature_fraction': 0.5;(8) 'bagging freq': 1;(9) 'num_leaves': 64; and(10) 'min_data_in_leaf : 20.

[0058] To support a variety of use cases, the machine learning method supports the following methods of splitting training data and testing data (i.e. sample data that the trained model can analyze to identify variants):(1) Training the model on data from a subset of chromosomes in a given dataset (e.g., chromosomes 1, 2, and 3), and using the model to produce variant calls on the remaining chromosomes;(2) Training the model on one or more full genome-wide datasets, and using the model for calling variants in a different genome-wide dataset;PATENTClient Reference No.: P39561-WO-I(3) Training the model using all available data and producing a trained model for future use; and(4) Using a pre-trained model on all available datasets to produce genome-wide variant calls.

[0059] Although the example provided above and further described herein uses GATK HaplotypeCaller as the starting short-read variant caller, other types of germline variant callers can be used instead as the starting point.Germline InDei Variant Calling

[0060] FIG. 3 illustrates an embodiment of a secondary analysis bioinformatics pipeline for calling germline InDei variants from sequencing data, and particularly from nanoporebased sequencing data. As above, this method adapts a germline variant caller, such as the GATK Haplot peCaller which has been optimized for use with Illumina sequencing data, for use with other types of sequencing data (e.g., nanopore sequencing data). The variant calls made by the GATK HaplotypeCaller can be filtered using a machine learning algorithm to remove false positives. The workflow for the germline Indel variant caller includes seven main steps as illustrated in FIG. 3:(1) Generating the consensus calls from raw read data 300;(2) Running an existing germline variant caller like GATK HaplotypeCaller on consensus reads to generate assembled haplotypes and realigned reads in active regions and generate an initial set of InDei calls 302;(3) Adding truthset information for training data (correct genotype or FP for false positives) by running a vcf evaluation tool 304 like rtg or hap.py in advance;(4) Extract a set of features 306 from realigned reads and / or from the vcf file generated by the variant caller;(5) Preprocessing of all called variants, computing several features 308 for each position that has a called variant;Y1PATENTClient Reference No.: P39561-WO-I(6) Training the machine learning model 310, such as gradient boosting model like LightGBM, based on training data; and(7) Applying the trained model on a new set of variant call data from the existing variant caller to filter out false positives and regenotype variants that were incorrectly genotyped and create a new set of InDei variant calls 312.

[0061] The workflow starts with generating consensus reads 300 from, for example, raw sequencing reads generated by a nanopore sequencer 110 or another type of sequencer. A typical format for receiving these reads is in a file using the FASTQ format. The consensus reads stored in the FASTQ file can be accessed and processed by the pipeline to perform an alignment of the consensus reads against a reference sequence. The aligned sequences can be stored in a SAM or BAM file that can be sent to a germline variant caller 204, such as GATK HaplotypeCaller, for example.

[0062] The germline variant caller 204 performs a variety of steps in order to call germline variants. First, a callable region can be identified from the aligned sequences as a region with a sufficient depth of read coverage to allow for variant calling. For example, the callable region can have a minimum depth of read coverage of at least 5, 10, 15, 20, 25, 30, 35, 40, 45, or 50. In some embodiments, the required depth of coverage can vary depending on a variety of factors, such as the types of bases in the region (i.e.. GC rich regions or other ty pe of motif which may have a higher error rate during sequencing would benefit from higher depth of coverage), and the quality scores of the bases in the region.

[0063] Next, within the callable regions, mismatches or gaps between the base calls of the reads and the reference sequence are used to identify active regions that contain potential SNVs (single nucleotide variants) and InDeis (insertions and deletions).

[0064] For each active region, an assembly graph, such as a Debruijn graph for example, is generated by starting with the reference sequence, which is decomposed into a series of k- mers (short sequences k bases long), with each successive k-mer overlapping the previous k- mer by k-1 bases. The k-mers can be represented as nodes that can be joined by lines called edges. The edges can be weighted to keep track of the number of k-mers found in the sample, with the weights initially set to zero. This forms a reference graph.PATENTClient Reference No.: P39561-WO-1

[0065] Next, each sequence read can similarly be decomposed into a series of k-mers, and can be matched to the reference graph. Each time two successive k-mers are matched to the graph, the weight of the edge joining the two k-mer nodes is incremented by one. If a k- mer cannot be matched to the graph, a new node and edge is added. This is repeated for all the k-mers in the sequence read. Under this method, the weights indicate the number of times a particular k-mer was found in the sequence reads. In addition, the weights can also be used to determine the most likely path through the graph.

[0066] In some embodiments, the assembly graph can then be optionally pruned by removing sections of the graph that are supported by an edge weight that is fewer than a threshold value, such as 2. A weight of 2, for example, means that 2 reads in the sample support that segment of the graph. The threshold can be increased, for more aggressive pruning, which will result in faster processing and higher specificity, but with lower sensitivity. Decreasing the threshold, means less pruning, which lowers the specificity, but increases the sensitivity. In some embodiments, the threshold can be set to zero, which means that the graph is not pruned for maximum sensitivity.

[0067] After the graph has been assembled and optionally pruned, haplotype sequences can be generated from the graph by traversing all paths in the graph, with a likelihood score calculated for each haploty pe sequence as the product of the transition probabilities of the path edges. The probability of an edge can be calculated as the weight of the edge divided by the sum of the weights of all the edges that share the same source node. The haplotypes with the highest likelihood scores can be used for candidate variant detection.

[0068] These haplot pes are then aligned to the reference sequence in order to identify candidate variants. The alignment can be done using a Smith-Waterman alignment (SWA), for example, and can be stored in a SAM or BAM file. The candidate variants can be determined by comparing the aligned sequence to the reference sequence and can be stored in a VCF (variant call format) file.

[0069] The steps described above can be performed, for example, using optimized settings for GATK HaplotypeCaller. an open source germline variant caller 204. The command line parameters for optimized settings is provided as an example: gatk haplotypecaller -mapping-quality-threshold-for-genotyping 1 -RF MappingQualityReadFilter -minimum-mapping-quality 1 —enable-dynamic-read-PATENTClient Reference No.: P39561-WO-1 disqualification-for-genotyping -min-base-quality-score 0 or 1 -adaptive-pruning —smithwaterman FASTEST AVAILABLE.

[0070] Using adaptive pruning and the fastest available version of Smith-Waterman alignment (SWA) reduces running time while minimum mapping and base qualities are optimized based on the chemistry (i.e., duplex / non-duplex, raw read quality, consensus base quality) and consensus method used.

[0071] Once a set of variant calls and realigned reads are provided by the germline variant caller 204 (e.g. by GATK Haplotype caller VCF and bamout), preprocessing is performed to extract features for each variant of interest 306. The collection of features includes: number of reads that support a variant observed in the BAM and number of reads that support the reference, along with other attributes such as mapping and base quality, distance of variant from the end of the read and whether the read is duplex or not.

[0072] For extracting the features, several considerations are taken into account during the process to reduce errors in InDeis, which often predominantly occur in regions with tandem repeats and homopolymers. To compute the features accurately within these regions, the set of reads supporting each allele and reference is refined, ty pically by the exclusion of ambiguous read supports. For example as shown in FIG. 4, reads not extending across the repeat region are ambiguous and discarded in counting effective REF / ALT support 400. Reads with low base quality (e.g. 0 or 18 for duplex reads or 0 to 13 for other types of reads and a deep learning consensus method or 0 to 38 for other types of reads and majority vote consensus method) at the site of the variant 402, and alternate allele (ALT) supporting reads with zero base quality in the inserted sequence are also disregarded or discarded 404.

[0073] Due to a typically higher InDei error rate in homopolymer and tandem repeats, a higher than normal rate of discordant base (with low Q score) can be observed. Therefore, discarding reads with low Q scores can lead to a significant drop in coverage. However, these reads can still be useful if leveraged appropriately. These reads can be utilized by incorporating raw read sequence information prior to intra-molecular consensus formation, referred to as Raw7read 1 (Rl) and Raw7read 2 (R2) in FIG. 5. During consensus formation, the precise sequence of the raw reads used to form the consensus sequence (e.g. the length of homopolymer supported by Rl and R2) can be stored and later restored by utilizing a lossless encoding technique. Lossless encoding can preserve raw read sequence information, using forPATENTClient Reference No.: P39561-WO-1 instance a Compact Idiosyncratic Gapped Alignment Report (CIGAR) type tag separately for each read in the BAM file, thereby allowing the reconstruction of R1 and R2 during variant calling, if needed. As shown in FIG. 5, a partial ALT / REF support count can be assigned to reads with low Q scores (discordant R1 and R2). Consequently, bases with a low quality score (e.g., 0) are not ignored and they can partially contribute to the support of REF or ALT alleles. The InDei model is capable of optionally incorporating these partial weights in the model as new' features and the increased support from partial weights can enhance the accuracy of the model, especially in homopolymer-rich regions.

[0074] The lossless encoding tags can be generated by comparing the read sequence to the reference sequence and then identifying differences between the two. For example, Table 1 below shows an example chart that can be used to encode a read that contains A, T, C, and G nucleotides using 14 possible encodings. The encodings shown in Table 1 are just an example, and can be modified. The sequence read may then be encoded into a text or a bit string using the encodings to create the CIGAR type tag mentioned above. The bit string or text that is encoded at the base level can then be compressed in later steps. The encodings include a match, the 4 substitutions, 4 soft clips (the end of a read is not aligned). 4 insertions, and a deletion.Table 1PATENTClient Reference No.: P39561-WO-I

[0075] The genomic location information in the reference sequence is substituted for at least a portion of the sequence that matches the reference sequence. For example, if a portion of the nucleotides in the beginning of a sequence matches with the reference sequence and then there is one or more mismatches, the nucleotides in the first portion can be replaced by a start location relative to the reference sequence, a number that shows the length of the portion, and the code that represents a mismatch. The one or more mismatches may then remain as encoded. Any portion of matching sequences may similarly be replaced (i.e., to compress the sequence data) by a start location corresponding to the position of a first matching nucleotide and a length of the portion of matching sequences. The code for a sequence match may or may not be included. A portion of the sequence that matches with a reference sequence may be 2 bases, 3 bases, 5 bases, 10 bases, 20 bases, 30 bases, 40 bases, 100 bases, 500 bases, or longer. The portion can then be substituted with, for example, only 3 numbers including a chromosome number, a start location for a location of the first nucleotide in the portion that matches with the reference sequence, and the length of the portion. In some embodiments, the length of the read must be stored as part of the location and identification of the matching bases, and may be used to decode the final compressed data.

[0076] Although the use of lossless encoding is described above in connection with the InDei variant calling model, it can also be optionally used in a similar manner for the SNV variant caller model as described above (see FIG. 2) to recover additional sequencingPATENTClient Reference No.: P39561-WO-1 information from the raw reads. The SNV model is capable of optionally incorporating these partial weights in the model as new features and the increased support from partial weights can enhance the accuracy of the model.

[0077] In the preprocessing step 308, the previously extracted features are processed and converted to numerical features. For this step, multiple methods for converting these features into numeric form can be used. For example, from the set of reads supporting each variant, the following features are extracted:(1) Number of reads supporting reference and allele (after removing read supports near a base with 0 base quality in repeats);(2) Sum of all mapping qualities and base qualities, and distance of variant from the end of the reads. For REF / ALT larger than one base, the average of the value of the feature for bases in the ALT or REF is computed for each read;(3) Maximum and minimum of the mapping qualities, and distance of variant from the end of the reads;(4) Proportion of reads with mapping quality’ below 20, 30, 40 and 60 for ALT and REF supporting reads;(5) Similar to 3 and 4, if base qualities have continuous values, there are minimum base quality’ of the list of base qualities, and proportions of reads with base qualities below certain thresholds (refers to bases of the reads that cover a position, i.e., if position 1 has a depth of 10 there will be 10 bases each with a quality’ score for that position. If there are multiple bases in REF or ALT, we compute the average of the value of the feature for the multiple bases for each 10 reads);(6) If unique molecular identifier (UMI) data is present, minimum, maximum, sum and proportions of reads with family size below7certain thresholds (e.g. 5 and 10, which can be dependent on chemistry and actual data);(7) Weighted Alt and REF support, which is the sum of mapping quality multiplied by base quality' for all reads (if continuous values of base qualities exists), normalized by maximum mapping quality and base quality;PATENT Client Reference No.: P39561-WO-1(8) Allele frequency based on the feature extraction output, also weighted by base quality and mapping quality (i.e., sum of mapping quality of ALT supporting reads divided by the sum of mapping quality of all reads); and(9) Number of reads partially supporting reference and ALT allele with low base quality bases in them (from lossless encoding)(note that this feature can also be extracted for use in the SNV model described in FIG. 2); and(10) Total depth covering the position.

[0078] Most of these features are counted for simplex reads as well as duplex read supports. In addition to the features above, additional features can be extracted from the VCF file provided by the germline variant caller 204:(11) Ref and ALT sequence, AD (allelic depths for the ref and alt alleles), QD (variant confidence / Quality by Depth), HAPCOMP (edit distances of each alt allele's most common supporting haplotype from closest germline haplotype, excluding differences at the site in question), HAPDOM (for each alt allele, fraction of read support that best fits the most-supported haplotype containing the allele); and(12) ADT and ADTL computed based on AD field of VCF and also from the corrected AD values from feature extraction above, where ADT and ADTL evaluate imbalances in calling the reference and the alternative allele (ADT), also depending on the amount of reads that map on the variant locus (ADTL). They are defined as ADT=|(AD1-AD2)| / (AD1+AD2), ADTL=logl O(AD1+AD2)*ADT where ADI and AD2 are the number of unfiltered reads calling the reference and alternate allele, respectively.

[0079] In addition, a set of features can be extracted from the reference about the region in which variant resides:(13) 2 base pair context before and after the variant;(14) Whether the variant is from a homopolymer or tandem repeat, and the repeat unit and the number of times repeat unit is observed in reference;(15) Complexity of reference sequence at that region by number of unique K-mers for K=3,4,5,6;PATENTClient Reference No.: P39561-WO-I(16) Number of variations observed in each position compared to reference; and(17) Length of InDei.

[0080] In addition, an additional feature is the allele frequency (AF) of each allele from a public dataset, such as a genomic database like the gnomAD database.

[0081] Once all features are computed for each allele, to make a prediction for a position that a call has been made, and a ML model is trained on the features. Using this setup, a machine learning model, such as a gradient boosting machine (e.g. LightGBM, etc.) or other models (such as, but not limited to, XGBoost, random forest, etc.) can be used. The parameter set used for training the model is optimized as follows for the current data (duplex chemistry and the first three chromosomes as training set):(1) 'objective': 'multiclass',(2) 'boosting_type': 'gbdf,(3) 'metric': 'multi logloss',(4) 'num_classes': 4,(5) 'numjeaves': 59,(6) 'lambda': 0,(7) 'bagging fraction': 0.880623.(8) 'num_iterations' :8000,(9) 'learning rate' :0.0420473,(10) feature_fraction: 0.902813.

[0082] These parameters are optimized based on models trained using six input BAM files.

[0083] Generally, training data is an independent dataset from an evaluation (test) set, e.g. other samples with different genetic makeup (unique genotype). However, in the case of lack of data availability, users can train the model on a subset of chromosomes that is thePATENTClient Reference No.: P39561-WO-1 same as the experimental setup. The size of the training data set can also affect the optimized parameter values.

[0084] Once the model is trained, the model can be used to perform germline InDei calling on test samples with the following steps which are mostly similar to training process except the last two steps:(1) Generating the consensus calls from raw data;(2) Running an existing germline variant caller like GATK HaplotypeCaller on consensus reads to generate assembled haplotypes and realigned reads in active region and initial set of variant calls;(3) Extract a set of features from realigned reads;(4) Preprocessing of all called variants, computing several features for each position that has a called variant;(5) Predict calls using the trained machine learning model (e.g.. LightGBM);(6) Post-process the vanants to merge multi-allelic calls with appropriate genotype (1 / 2); and(7) Create a standard VCF from the model predictions.

[0085] After completing all the necessary preparation steps for the test / sample data, the ML model is deployed 312 to generate predictions for each position with observed variants. These predictions categorize each call as either a false positive, indicating that the call should be filtered out (with 0 / 0 genotype), or as one of the possible genotypes (1 / 1, 0 / 1, and 1 / 2). Following the generation of variant predictions for the test / sample data, a standard VCF file is produced. This file encompasses correct alleles and genotypes, plus essential information that influenced the Machine Learning model's predictions like sum of the mapping qualities for REF and ALT, length of homopolymer if the variant resides in a homopolymer, and average distance of variant to the ends of the read. Consequently, the output of this module is compatible with existing VCF evaluation tools. This compatibility allows for an accurate reporting of the accuracy of the final set of calls, satisfying the need for robust evaluation in variant calling tasks. Our benchmarks with datasets with known truth sets like genome in a bottle project, showed very high sensitivity and specificity for the final calls made by this toolPATENTClient Reference No.: P39561-WO-1 and comparable accuracy with Illumina data, which was not possible with any of the existing tools.Computer-Implemented Methods for SNV Calling

[0086] FIGS. 6A & 6B are flow charts of a method 600 for performing SNV calling, in accordance with at least some embodiments. The method 600 can be performed, at least in part, by at least one processor. In some embodiments, the at least one processor comprises a general programmable processor configured to execute a set of instructions stored on a computer-readable medium. In other embodiments, the at least one processor may include any combination of hardware logic circuits, such as implemented on an application specific integrated circuit (ASIC) or field programmable gate array (FPGA), or general programmable processors executing any combination of firmware and / or software. It will be appreciated that any system capable of implementing one or more steps of the method 600 is within the scope of the present disclosure.

[0087] At 602. sequencing data is received from a sequencer instrument. In an embodiment, the sequencer instrument is a nanopore-based sequencer instrument configured to perform a sequencing assay based on a sequencing-by-expansion (SBX) chemistry. The sequencing data can be provided in a FASTQ file.

[0088] At 604, a plurality of consensus reads are generated from the sequencing data, and, at 606, the plurality of consensus reads are aligned to a reference sequence. In an embodiment, the consensus reads are aligned to a reference sequence and then the aligned reads are used to identify one or more callable regions in the genome.

[0089] At 608, a first set of germline variants are called using a germline variant caller. For each callable region, haploty pe information is obtained using a graph-based approach. The haplotype information is then used to identify a number of candidate variants. The candidate variants may be filtered using various criteria to generate a final set of variant calls, which may be provided in a VCF file. In an embodiment, the germline variant caller is the GATK HaplotypeCaller.

[0090] At 610, a set of features is extracted from the variant information generated by the germline variant caller. The features can include read alignment features, variant call information features, reference sequence context features, and / or public dataset features.PATENTClient Reference No.: P39561-WO-1

[0091] At 612, the extracted set of features for a number of called variants are used to train a machine learning model to predict corrected genotype information for the called variants. In an embodiment, each variant may be identified as either: 0 / 0, 0 / 1, 1 / 1, 1 / 2, and so forth, with 0 / 0 indicating a false positive in the called variant information. The ML model may be, e.g., a GBM such as LightGBM.

[0092] The method 600 can continue with the steps of method 650, which can be performed once the ML model is trained using the steps of method 600 for a number of different samples in a training dataset.

[0093] At 652, the trained ML model is applied to a second set of variant information generated by the germline variant caller to predict updated genotype information for the observed variants in the second set of variant information. The second set of variant information can be selected from a subset of variant information in the same sample used to train the model (e.g., variants from different chromosomes of the same sample), or from different samples from different patients, for example.

[0094] At 654, a set of filtered germline variant calls are generated from the second set of germline variants based on the updated genotype information predicted by the ML model. The set of filtered germline variant calls can exclude (i. e.. remove) those variants in the second set of variant information that are categorized as false positives (i.e., predicted genotype is 0 / 0).

[0095] At 656, an updated genotype is generated based on the set of filtered germline variant calls. In an embodiment, a new VCF file can be generated that includes the set of filtered germline variant calls, omitting false positives, and can include information about the corrected genotype data predicted by the ML model.Exemplary Computer System

[0096] FIG. 7 illustrates an example computer system 700 that may be utilized to implement techniques disclosed herein. Any of the computer systems mentioned herein, such as for hosting the systems and implementing the processes described for calculating variants, may utilize any suitable number of subsystems. Examples of such subsystems are shown in FIG. 7 as included in computer system 700. In some embodiments, a computer system includes a single computer apparatus, where the subsystems can be the components of thePATENTClient Reference No.: P39561-WO-I computer apparatus. In other embodiments, a computer system can include multiple computer apparatuses, each being a subsystem, with internal components. A computer system can include desktop and laptop computers, tablets, mobile phones, telecommunication devices or other mobile devices. In some embodiments, a cloud infrastructure (e.g., Amazon Web Services), a graphical processing unit (GPU), etc., can be used to implement the disclosed techniques.

[0097] The subsystems shown in FIG. 7 are interconnected via a system bus 75. Additional subsystems such as a printer 74, keyboard 78, storage device(s) 79, monitor 76, which is coupled to display adapter 82, and others are shown. Peripherals and input / output (I / O) devices, which couple to I / O controller 71, can be connected to the computer system by any number of means known in the art such as input / output (I / O) port 77 (e.g., USB, FireWire®). For example, I / O port 77 or external interface 81 (e.g. Ethernet, Wi-Fi, etc.) can be used to connect computer system 700 to a wide area network such as the Internet, a mouse input device, or a scanner.

[0098] The interconnection via system bus 75 allows the central processor 73 to communicate with each subsystem and to control the execution of a plurality of instructions from system memory 72 or the storage device(s) 79 (e.g., a fixed disk, such as a hard drive, or optical disk), as well as the exchange of information between subsystems. The system memory 72 and / or the storage device(s) 79 may embody a computer readable medium. Another subsystem is a data collection device 85, such as a camera, microphone, accelerometer, and the like. Any of the data mentioned herein can be output from one component to another component and can be output to the user.

[0099] A sequencer instrument 90 (e.g., a nanopore sequencer 110), which may sometimes be referred to as a sequencing device, is connected through external interface 81 for providing sequencing data to a data collection device 85 and / or storage devices 79.

[0100] A computer system 700 can include a plurality of the same components or subsystems, e.g., connected together by external interface 81 or by an internal interface such as the bus 75. In some embodiments, computer systems, subsystems, or apparatuses can communicate over a network. In such instances, one computer can be considered a client and another computer a server, where each can be part of the same computer system. A client and a server can each include multiple systems, subsystems, or components.PATENTClient Reference No.: P39561-WO-1

[0101] Aspects of embodiments can be implemented in the form of control logic using hardware (e.g. an application specific integrated circuit or field programmable gate array) and / or using computer software with a generally programmable processor in a modular or integrated manner. As used herein, a processor includes a single-core processor, multi-core processor on a same integrated chip, or multiple processing units on a single circuit board or networked. Based on the disclosure and teachings provided herein, a person of ordinary skill in the art will know and appreciate other ways and / or methods to implement embodiments of the present invention using hardware and a combination of hardware and software.

[0102] Machine learning models utilized herein may include one or more of a Naive Bayes (NB) model, a logistic regression (LR) model, a random forest (RF) model, a support vector machine (SVM) model, an artificial neural network model, a multilayer perceptron (MLP) model, a convolutional neural network (CNN), a Large Language model (LLM), and / or other machine learning or deep learning models, etc. The machine learning models can be updated / trained using a supervised learning technique, an unsupervised learning technique, etc.

[0103] Any of the software components or functions described in this application may be implemented as software code to be executed by a processor using any suitable computer language such as, for example, Java. C, C++, C#, Objective-C, Swift, or scripting language such as Perl or Python using, for example, conventional or object-oriented techniques. The software code may be stored as a series of instructions or commands on a computer readable medium for storage and / or transmission. A suitable non-transitory computer readable medium can include random access memory (RAM), a read only memory (ROM), a magnetic medium such as a hard-drive or a floppy disk, or an optical medium such as a compact disk (CD) or DVD (digital versatile disk), flash memory, and the like. The computer readable medium may be any combination of such storage or transmission devices.

[0104] Such programs may also be encoded and transmitted using carrier signals adapted for transmission via wired, optical, and / or wireless networks conforming to a variety of protocols, including the Internet. As such, a computer readable medium may be created using a data signal encoded with such programs. Computer readable media encoded with the program code may be packaged with a compatible device or provided separately from other devices (e.g.. via Internet download). Any such computer readable medium may reside on or within a single computer product (e.g. a hard drive, a CD, or an entire computer system), andPATENTClient Reference No.: P39561-WO-1 may be present on or within different computer products within a system or network. A computer system may include a monitor, printer, or other suitable display for providing any of the results mentioned herein to a user.

[0105] Any of the methods described herein may be totally or partially performed with a computer system including one or more processors, which can be configured to perform the steps. Thus, embodiments can be directed to computer systems configured to perform the steps of any of the methods described herein, potentially with different components performing a respective step or a respective group of steps. Although presented as numbered steps, steps of methods herein can be performed at the same time or in a different order. Additionally, portions of these steps may be used with portions of other steps from other methods. Also, all or portions of a step may be optional. Additionally, any of the steps of any of the methods can be performed with modules, units, circuits, or other means for performing these steps.

[0106] When a feature or element is herein referred to as being “on” another feature or element, it can be directly on the other feature or element or intervening features and / or elements may also be present. In contrast, when a feature or element is referred to as being “directly on” another feature or element, there are no intervening features or elements present. It will also be understood that, when a feature or element is referred to as being “connected”, “attached” or “coupled” to another feature or element, it can be directly connected, attached or coupled to the other feature or element or intervening features or elements may be present. In contrast, when a feature or element is referred to as being “directly connected”, “directly attached” or “directly coupled” to another feature or element, there are no intervening features or elements present. Although described or shown with respect to one embodiment, the features and elements so descnbed or shown can apply to other embodiments. It will also be appreciated by those of skill in the art that references to a structure or feature that is disposed “adjacent” another feature may have portions that overlap or underlie the adjacent feature.

[0107] Terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. For example, as used herein, the singular forms “a”, “an” and “the” are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms “comprises” and / or “comprising,” when used in this specification, specify the presence of stated features,PATENTClient Reference No.: P39561-WO-1 steps, operations, elements, and / or components, but do not preclude the presence or addition of one or more other features, steps, operations, elements, components, and / or groups thereof. As used herein, the term “and / or” includes any and all combinations of one or more of the associated listed items and may be abbreviated as “ / ”.

[0108] Spatially relative terms, such as “under”, “below”, “lower”, “over”, “upper” and the like, may be used herein for ease of description to describe one element or feature’s relationship to another element(s) or feature(s) as illustrated in the figures. It will be understood that the spatially relative terms are intended to encompass different orientations of the device in use or operation in addition to the orientation depicted in the figures. For example, if a device in the figures is inverted, elements described as “under” or “beneath” other elements or features would then be oriented “over” the other elements or features. Thus, the exemplary term “under” can encompass both an orientation of over and under. The device may be otherwise oriented (rotated 90 degrees or at other orientations) and the spatially relative descriptors used herein interpreted accordingly. Similarly, the terms “upwardly”, “downwardly”, “vertical”, “horizontal” and the like are used herein for the purpose of explanation only unless specifically indicated otherwise.

[0109] Although the terms “first” and “second” may be used herein to describe various features / elements (including steps), these features / elements should not be limited by these terms, unless the context indicates otherwise. These terms may be used to distinguish one feature / element from another feature / element. Thus, a first feature / element discussed below could be termed a second feature / element, and similarly, a second feature / element discussed below could be termed a first feature / element without departing from the teachings of the present invention.

[0110] Throughout this specification and the claims which follow, unless the context requires otherwise, the w ord “comprise”, and variations such as “comprises” and “comprising” means various components can be co-jointly employed in the methods and articles (e.g.. compositions and apparatuses including device and methods). For example, the term “comprising” will be understood to imply the inclusion of any stated elements or steps but not the exclusion of any other elements or steps.

[0111] As used herein in the specification and claims, including as used in the examples and unless otherwise expressly specified, all numbers may be read as if prefaced by the wordPATENTClient Reference No.: P39561-WO-1“about” or “approximately,” even if the term does not expressly appear. The phrase “about” or “approximately” may be used when describing magnitude and / or position to indicate that the value and / or position described is within a reasonable expected range of values and / or positions. For example, a numeric value may have a value that is + / - 0. 1% of the stated value (or range of values), + / - 1% of the stated value (or range of values), + / - 2% of the stated value (or range of values), + / - 5% of the stated value (or range of values), + / - 10% of the stated value (or range of values), etc . Any numerical values given herein should also be understood to include about or approximately that value, unless the context indicates otherwise. For example, if the value “10” is disclosed, then “about 10” is also disclosed. Any numerical range recited herein is intended to include all sub-ranges subsumed therein. It is also understood that when a value is disclosed that “less than or equal to” the value, “greater than or equal to the value” and possible ranges between values are also disclosed, as appropriately understood by the skilled artisan. For example, if the value “X” is disclosed the “less than or equal to X” as well as “greater than or equal to X” (e.g., where X is a numerical value) is also disclosed. It is also understood that the throughout the application, data is provided in a number of different formats, and that this data, represents endpoints and starting points, and ranges for any combination of the data points. For example, if a particular data point “10” and a particular data point “15” are disclosed, it is understood that greater than, greater than or equal to, less than, less than or equal to, and equal to 10 and 15 are considered disclosed as well as between 10 and 15. It is also understood that each unit between two particular units are also disclosed. For example, if 10 and 15 are disclosed, then 11, 12. 13. and 14 are also disclosed.

[0112] Although various illustrative embodiments are described above, any of a number of changes may be made to various embodiments without departing from the scope of the invention as described by the claims. For example, the order in which various described method steps are performed may often be changed in alternative embodiments, and in other alternative embodiments one or more method steps may be skipped altogether. Optional features of various device and system embodiments may be included in some embodiments and not in others. Therefore, the foregoing description is provided primarily for exemplary purposes and should not be interpreted to limit the scope of the invention as it is set forth in the claims.PATENTClient Reference No.: P39561-WO-1

[0113] The examples and illustrations included herein show, by way of illustration and not of limitation, specific embodiments in which the subject matter may be practiced. As mentioned, other embodiments may be utilized and derived there from, such that structural and logical substitutions and changes may be made without departing from the scope of this disclosure. Such embodiments of the inventive subject matter may be referred to herein individually or collectively by the term “invention” merely for convenience and without intending to voluntarily limit the scope of this application to any single invention or inventive concept, if more than one is, in fact, disclosed. Thus, although specific embodiments have been illustrated and described herein, any arrangement calculated to achieve the same purpose may be substituted for the specific embodiments shown. This disclosure is intended to cover any and all adaptations or variations of various embodiments. Combinations of the above embodiments, and other embodiments not specifically described herein, will be apparent to those of skill in the art upon reviewing the above description.

Claims

PATENTClient Reference No.: P39561-WO-1CLAIMS1. A method for germline variant calling, the method comprising: receiving sequencing data generated by a sequencer instrument performing a sequencing assay; generating a plurality of consensus reads from the sequencing data; aligning the plurality of consensus reads to a reference sequence to generate alignment data; identifying at least one callable region of the sequencing data; identify ing, for each callable region in the at least one callable region using a germline variant caller, variant information for a plurality of observed variants in the callable region based on the alignment data; extracting a plurality of features from the vanant information; and training a machine learning (ML) model to predict a corrected genotype for at least one variant of the plurality of observed variants based on the extracted plurality of features.

2. The method of claim 1, wherein the germline variant caller comprises a GATK Haplotype Caller.

3. The method of claim 2, wherein optimized settings for the GATK Haplotype Caller comprise at least one of: mapping-quality -threshold-for-genotyping 1 ; minimum-mapping-quality' 1; enable-dynamic-read-disqualification-for-genotyping; min-base-quality-score 0; adaptive-pruning true; and smith-waterman FASTEST AV AIL ABLE.

4. The method of claim 1 , wherein the plurality of features includes at least one feature of each of the following: one or more read alignment features from the alignment data; one or more variant call features from the variant information; one or more reference context features based on the reference sequence; andPATENTClient Reference No.: P39561-WO-1 one or more features extracted from a public dataset.

5. The method of claim 4, wherein the public dataset comprises the gnomAD public database.

6. The method of claim 1, wherein the plurality of features for an observed variant at a specific variant site comprise at least one of: distribution of quality scores among reads that support a reference allele at the variant site; distribution of quality scores among reads that support an alternative allele at the variant site; a type of the observed variant; a count of reads that support at least one of the reference allele or the alternative allele; a distance of the variant site from the end of each read, or a mean, minimum, maximum, or sum of distances of the variant site from a plurality of reads that support a particular allele; a number of other variant calls within a given distance of the variant site; an allele frequency of the alternative allele based on the public dataset; a number of bases immediately preceding or immediately following the variant site; and / or a depth of coverage for the alternative allele.

7. The method of claim 1, wherein the ML model is a gradient boosting machine (GBM) model.PATENTClient Reference No.: P39561-WO-18. The method of claim 7, wherein the ML model is lightGBM, and wherein the ML model is configured with the following parameters:(1) 'objective': 'multiclass';(2) 'boosting_type': 'gbdt';(3) 'metric': 'multi logloss';(4) 'leaming_rate': 0.01 ;(5) 'mounds': 2000;(6) 'bagging_fraction': 0.9;(7) eature_fraction': 0.5;(8) 'bagging_freq': 1;(9) 'num_leaves': 64; and(10) 'min_data_in_leaf : 20.

9. The method of claim 1 , the method further comprising: applying the trained machine learning model to a second set of variant information generated by the germline variant caller to re-genotype observed variants in the second set of variant information that were incorrectly genotyped by the germline variant caller; generating a set of filtered germline / regenotyped variant calls from the second set of variant information; and generating an updated genotype for the sample based on the set of filtered germline / regenotyped variant calls.

10. The method of claim 9, wherein the set of filtered germline / regentotyped variant calls comprises at least one single nucleotide variant (SNV).

11. The method of claim 9, wherein the observed variants comprise at least one InDei.

12. The method of claim 1 , wherein the sequencing data comprises sequencing-by- expansion (SBX) sequencing data measured by a nanopOre-based sequencer instrument.

13. The method of claim 12, wherein the germline variant caller is optimized for non- SBX sequencing data.PATENTClient Reference No.: P39561-WO-114. The method of claim 13, wherein the germline variant caller is optimized for Illumina-based sequencing data.

15. A system for performing germline variant calling, the system comprising: a memory storing sequencing data generated by a sequencer instrument performing a sequencing assay; and at least one processor configured to: generate a plurality of consensus reads from the sequencing data, align the plurality of consensus reads to a reference sequence to generate alignment data, identify at least one callable region of the sequencing data, identify, for each callable region in the at least one callable region using a germline variant caller, variant information for a plurality of observed variants in the callable region based on the alignment data, extract a plurality of features from the variant information, and train a machine learning (ML) model to predict a corrected genoty pe for at least one variant of the plurality of observed variants based on the extracted plurality of features.

16. The system of claim 15, further comprising the sequencer instrument communicatively coupled to the at least one processor, wherein the sequencer instrument is a nanopore-based sequencer instrument configured to perform a sequencing-by-expansion (SBX) sequencing assay.

17. The system of claim 15, wherein the plurality7of features includes at least one feature of each of the following: one or more read alignment features from the alignment data; one or more variant call features from the variant information; one or more reference context features based on the reference sequence; and one or more features extracted from a public dataset.PATENTClient Reference No.: P39561-WO-118. The system of claim 15, wherein the one or more features for an observed variant at a specific variant site comprise at least one of: distribution of quality scores among reads that support a reference allele at the variant site; distribution of quality scores among reads that support an alternative allele at the variant site; a type of the observed variant; a count of reads that support at least one of the reference allele or the alternative allele; a distance of the variant site from the end of each read, or a mean, minimum, maximum, or sum of distances of the variant site from a plurality of reads that support a particular allele; a number of other variant calls within a given distance of the variant site; an allele frequency of the alternative allele based on the public dataset; a number of bases immediately preceding or immediately following the variant site; and / or a depth of coverage for the alternative allele.

19. The system of claim 15, wherein the ML model is a gradient boosting machine (GBM) model.

20. A non-transitory computer-readable medium storing instructions that, executed by at least one processor, cause a computer device to: receive sequencing data generated by a sequencer instrument performing a sequencing assay; generate a plurality of consensus reads from the sequencing data; align the plurality of consensus reads to a reference sequence to generate alignment data;PATENTClient Reference No.: P39561-WO-1 identify at least one callable region of the sequencing data; identify, for each callable region in the at least one callable region using a germline variant caller, variant information for a plurality of observed variants in the callable region based on the alignment data; extract a plurality of features from the variant information; and train a machine learning (ML) model to predict a corrected genotype for at least one variant of the plurality of observed variants based on the extracted plurality of features.

Images