Machine learning based filtering for population level joint calling quality control

Abstract

A machine learning (ML) model may be trained and / or implemented to assist in quality control of variant calls in cohort level sequencing data. A computing device may receive a variant call file comprising cohort level sequencing data during training. Training data may be identified based on a predefined set of features related to genotyping rate or quality distribution. The training data may include positive labeled training data from at least one high-confidence region of the genome for variant calling and / or negative labeled training data from the variants identified from at least one low-confidence region of the genome for variant calling. The training data may be used to train the ML model to predict a machine learning site quality (MLSQ) score. The trained ML model may predict the MLSQ score for each variant. The MLSQ score may be used to filter the variants identified in the cohort level sequencing data.

Description

MACHINE LEARNING BASED FILTERING FOR POPULATION LEVEL JOINTCALLING QUALITY CONTROLCROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims priority to U.S. Provisional Patent Application No. 63 / 696,602, filed September 19, 2024 and entitled “MACHINE LEAKING BASED FILTERING FOR POPULATION LEVEL JOINT CALLING QUALITY CONTROL,” the entirety of which is incorporated by reference herein in its entirety.BACKGROUND

[0002] Population-scale genomics experiments can include aggregating and / or merging data associated with variants from a large number (e.g., hundreds of thousands to millions) of samples. In large population-scale analysis, quality control of genetic variants and genotypes may significantly impact the results of downstream analysis. However, the truth genotype data that is publicly available is limited (e.g., 7 total samples). This imposes challenges in assessing the quality of variants discovered in cohort level sequencing data with hundreds of thousands to millions of samples, and in comparing call sets generated from different projects, sites, and / or pipelines that are aggregated into a single dataset for common analysis.

[0003] The existing sequencing analyses may use a genome variant call format (gVCF) file. The gVCF file stores sequencing information for both variant and non-variant positions. The gVCF file may enable representation of genotype, annotation, and other information across all sites in the genome. gVCF genotypers may include population-based analysis tools that jointly analyze variants from unrelated individuals. However, there can be variation in the number of called variants, genotyping rates, and / or genotyping consistency between gVCF genotypers, particularly in low-confidence regions of the genome in which detecting variants may be a known challenge.SUMMARY

[0004] Systems, methods, and apparatus are described herein for training and / or implementing a machine learning (ML) model to assist in quality control of variant calls in cohort level sequencing data. As described herein, one or more computing devices may be configured to receive at least one variant call file comprising cohort level sequencing data during training. The cohort level sequencing data received in the at least one variant call file may indicate possible variants identified in the cohort level sequencing data. The training data may be identified from the at least one variant call file based on a predefined set of features related to genotyping rate or distribution quality metrics. For example, the predefined set of features may include at least one of a variant type, a number of samples with missing genotype, a number of samples with no coverage, a variant genotyping quality (GQ) histogram, or a reference GQ histogram. The training data may include positive labeled training data and negative labeled training data from the variants identified in the cohort level sequencing data. The positive labeled training data may include the variants with a genotyping rate above a genotyping threshold from at least one high- confidence region of the genome for variant calling. The negative labeled training data may include the variants with a genotyping rate below the genotyping threshold from at least one low- confidence region of the genome for variant calling. The training data may be used to train the ML model to predict cohort that may be used to indicate a genotyping rate or a quality distribution. The trained parameters may be stored for implementation.

[0005] During implementation, the ML model may receive at least one variant call file comprising cohort level sequencing data. The trained ML model may be implemented to predict a machine learning site quality (MLSQ) score for each variant in the at least one variant call file based on the input data from the at least one variant call file. The predicted MLSQ score may be stored for each variant in the memory. The MLSQ score may be used to filter the variants identified in the cohort level sequencing data in the at least one variant call file having the MLSQ score below a filter threshold. The filtered variants may be stored in an updated variant call file in memory. The filtered variants may be stored in the updated variant call file with an indication of being filtered based on the MLSQ score being below the filter threshold.

[0006] The ML model may be trained based on cohort level sequencing data related to monozygotic twins and / or trios. The ML model may be trained based on a genotype consistency in trios above a consistency threshold in the at least one high-confidence region in the genomeand / or below a consistency threshold in the at least one low-confidence region in the genome. The ML model may be trained based on a genotype consistency in monozygotic twins above a consistency threshold in the at least one high-confidence region in the genome and / or below a consistency threshold in the at least one low-confidence region in the genome.

[0007] The at least one variant call file may include hundreds of thousands to millions of samples. The at least one variant call file comprising the cohort level sequencing data may be partitioned into a plurality of variant call files based on regions of the genome. Each variant call file of the variant call files may be processed sequentially or in parallel based on a number of threads or instances operating the ML model. Different variant call files corresponding to specific regions of the genome may be processed at different threads, computing devices, or other processing entities.

[0008] Updates may be performed to the ML model. For example, updates to the ML model may be performed to re-train the ML model. The collected feature values from the training data may be stored in a separate file that may be utilized during an updated training procedure. The MLSQ score may be reconfigured on different features or may be re-trained on a specific feature by updating the training process for a particular feature. Updated features and / or feature values may be implemented and the ML model may be retrained on the updated features and / or feature values. For example, the feature list may be updated or may be the same and the feature values may be updated based on updated training data. In another example, the training data may be updated by updating the labels on the positive and / or negative labeled training data. A lookup may be performed on the previous training data to perform updated training on a subset of the training data.

[0009] The ML model may be trained on cohort level sequencing data (e.g., multi-sample variant call format (msVCF) files, databases, etc.) that is generated from different versions of a genotyper configured to generate the cohort level sequencing data (e.g., msVCF files, databases, etc.). For example, different versions of the genotyper may generate different forms of cohort level sequencing data, which may include different data related to genotyping including genotyping rates and / or distributions of quality. Different versions of the genotyper may generate different feature values for each of the features on which the ML model is trained. The ML model may be updated for different versions of the iterative genotyper by identifyingdifferent feature values in the cohort level sequencing data generated by different versions of the genotyper and may re-train the ML model based on the updated feature values.

[0010] Different ML models may be trained for specific cohorts using training data specific to the cohort. For example, the ML models may be customized for a specific cohort on which secondary and / or tertiary analysis may be performed. The ML model may be updated for improving the predicted MLSQ scores for a particular cohort and / or to add / remove individuals from the cohort.BRIEF DESCRIPTION OF DRAWINGS

[0011] FIG. 1 illustrates a schematic diagram of a system environment (or “environment”).

[0012] FIG. 2 illustrates a block diagram of an example computing device.

[0013] FIG. 3 depicts an example process that may be implemented by one or more computing devices to perform iterative gVCF genotyping, as described herein.

[0014] FIG. 4 is a flowchart depicting an example procedure for training an ML model to predict a machine learning site quality (MLSQ) score based on feature values from high and / or low confidence regions in cohort level sequencing data.

[0015] FIG. 5 illustrates an example system environment and / or process for training and / or implementing a supervised ML model.

[0016] FIG. 6 is a flowchart depicting an example procedure for implementing an ML model for predicting an MLSQ score and applying a genotyping filter for assisting with maintaining genotyping rate and / or genotyping consistency in variant calls in cohort level sequencing data.DETAILED DESCRIPTION

[0017] FIG. 1 illustrates a schematic diagram of a system environment (or “environment”) 100 in which an iterative genome variant call format (gVCF) genotyper may operate, as described herein. As illustrated, the environment 100 includes one or more server device(s) 102 connected to a client device 108 and a sequencing device 114 via a network 112.

[0018] As shown in FIG. 1, the server device(s) 102, the client device 108, and the sequencing device 114 may communicate with each other via the network 112. The network 112 may comprise any suitable network over which computing devices can communicate. The network 112 may include a wired and / or wireless communication network. Example wirelesscommunication networks may be comprised of one or more types of RF communication signals using one or more wireless communication protocols, such as a cellular communication protocol, a WIFI communication protocol, and / or another wireless communication protocol. In addition, or in the alternative to communicating across the network 112, the server device(s) 102, the client device 108, and / or the sequencing device 114 may bypass the network 112 and may communicate directly with one another.

[0019] As indicated by FIG. 1, the sequencing device 114 may comprise a device for sequencing a biological sample. The biological sample may include human and non-human deoxyribonucleic acid (DNA) and / or ribonucleic acid (RNA) to determine individual nucleotide bases of nucleic-acid sequences (e.g., sequencing by synthesis). The sequencing device 114 may analyze nucleic-acid segments and / or oligonucleotides extracted from samples to generate nucleotide reads and / or other data utilizing computer implemented methods and systems described herein either directly or indirectly on the sequencing device 114. More particularly, the sequencing device 114 may receive and analyze, within nucleotide-sample slides (e.g., flow cells), nucleic-acid sequences extracted from samples. The sequencing device 114 may utilize SBS to sequence nucleic-acid segments into nucleotide reads.

[0020] As further indicated by FIG. 1, the server device(s) 102 may generate, receive, analyze, store, and / or transmit digital data, such as data for determining nucleotide-base calls or sequencing nucleic-acid polymers. As shown in FIG. 1, the sequencing device 114 may send (and the server device(s) 102 may receive) call data from the sequencing device 114. The server device(s) 102 may also communicate with the client device 108. In particular, the server device(s) 102 may send data to the client device 108, including a variant call file or other information indicating nucleotide-base calls, sequencing metrics, error data, and / or other metrics associated with a nucleotide-base call, such as a call quality, a genotype, and / or a genotype quality.

[0021] The server device(s) 102 may comprise a distributed collection of servers where the server device(s) 102 include a number of server devices distributed across the network 112 and located in the same or different physical locations. Further, the server device(s) 102 may comprise a content server, an application server, a communication server, a web-hosting server, or another type of server.

[0022] As further shown in FIG. 1 , the server device(s) 102 may include a bioinformatics subsystem 104. The bioinformatics subsystem 104 may analyze base call data, such as sequencing metrics received from the sequencing device 114, to determine nucleotide base sequences for nucleic-acid polymers. For example, the bioinformatics subsystem 104 may receive raw data from the sequencing device 114 and may determine a nucleotide base sequence for a nucleic-acid segment. The bioinformatics subsystem 104 may determine the sequences of nucleotide bases in DNA and / or RNA segments or oligonucleotides. In addition to processing and determining sequences for biological samples, the bioinformatics subsystem 104 may generate a variant call file indicating one or more nucleotide-base calls for one or more genomic coordinates. The bioinformatics subsystem 104 may comprise one or more iterative genotypers capable of performing analysis, file generation, data aggregation, compression, and / or serialization, as described herein. The iterative genotypers may also be distributed across the bioinformatics subsystem 104, the client subsystem 110, and / or the database 116 to perform as described herein. The types of analysis performed by the bioinformatics subsystem 104 may include primary analysis, secondary analysis, and / or tertiary analysis. Primary analysis may include the detection and analysis of raw data (e.g., raw data in a file format, such as FASTA or FASTQ, comprising reads generated from the sequencing device 114). Secondary analysis may include analysis from reads (e.g., raw data in a file format, such as FASTA or FASTQ, comprising reads generated from the sequencing device 114) to VCF files (e.g., gVCFs) that are aligned against a reference genome (or denovo assembled), including graph genome mapping and aligning, and / or the performance of the calling, including germline variant calling functionality. The tertiary analysis may include variant annotation, variant filtering, prioritization, data visualization, and / or clinical report generation, which may be based on variant analysis from secondary analysis. Some may refer to variant calling, based on primary and / or secondary analysis, and / or other forms of downstream analysis for which conclusions may be made relating to genetic features, expression, or mutations as tertiary analysis. The primary analysis, the second analysis, and / or the tertiary analysis may be performed on the sequencing device 114 or off of the sequencing device 114 (e.g., at the server device(s) 110).

[0023] The client device 108 may generate, store, receive, and / or send digital data. In particular, the client device 108 may receive sequencing metrics from the sequencing device 114. Furthermore, the client device 108 may communicate with the server device(s) 102 to receive avariant call file comprising nucleotide base calls and / or other metrics, such as a call-quality, a genotype indication, and / or a genotype quality. The client device 108 may present and / or display information pertaining to the nucleotide-base call within a graphical user interface to a user associated with the client device 108. For example, the client device 108 may present an interface that includes a visualization for operating a genotyper and / or the results of analysis generated therefrom.

[0024] The client device 108 illustrated in FIG. 1 may comprise various types of client devices. In examples, the client device 108 may include non-mobile devices, such as desktop computers or servers, or other types of client devices. In other examples, the client device 108 may include mobile devices, such as laptops, tablets, mobile telephones, or smartphones.

[0025] As further illustrated in FIG. 1, the client device 108 may include a client subsystem 110. The client subsystem 110 may be a web application or a native application stored and executed on the client device 108 (e.g., a mobile application, desktop application). The client subsystem 110 may include instructions that (when executed) cause the client device 108 to receive data from the sequencing device 114 and present, for display at the client device 108, data from a variant call fde. Furthermore, the client subsystem 110 may instruct the client device 108 to display a visualization of contribution measures for sequencing metrics of a nucleotide- base call.

[0026] As further illustrated in FIG. 1, the environment 100 may include a database 116. The database 116 can store information such as variant call files, sample nucleotide sequences, nucleotide reads, nucleotide-base calls, and sequencing metrics. The server device(s) 102, the client device 108, and / or the sequencing device 114 may communicate with the database 116 (e.g., via the network 112) to store and / or access information, such as variant call files, sample nucleotide sequences, nucleotide reads, nucleotide-base calls, and / or sequencing metrics. The database 116 may also store one or more machine learning (ML) models and / or datasets associated with those ML models (e.g., feature sets, output files, etc.).

[0027] The environment 100 may be included in a local network or local high-performance computing (HPC) system. For example, the iterative genotypers described herein may be executed on one or more client devices 108 (e.g., as a part of the client subsystem 110), one or more server devices 102 (e.g., as a part of the bioinformatics subsystem 104), and / or one or more sequencing devices 114 in a local network or HPC system. The environment 100 may beincluded in a cloud computing environment comprising a plurality of server devices, such as server device(s) 102, having software and / or data distributed thereon. For example, the iterative genotypers described herein may be executed on one or more server devices 102 in a cloud computing environment. The bioinformatics subsystem 104 may be implemented to execute the gVCF genotypers described herein, and may be distributed across server devices 102 having access to the database 116 via the network 112 in a cloud-based computing system.

[0028] Though FIG. 1 illustrates the components of environment 100 communicating via the network 112, it will be appreciated that the components of environment 100 may communicate directly with each other, for example, bypassing the network 112. For example, the client device 108 may communicate directly with the sequencing device 114.

[0029] FIG. 2 illustrates a block diagram of an example computing device 200. One or more computing devices such as the computing device 200 may implement one or more features of the iterative genotyper described herein and / or the bioinformatics subsystem 104. One or more computing devices such as the computing device 200 may operate as a client device, server device, and / or a sequencing device, as described herein. As shown by FIG. 2, the computing device 200 may comprise a processor 202, a memory 204, a storage device 206, an I / O interface 208, and / or a communication interface 210, which may be communicatively coupled by way of a communication infrastructure 212. It will be appreciated that the computing device 200 may include fewer or more components than those shown in FIG. 2.

[0030] The processor 202 may include hardware for executing instructions, such as those making up a computer program. In examples, to execute instructions for dynamically modifying workflows, the processor 202 may retrieve (or fetch) the instructions from an internal register, an internal cache, the memory 204, or the storage device 206 and decode and execute the instructions. The memory 204 may be a volatile or non-volatile memory used for storing data, metadata, and programs for execution by the processor(s). The storage device 206 may include storage, such as a hard disk, flash disk drive, or other digital storage device, for storing data or instructions for performing the methods described herein. The memory 204 may have stored thereon computer-readable or machine-readable instructions for performing one or more processes or methods described herein.

[0031] The I / O interface 208 may allow a user to provide input to, receive output from, and / or otherwise transfer data to and receive data from the computing device 200. The I / O interface208 may include a mouse, a keypad or a keyboard, a touch screen, a camera, an optical scanner, network interface, modem, other known I / O devices or a combination of such I / O interfaces.The I / O interface 208 may include one or more devices for presenting output to a user, including, but not limited to, a graphics engine, a display (e.g., a display screen), one or more output drivers (e.g., display drivers), one or more audio speakers, and one or more audio drivers. The I / O interface 208 may be configured to provide graphical data to a display for presentation to a user. The graphical data may be representative of one or more graphical user interfaces and / or any other graphical content.

[0032] The communication interface 210 may include hardware, software, or both. In any event, the communication interface 210 may provide one or more interfaces for communication (such as, for example, packet-based communication) between the computing device 200 and one or more other computing devices or networks. As an example, and not by way of limitation, the communication interface 210 may include a network interface controller (NIC) or network adapter for communicating with an Ethernet or other wire-based network or a wireless NIC (WNIC) or wireless adapter for communicating with a wireless network, such as a WI-FI.

[0033] Additionally, the communication interface 210 may facilitate communications with various types of wired or wireless networks. The communication interface 210 may also facilitate communications using various communication protocols. The communication interface 210 may also include hardware, software, or both that couples components of the computing device 200 to each other. For example, the communication interface 210 may use one or more networks and / or protocols to enable a plurality of computing devices connected by a particular infrastructure to communicate with each other to perform one or more aspects of the processes described herein. To illustrate, the sequencing process may allow a plurality of devices (e.g., a client device, sequencing device, and server device(s)) to exchange information such as sequencing data and error notifications.

[0034] The computing devices, systems, and portions thereof described herein may be implemented to assist in genome sequencing, including genomic variant calling and processing of genotyping files, as described herein. In one example, the computing devices, systems, and portions thereof described herein may be used for processing cohort level sequencing data. For example, one or more computing devices may be configured to receive one or more genomic variant call files associated with one or more samples. An example of the genomic variant callfiles may be a genome variant call format (gVCF) file. gVCF genotypers may be implemented by one or more computing devices to perform sequencing analysis on the genomic variant call files, as described herein. For example, a gVCF genotyper may be implemented at one or more computing devices and may receive sample data from one or more sequencing devices to generate a gVCF file. The gVCF file may be a digital file generated in a publicly available standard text format that includes a number of predefined fields of summary information related to a sample, such as genome variant data, related to the sample to which the gVCF file corresponds. The summary information in the gVCF file may include genome variant data about the variants and non-variant genomic blocks at specific genomic coordinates, including metainformation lines, a header line, and data lines where each data line contains information about a single nucleotide-base call (e.g., a single variant). The genome variant data in the gVCF file may include one or more nucleotide-base calls (e.g., variant calls) along with other information pertaining to the nucleotide-base calls e.g., variant calls, quality, mapping alignment, and other metrics). The variant calls may be per allele, such that each variant allele is identified.

[0035] As gVCF genotypers may be focused on per-sample high performance sequencing analysis, the genome variant data in each gVCF file may include information related to a single sample from a single sequencing run, a sequencing cycle, or multiple sequencing runs at a sequencing device. A gVCF genotyper may take the summary information from a batch of multiple gVCF files, which each correspond to a single sample from a sequencing device, and analyze the information in an attempt to identify aggregate genome variant data and / or other genotype data from the batch of gVCF files. Due to the size of the gVCF files and large number of gVCFs to aggregate, the storage of the gVCF files may utilize large amounts of memory resources and analysis of the gVCF files may utilize large amounts of processing resources. In one example, the gVCF files may utilize 48 threads and / or 250 GB of memory for analysis on a local server or other computing device. In an example cloud instance (e.g., on an AWS c5d,18xlarge), the gVCF files may utilize 72 threads and / or 144 GB of memory for analysis. In another example of the processing power that may be utilized for analyzing a single gVCF file, the gVCF file may utilize 1 thread and / or 4 GB of memory, though additional resources may be implemented.

[0036] It may be desirable to analyze samples that have been taken by multiple sequencing devices across multiple sites and / or pipelines to identify aggregate genome variant data and / orother genotype data across sequencing devices to build a dataset having a larger aggregate dataset of genome variant data and other genotype data from large cohorts. Aggregating and joining genome variant data and other genotype data in batches of gVCF files from multiple sequencing devices at multiple sites may scale beyond several thousand samples. The size of a single gVCF file may be relatively large when compared to other file types that may be used for storing and analyzing data, as the gVCF files are comprised of text fields that may occupy a greater amount of storage than other field types. For example, a gVCF file with a 3 Ox sequencing depth may be 4 to 5 GB. The number of gVCF files may include thousands or hundreds of thousands of gVCF files.

[0037] The standard gVCF files that are used for sequencing analysis may include a plurality of fields (e.g., which comprise genome variant data) that are utilized in the standard format. For example, the plurality of fields in the gVCF file(s) may include a genotype (GT) field, a genotype quality (GQ) field, a minimum of genotype quality (GQX) field, a filtered base call depth (DP) field, a base calls filtered from input (DPF) field, an allelic depth (AD) field, a read depth associated with Indel (DPI) field, a mapping qualities (MQ) field, a filter (FT) field, a quality (QL) field, a phred-scaled genotype likelihood (PL) field, and a reference allele, one or more alternate alleles + genotype (GT) field, a contig name (CHROM), the start and end position of the record (POS, END), the reference allele sequence (REF), and / or the sequence of one or more alternate alleles (ALT).

[0038] The GT field may be encoded as allele values separated by a delimiter (e.g., either of / or |). The allele values may be zero for the reference allele (e.g., what is in the REF field), one (1) for the first allele listed in ALT, two (2) for the second allele list in ALT, and so on. For diploid calls, allele value examples may include 0 / 1, 1 | 0, or 1 / 2, etc. For haploid calls, e.g., on Y, male nonpseudoautosomal X, or mitochondrion, one (e.g., only one) allele value may be given; a triploid call may be 0 / 0 / 1. If a call (e.g., allele call) cannot be made for a sample at a given locus, ‘ .’ may be specified for each missing allele 5 in the GT field (for example ‘. / .’ for a diploid genotype and for haploid genotype). The separator ° I : may represent an unphased genotype. The separator ° | : may represent a phased genotype. The REF field may indicate one or more reference bases (e.g., A, C, G, T, N). The ALT field may be a comma separated list of alternate non-reference alleles called on at least one of the samples.

[0039] The GQ field may indicate conditional genotype quality encoded as a phred quality -lOlo lO probability of the genotype call being wrong (e.g., the error rate of genotyping), conditioned on the site’s being variant (Integer). As indicated, the GQ field may store a log transformed probability of an error that indicates whether a genotype call is correct or incorrect.

[0040] The GQX field may indicate the genotype quality assuming variant position or assuming non-variant position.

[0041] The filtered base call depth may be used for site genotyping. The DP field may be an integer value.

[0042] The AD field may indicate allelic depths for the ref and alt alleles, for example, in the order listed. For Indels, this value may include reads that confidently support each allele, for example, having a posterior probability of 0.999 or higher that read contains indicated allele vs. other (e.g., all other) intersecting Indel alleles.

[0043] The DPI field may be taken from a site preceding the Indel.

[0044] The MQ field may indicate RMS mapping quality. The MQ field may be an integer value. The system may identify or generate mapping quality scores for nucleobase calls at genomic coordinates, where a MAPQ score represents -lOloglO probability of the read mapping position being wrong e.g., the error rate of read mapping position), rounded to the nearest integer. As indicated, the MAPQ score may include a log transformed probability of an error that indicates whether a read mapping position is correct or incorrect. Additionally or alternatively, the system may determine soft-clipping metrics for sample nucleic-acid sequences by, for example, determining a total number of soft-clipped nucleobases spanning a genomic coordinate. If any of the fields is missing, it is replaced with the missing value. For example, if the FORMAT is GT:GQ:DP:HQ then 0 0 : . : 23 : 23, 34 indicates that GQ is missing. Trailing fields can be dropped (with the exception of the GT field, which should always be present if specified in the FORMAT field). See below for additional genotype fields used to encode structural variants. Additional Genotype fields can be defined in the meta-information. However, software support for such fields is not guaranteed.

[0045] The FT field may include a sample genotype filter indicating if this genotype was “called” (e.g., similar in concept to the FILTER field). PASS may be used in the FT field to indicate that all filters have been passed. A semicolon-separated list of codes may be used in the FT field to indicate one or more filters that fail. A period may be used in the FT field toindicate that filters have not been applied. These FT field values may be described in the metainformation in the same way as FILTERS. The FT field may be a string with no whitespace or semicolons permitted.

[0046] The QL field may indicate a phred-scaled quality score (e.g., for the assertion made in the ALT field - lOloglO probability of the call in ALT being wrong (e.g., the error rate of ALT)). The phred-scaled quality score may be a measure of base calling accuracy. The phred- scaled quality score may be used to assess the accuracy of a sequencing platform. The phred- scaled quality score may indicate the probability that a given base is called incorrectly by the sequencer. If the ALT field indicates (no variant) then the phred-scaled quality score may be calculated as -lOloglO prob (variant), and if ALT is not the phred-scaled quality score may be calculated as -lOloglO prob(no variant). If unknown, the MISSING value must be specified.

[0047] The PL field may indicate the phred-scaled genotype likelihoods rounded to the closest integer (e.g., and otherwise defined precisely as the GL field). The PL field may be an integer value.

[0048] The RAGT (reference alternate genotype) field may be a combination of “Reference allele + Alternate allele + Genotype” from one sample in one genomic region or at one genomic position. The RAGT field may be critical in calculation of allele frequency, allele count, and / or normalization of different variant alleles across large number of samples. The RAGT statistics (e.g., in the RAGT field) may indicate how many samples with what kinds of alleles are at each genomic position. The RAGT statistics (e.g., in the RAGT field) may indicate whether a particular region in the genome is difficult to sequence and / or genotype. The RAGT values may include key information on allele frequency and variants. The RAGT may hash the reference allele, the alternate allele, and the genotype of one sample at one genomic position or registration into one key. The key may be population-genotype-specific.

[0049] An iterative genotyper is described herein that may be implemented on one or more computing devices (e.g., one or more server devices 102 and / or client device 108 shown in FIG. 1) to iteratively aggregate variant data across a plurality (e.g., hundreds of thousands or millions) of samples from different projects, sites, and / or pipelines. The iterative genotyper may be implemented in a joint genotyping pipeline configured to efficiently aggregate variants from hundreds of thousands to millions of samples. FIG. 3 depicts an example process 300 that may be implemented by an iterative genotyper on one or more computing devices (e.g., such as theserver device 102, the client device 108, and / or the sequencing device 114 shown in FIG. 1) to perform iterative gVCF genotyping, as described herein. The iterative genotyper may be implemented by at least one processor via computer-readable or machine-readable instructions stored in and accessed from memory. For example, the process 300, or portions thereof, may be implemented in computer-executable instructions that are stored in memory and executed by a processor at the one or more computing devices. The process 300, or portions thereof, may be performed to efficiently aggregate and / or store variant data. The process 300, or portions thereof, may be performed at a single computing device or may be distributed across multiple computing devices (e.g., multiple server(s), sequencing device(s) and / or client computing device(s)). The process 300, or portions thereof, may reduce the amount of processing resources and / or storage space (e.g., memory) used by the computing device(s) during a variant analysis. The process 300 may enable high throughput for sequencing analysis of gVCF files.

[0050] The iterative genotyper may receive per sample gVCF data as an input. For example, the iterative genotyper may receive a gVCF file for each sample. As illustrated in FIG. 3, at 305, one or more batches of sample data may be received in batches of gVCF files, such as gVCF files 302A, 302B, 302C, 302D in batch 1 and / or gVCF files 304A, 304B, 304C, 304D in batch 2. Each of the one or more batches may have a batch size of thousands or hundreds of thousands of individuals (or samples). Each gVCF file 302A, 302B, 302C, 302D in batch 1 may include a summary of genome variant data that relates to a single sample. The gVCF files 302A, 302B, 302C, 302D may be received from a first sequencing device at a first site. Each gVCF file 304A, 304B, 304C, 304D in batch 2 may include a summary of genome variant data that relates to a single sample. The gVCF files 304A, 304B, 304C, 304D may be received from a second sequencing device at a second site. Genome variant data may be used interchangeably with genotype variant data. Each of the batches (e.g., batch 1 and batch 2) may be generated for different projects and / or from different sites and / or pipelines. The genome variant data collected for a cohort or population may be included in cohort level sequencing data.

[0051] The genome variant data may include a NON REF field that indicates whether there are any possible alternative allele at a respective location. The genome variant data may include a QU AL field that may include a phred-scaled probability that a site has no variant and that may be used to determine which variants may be given a PASS / FAIL call status. The genome variant data may include a FILTER field that may include an indication of applied filters or alack of applied filters. The FILTER field may indicate the filter(s) that the variant fails to pass, or the value pass if each of the variant passed each filter. The genome variant data may include a LowQUAL and / or LowGQ field that indicates filter values indicating the quality of hom-ref or variant genotype calls in the gVCF file. The genome variant data may include an allelic depth (AD) field that indicates allelic depths for the ref and alt alleles in the order listed. The genome variant data may include a DP field that indicates an approximate read depth (e.g., reads with MQ=255 or with bad mates may be filtered). The genome variant data may include a GQ field that indicates a genotype quality associated with the record. The genome variant data may include a GT field that indicates a genotype associated with the record. As described herein, the GT field may also indicate a missing genotype or that a call cannot be made for a sample at a given locus. The genome variant data may include a MIN DP value that indicates a minimum DP observed within the GVCF block. The genome variant data may include a PGT value that indicates physical phasing haplotype information that describes how the alternate alleles are phased in relation to one another. The genome variant data may include a PID value that indicates physical phasing ID information, where each unique ID within a given sample (e.g., but not across samples) connects records within a phasing group. The genome variant data may include a PL value that indicates normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification. The genome variant data may include an SB value that indicates persample component statistics which comprise the Fisher's Exact Test to detect strand bias.

[0052] The genome variant data may include a minimum amount of coverage observed at any one site within a block of records. The genome variant data may include a BaseQRankSum value that indicates a Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities. The genome variant data may include a ClippingRankSum value that indicates a Z-score From Wilcoxon rank sum test of Alt vs. Ref number of hard clipped bases. The genome variant data may include a DP value that indicates an approximate read depth (e.g., some reads may have been filtered). The genome variant data may include a DS value that indicates which samples were downsampled. The genome variant data may include an END value that indicates a stop position of the interval. The genome variant data may include an ExcessHet value that indicates a Phred-scaled p-value for exact test of excess heterozygosity. The genome variant data may include an InbreedingCoeff value that indicates an inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hardy-Weinberg expectation. Thegenome variant data may include a maximum likelihood expectation allele count (MLEAC) value that indicates a maximum likelihood expectation (MLE) for the allele counts (e.g., not necessarily the same as the AC), for each ALT allele, in the same order as listed. The genome variant data may include a maximum likelihood expectation allele frequency (MLEAF) that indicates a maximum likelihood expectation (MLE) for the allele frequency (e.g., not necessarily the same as the AF), for each ALT allele, in the same order as listed. The genome variant data may include an MQ value that indicates an RMS Mapping Quality. The genome variant data may include an MQRankSum value that indicates a Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities. The genome variant data may include a RAW value that indicates raw data for RMS mapping quality. The genome variant data may include a ReadPosRankSum value that indicates a Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias.

[0053] The iterative genotyper may aggregate the gVCF fdes over an iterative process into cohort files 311, 312 that include cohort data for each of the gVCF files in batch of data, followed by generation of an msVCF file from cohort files of the same batch. In one example, the VCF, gVCF, and / or msVCF files may be generated using the ILLUMINA DRAGEN iterative genotyper for performing secondary analysis. The cohort files 311, 312 may be iteratively generated and aggregated as batches of gVCF files are received. The gVCF genotyper may also, or alternatively, write a multi-sample VCF (msVCF) file that includes cohort level sequencing data without writing separate cohort files iteratively.

[0054] As shown in FIG. 3, the iterative genotyper may use the process 300 to aggregate a batch of gVCF files into msVCF files that include cohort level sequencing data for batches of samples at different sites. At 310, the batch of gVCF files 302A, 302B, 302C, 302D may be converted into a cohort file 312 that includes cohort level sequencing data generated from the batch of gVCF files 302A, 302B, 302C, 302D. At 315, the batch of gVCF files 304A, 304B, 304C, 304D may be converted into a cohort file 311 that includes cohort level sequencing data generated from the batch of gVCF files 304A, 304B, 304C, 304D.

[0055] The cohort files 311, 312 may include cohort level sequencing data from the gVCF files from which they are generated. The cohort file 311 may include cohort level sequencing data from the respective batches of gVCF files 304A, 304B, 304C, 304D from which it is generated. The cohort file 312 may include cohort level sequencing data from the respectivebatches of gVCF files 302A, 302B, 302C, 302D from which it is generated. The cohort level sequencing data in each of the cohort files 311, 312 may include data from one or more fields of the gVCF files from which the cohort file is generated. For example, the cohort files 311, 312 may each include at least the AD field, the GT field, the GQ field, the FT field, the QL field, the PL field, DP (e.g., total depths of a ref, alt, and / or non ref allele), Mapping Quality Ranksum (MQR), and / or the RAGT field from the batches of gVCF files from which they are generated. The iterative genotyper may be configured to store single sample or multi-sample level variant data in the cohort files 311, 312. For example, the iterative genotyper may generate a cohort file 311, 312 for each batch of samples (e.g., batch 1, batch 2, etc.) by aggregating data from gVCF files associated with each of the samples in a respective batch for the subset of fields in the cohort files 311, 312.

[0056] The iterative genotyper may generate at 330, 335 a respective multi-sample VCF (msVCF) file 331, 332 for each batch of samples. The msVCF files 331, 332 may be generated using a respective one of the cohort files 311, 312. For example, the msVCF file 332 may be generated (e.g., for batch 1) using the cohort file 312. The msVCF file 332 may be generated (e.g., for batch 2) using the cohort file 311. The msVCF files 331, 332 may include batchspecific genome variant data and / or other genotype data. The specific genome variant data in each msVCF file 331, 332 may identify variant types (e.g., single nucleotide variants (SNVs), insertions or deletions (Indels), and / or other variant types) for each variant allele in the msVCF file 331, 332. The genome variant data in each of the msVCF files 331, 332 may include a GT field that indicates a genotype associated with the record. As described herein, the GT field may also indicate a missing genotype or that a call cannot be made for a sample at a given locus.

[0057] The msVCF files 331, 332 may include other genotype data, such as global census data that includes genome variant data and / or other genotype data identified from samples taken at other sites Examples of the global census data may include the global set of variant sites and the variant alleles (e.g., variant types for each variant allele), genotypes called for each sample, read depth for variant alleles and hom-ref positions, and / or likelihoods and quality scores for alleles. Different msVCF files 331, 332 may be generated for respective batches (e.g., batch 1 and batch 2) of samples for processing and scalability. For example, the processing of the data received in the gVCF files (e.g, the gVCF files 302A, 302B, 302C, 302D, 304A, 304B, 304C, 304D) may be executed in parallel across multiple samples. For example, the cohort file 312may be generated from batch 1 of the gVCF files 302A, 302B, 302C, 02D, while the cohort file 311 may be generated from batch 2 of the gVCF files 304 A, 304B, 304C, 304D using parallel processing as described herein.

[0058] The iterative genotyper may merge the respective msVCF files 331, 332 for each batch into an aggregate msVCF file 333 with each of the samples and data of the respective msVCF files 331, 332. It will be appreciated that although the example shown in FIG. 3 has two batches, the process 300 is not limited to two batches. Instead, the process 300 may receive any number of batches of gVCF files from sequencing devices at any number of sites. The iterative genotyper may maintain and store in the aggregate msVCF file 333 summary statistics for variant sites and non-variant genomic blocks for a single sample and / or across multiple samples (e.g., multiple samples on a single sequencer). The iterative genotyper may maintain the RAGT statistics and the variant statistics from the each of the batches of gVCF files, the RAGT statistics and the variant statistics from the batches of gVCF files (e.g., for samples included in the per batch msVCF file), and other metrics e.g., non-RAGT metrics) from other fields, such as FT, GQ, AD, PL (with sample identifiers). In one example, the iterative genotyper may keep track of the RAGT statistics for different alleles that have been identified from the different sample data received from the batches of samples at different sites. The RAGT field may be a summary field used to determine each of the variant statistics (e.g., POS, REF, ALT, INFO fields) in summary fields in other filed (e.g., in the msVCF file described herein). From RAGT statistics, site statistics may be stored that may include the number of samples with genotype or GT, a number of samples without genotype (but with coverage), a number of samples without coverage, the total number of samples, the alternative allele count, and allele order, and / or mapping between global alternate alleles and per sample ALT alleles as stored in the RAGT. The iterative genotyper may maintain and store in the aggregate msVCF file 333 one or more secondary metrics, for example, such as allele counts, percentage of samples without sequencing coverage, percentage of samples without reliable genotypes, and / or the like. The aggregate msVCF file 333 may aggregate a large number (e.g., thousands) of batches. For example, the aggregate msVCF file 333 may include hundreds of thousands of samples to millions of samples and data associated therewith.

[0059] In an example, the msVCF files 331, 332 and the aggregate msVCF file 333 may include one or more global statistics, for example, allele frequencies, a read depth for variantalleles and / or hom-ref positions, a number of samples with or without genotypes, and / or a number of samples without coverage. The msVCF fdes 331, 332 and the aggregate msVCF fde 333 may include a likelihood and / or quality score for alleles. Similar statistics among samples in the batch are also included. An msVCF file may be generated for each batch of samples, to facilitate data read / write efficiency, data transfer, storage, and / or querying.

[0060] The fields in each of the msVCF files 331, 332 and the aggregate msVCF file 333 may include a summary of genotypes and / or samples. For example, the fields in the msVCF files 331, 332 and the aggregate msVCF file 333 may include an allele count in genotypes, a total number of alleles in called genotypes, a total number of samples, a total number of samples with called genotypes, a total number of samples with unknown genotypes, and / or a total number of samples with no coverage. The fields in each of the msVCF files 331, 332 and the aggregate msVCF file 333 may also include INFO fields, which may include features such as inbreeding count (IC), allelic balance metrics, Hardy-Weinberg Equilibrium P-value (HWE), and / or other features.

[0061] As will be understood, the DNA sequencing process may have a random component. A sequencing run may not yield any information (e.g., coverage from sequencing reads) for a set of positions of the genome. This may be captured by a “no coverage” metric. The number of these positions may be used to distinguish from other positions. For example, the variant information or statistics may include variant data indicating that a position has a DNA mutation, sequence reference materials (e.g., hom-ref data) indicating that a position has no DNA mutation, or no-coverage data that indicates that the variant information of the position is unknown. The gVCF variant files may include this information for a specific sample. The gVCF file may include, for each position in the genome, the variant information indicating that the position has a DNA mutation, hom-ref data indicating that the position has no DNA mutation, or no-coverage data that indicates that the variant information of the position is unknown.

[0062] The msVCF file may be a standardized public data format which enables users to run down-stream analysis, for example, genome- wide association studies (or GWAS), imputation and phasing, gene burden analysis, rare variant discovery, population specific allele frequencies, population substructure analysis, and estimating pathogenicity of mutation and classifying deleterious / benign variants in clinical analysis. The msVCF file may enable access to andquerying of aggregated variant data (e.g., in a batch). The msVCF file may include variant site level information both within the current batch and global across all batches of samples from sequencing devices at each site. The variant site level information may include, allele count, allele frequency, the total number of samples, the number of samples without coverage, the number of samples with coverage but not reliable genotypes and the number of samples with genotype. A user may prefer to use the msVCF file to access and query the aggregated variant data, for example, instead of other file formats. An example msVCF file of a batch may be represented by the following,

[0063] The iterative genotyper may continue to update the aggregate msVCF file 333 for each batch of samples to aggregate cohort level sequencing data in the aggregate msVCF file 333. The aggregate msVCF file 333 may include the variants and alleles discovered in each of the samples from the each of the batches received from genotyping devices across sites. Using the aggregate msVCF file 333, the iterative genotyper may provide an iterative population-based analysis option to jointly analyze samples from unrelated individuals.

[0064] The iterative genotyper may use the example process 300 to aggregate cohort level sample data from multiple sequencing devices at multiple sites into an existing data file, such as the aggregate msVCF file 333. The iterative genotyper may use the example process 300 toincrementally aggregate newly available genome variant data and other genotype data from batches of sample data with genome variant data and other genotype data available in prior batches. For example, the iterative genotyper may incrementally aggregate the genome variant data and / or other genotype data as it becomes available. The iterative genotyper may be scalable and may be executed on a plurality of computing platforms, for example, such as a cloud platform, a high performance cluster, and / or a single server. Though an example is provided for generating and storing cohort level sequencing data in a file type, such as the msVCF file, the same or similar data may be stored in another format, such as a database or other file structure, multiple files (e.g., a file for each genomic region), and / or a distributed file format. The use of the database or similar file structure may utilize queries to access the desired information from the file.

[0065] The aggregate msVCF file 333 may be used to jointly perform secondary analysis or tertiary analysis (e.g., analyze variants) on cohort level sequencing data from unrelated individuals. In next-generation sequencing (NGS)-based platforms, it is important to ensure quality control in the underlying data on which the analysis is being performed. For example, in performing secondary analysis or tertiary analysis on the cohort level sequencing data in the aggregate msVCF file 333, which may include hundreds of thousand samples to millions of samples and data associated therewith, there may be challenges in ensuring the quality of the results of such analysis, and particularly in comparing call sets generated from batches of data aggregated from different projects, sites, and / or pipelines. Data having a lower threshold of quality in the underlying dataset has been known to create problems in the results that are generated therefrom, such as false positives and false negatives e.g, creating too many or too few genotype calls), low genotyping consistency, and / or low genotyping rate or genotyping completion rate, for example. The genotype completion rate or genotyping rate may refer to the percentage of samples that are called (as either homozygous major, minor or heterozygous) at a given locus. Genotyping consistency may be particularly important in performing analysis on cohort level sequencing data related to trios and / or monozygotic twins. Sequencing data related to twins and trios may have biological constraints that need to be accounted for in the cohort level sequencing data.

[0066] The genotyping rate and / or genotyping consistency may vary due to a number of factors. For example, when aggregating cohort level sequencing data into msVCF files, somevariant types may have a consistent homogeneous genotype, while other variant types (e.g, Idels, etc.) may have a large percentage of samples (e.g., up to about half) with a missing genotype in the cohort level sequencing data. Even though there may be missing genotypes in regions of the sequencing data, the genotyper may make genotyping calls based on an estimated level of confidence using other portions of the cohort level sequencing data. For example, if an amount of sequencing coverage by the sequence reference materials (e.g., hom-ref data) and / or read depth is above a threshold, the genotyper may make a variant call in one sample even though there may be a missing genotype rate above a threshold estimate across all samples. As certain variant types (e.g., Indels) may be affected by higher missing genotype rates than other variant types (e.g., SNPs), certain variant types may be more difficult to genotype accurately than others. In smaller datasets, this may have less of an effect on the results. However, in msVCF files with cohort level sequencing data comprising hundreds of thousands to millions of samples, this may have a more significant impact. As the variant type (e.g, Indels) may indicate lower level of confidence in the genotyping rate than other variant types (e.g, SNPs), the variant type may be considered in determining the level of quality associated with a given variant call. When analyzing cohort level sequencing data, it can be difficult to ensure quality control for specific types of variants, particularly in different mapped regions of the genome due to noise or low confidence of mapping quality. It may also be difficult to distinguish between relatively good calls and relatively bad calls, given the location on the genome at which the calls are being made.

[0067] Additional factors may influence genotyping rate and / or genotyping consistency on cohort level sequencing data. For example, there are regions of the genome that are known to be difficult regions for making variant calls, as they may be difficult to map and / or may be inconsistent when genotyping. When performing analysis on sequencing data, one option to ensure quality control in variant calling is to obtain and sequence reference materials (e.g., hom- ref) and use them as a truth data set to measure variant calling performance against a truth dataset for the genotype data. The National Institute of Standards and Technology (NIST) has developed a publicly available set of sequencing data for seven human genomes that is characterized to allow a high threshold level of confidence when performing variant calling. An alternative or complementary approach may be to develop a predictive method that estimates sequencing performance for each sample. Variant calling may be impacted by both stochasticand systematic variability. For small variants, such as single nucleotide variants (SNVs) and short Indels, stochastic variability has a very minimal impact on variant calling performance in most of a genome (e.g., a genome that is sequenced to >30x coverage). However, genomic regions that are systematically impacted by lower quality - such as elevated error rates, low mapping quality or depth anomalies — may fail to deliver consistently accurate variant calls even for SNVs and Indels.

[0068] Due to a lack of truth data that can accurately be relied on and the number of variables that may influence the quality level of the results of the secondary or tertiary analysis performed, the threshold for identifying high confidence regions and low confidence regions of the genome may be difficult to identify and may not be particularly well defined, particularly when performing population-based analysis that jointly performs analysis on cohort level sequencing data from unrelated individuals. The cohort level sequencing data in the regions of the genome for which sequence reference materials are available as truth data for variant calling or other downstream analysis may be high confidence regions. As sequence reference materials are only available as truth data for a limited number ethnically homogenous samples e.g., 7 samples) and may not necessarily predict genome-wide performance in future samples, the cohort level sequencing data in the regions of the genome for which sequence reference materials are unavailable as a truth data may include relatively lower confidence regions.

[0069] Excluding variants based solely on reference features may eliminate many high-quality variant calls in relatively lower confidence regions. However, relying on the wrong data in lower confidence regions of the genome when performing secondary or tertiary analysis may lead to systematic errors. In one example, Genome in a Bottle (e.g., a consortium hosted by NIST and dedicated to characterization of benchmark human genomes) defines several classes of difficult regions of the genome based on additional characteristics of the genome, such as low mapability, segmental duplications, long tandem repeats or highly repetitive regions, and / or regions with extreme GC content. Highly repetitive regions having poor mapping quality and / or homopolymers may result in low base accuracy. This knowledge of highly repetitive regions has been used to classify the genome into high and low confidence regions. Though these classifications can be helpful in classifying regions of the genome, they are not a perfect representation of the actual performance within these regions. For example, a large segmentalduplication may be comprised of regions of high and low similarity, leading to very different variant calling accuracy.

[0070] Due to the lack of confidence in performing secondary or tertiary analysis in certain regions of the genome, there may be variation in the results (e.g., genotyping rate, consistency, etc.) based on different types of analysis performed on the same or similar cohort level data sets. This may be due to a difficulty in correctly mapping the reads in these regions or making high confidence calls of variants and genotypes, that are aggregated by genotyper. There may be different algorithms run on genotypers that may result in fewer variant calls in low confidence regions of the genome, while other forms of analysis may identify a greater number of variant calls in low confidence regions of the genome. One example algorithm run on a genotyper may include the Genome Analysis Toolkit (GATK), which may calls raw variants for each sample read. The GATK may analyze the variants against known variants, and apply a calibration procedure to compute a false discovery rate for each variant. Variants may be flagged as homozygous (1 / 1) or heterozygous (0 / 1) in a VCF file sample column. Other types of algorithms may be similarly run on the genotyper for making variant calls. When making these variant calls, it is also important to be consistent in the genotyping rate and / or genotyping consistency to ensure that the genotyping rate and / or consistency does not vary by more than a threshold in low confidence regions to high confidence regions.

[0071] Embodiments are described herein for implementing an ML model to assist with quality control of secondary and / or tertiary analysis (e.g., variant calling) performed by genotypers on cohort level sequencing data aggregated from different projects, sites, and / or pipelines. The ML model may be implemented to reduce the effects of missing genotypes in cohort level sequencing data. For example, the ML model may be trained on features of the cohort level data that relate to missing genotypes above a threshold. The ML model may also, or alternatively, be trained on features that related to different distributions of qualities across samples. For example, the ML model may be trained on cohort level sequencing data in certain regions of the genome which may be identified as high confidence regions and / or low confidence regions for mapping and / or variant calling.

[0072] FIG. 4 is a flowchart depicting an example procedure 400 for training an ML model to predict a machine learning site quality (MLSQ) score based on the feature values from the high and / or low confidence regions in cohort level sequencing data. The one or more portions of theprocedure 400 may be performed by one or more computing devices. For example, the one or more portions of the procedure 400 may be performed by one or more server devices (e.g., in a cloud-computing environment). One or more portions of the procedure 400 may be stored in memory as computer-readable or machine-readable instructions that may be executed by a processor of the one or more computing devices. One or more portions of the procedure 400 may be performed by one or more subsystems operating on one or more server devices. Though portions of the procedure 400 may be described herein as being performed by a computing device, one or more portions of the procedure 400 may be distributed across multiple computing devices, such as one or more client computing devices and / or one or more server computing devices.

[0073] The procedure 400 may begin at 402. As shown in FIG. 4, at 402 the computing device may receive cohort level sequencing data. The cohort level sequencing data may be received in one or more variant call fdes, in response to a query on a database, or from another storage medium. The cohort level sequencing data in the one or more variant fdes may indicate possible variants identified in the cohort level sequencing data. The variant call file may include one or more gVCF files (e.g, a batch of gVCF files) that may be used to generate an aggregate msVCF file that includes the cohort level sequencing data for performing population-scale variant analysis on sequencing data from different projects, sites, and / or pipelines. The cohort level sequencing data may be received in an aggregate msVCF file that has already been aggregated from msVCF files generated from respective batches of data. The cohort level sequencing data may include sequencing data for bi-allelic and / or multi-allelic sites or regions of the genome. Bi-allelic sites or regions may relate to or affect both alleles of a gene. Multi-allelic sites or regions may include one or more specific locus in a genome that includes three or more observed alleles (e.g, including the reference and two or more variant alleles).

[0074] The cohort level sequencing data may include hundreds of thousands of samples (e.g, 200,000 or more) to millions of samples (e.g., one to five million samples). The cohort level sequencing data may be included in a single file or resized into regions of the genome to effectively reduce the size for efficient handling / processing. If the cohort level sequencing data is separated into different files that include different regions of the genome, each file may be processed in parallel based on a number of available processing threads and / or computing devices implemented in processing.

[0075] The cohort level sequencing data may include sequencing data related to trios and / or monozygotic twins. Sequencing data related to twins and trios may have biological constraint about genotype quality that may assist with looking for consistency between samples, if the ML model is trained on such sequencing data.

[0076] At 404, the computing device may determine positive labeled training data and negative labeled training data from the cohort level sequencing data. The positive labeled training data may include the variants with a genotyping rate above a genotyping threshold from at least one high-confidence region in the genome. The positive labeled training data may identify regions of the genome with which there is a high-confidence (e.g., above a threshold) of mapping, genotyping rate, and / or distribution of other quality metrics in the cohort level sequencing data. The negative labeled training data may be selected from the variants with a genotyping rate below the genotyping threshold from at least one low-confidence region in the genome for variant calling. The negative labeled training data may identify regions of the genome with which there is a low-confidence (e.g., below a threshold) of mapping, genotyping rate, and / or distribution of quality metrics in the cohort level sequencing data. In one example, the positive labeled training data may include variants with a genotyping rate above a genotyping threshold and the negative labeled training data may include variants with a genotyping rate below the genotyping threshold. The genotyping threshold may be set to a relatively high threshold value (e.g., such as a 90% genotyping rate or higher, including a 95% genotyping rate, a 99% genotyping rate or 99.5% genotyping rate or higher). The genotyping threshold may depend on training data and / or the target accuracy of the prediction. A lower genotyping rate threshold may include more data for training, but the accuracy may be relatively lower.

[0077] The computing device may identify a joint high confidence truth region for the positive labeled training data and / or a joint low confidence truth region for the negative labeled training data in the cohort level sequencing data received at 402. The joint low confidence truth region may capture the majority of low confidence regions in genome for the cohort level sequencing data. The joint high confidence truth region may capture the majority of high confidence regions in genome for the cohort level sequencing data. The joint high confidence and joint low confidence regions may be selected for identifying respective positive and negative labeled training data for training the ML model. The number of samples selected for the positive labeled training data and the negative labeled training data may be maintained within a threshold toensure the ML model is not biased toward one dataset or the other. For example, training the model on training data selected from both regions may help maintain a consistent genotyping rate and ensure that the genotyping rate does not vary by more than a threshold.

[0078] High confidence regions may be regions that map better (e.g., have a mapping quality above a threshold) than other regions using the sequence reference materials. Low confidence regions in the cohort level sequencing data may be regions that map worse (e.g., have a mapping quality below a threshold) than other regions using the sequence reference materials. High confidence regions and / or low confidence regions in the cohort level sequencing data may be determined using sequence reference materials as a truth data set. For example, one or more of the 7 sequence reference samples may be used as a reference. The high-confidence regions of the cohort level sequencing data that are labeled as positive training data may include regions that may be mapped and / or called directly from the one or more of the 7 sequence reference samples. One or more of the sequence reference samples having issues when being used to map and / or call a particular region of the cohort sequencing data may cause the region to be identified as being a relatively lower confidence region and / or be labeled as negative training data. High confidence regions in the cohort level sequencing data may be regions that are easier to make a genotyping call (e.g., having a genotyping quality above a threshold) using the sequence reference materials. Low confidence regions in the cohort level sequencing data may be regions that are harder to make a genotyping call (e.g., having a genotyping quality below a threshold) using the sequence reference materials. For example, a high confidence region of the cohort level sequencing data may include regions that have a mapping quality score (e.g., MAPQ score) above a defined threshold and / or a genotyping quality score (e.g., QL score) above a defined threshold. A low confidence region of the cohort level sequencing data may include regions that have a mapping quality score (e.g., MAPQ score) below a defined threshold and / or a genotyping quality score (e.g., QL score) below a defined threshold. A high confidence region of the cohort level sequencing data may include regions that have a threshold read depth of coverage. A low confidence region of the cohort level sequencing data may include regions that have a threshold read depth of coverage. Low confidence regions may produce inconsistent results if sequenced multiple times. The computing device may identify non-uniform confidence (e.g., high and low confidence regions) across samples for genotyping and variant qualities. Low confidence regions may include regions having variant calls with relatively lower GQ and / or lower QUALvalues or GQ and / or QU AL values below a threshold. Low confidence regions may lead to inconsistent calls when sequencing related samples, using different variant detection algorithms or sequencing technologies.

[0079] High confidence regions may include a truth data set that is considered reliable, while low confidence regions may include truth data that is just a reference and / or may not be fully trusted as truth (e.g., both variant calls and genotype calls). The difference between high confidence regions and low confidence regions may be due to sequencing read length (e.g., short sequencing read length). Due to lack of truth data, it may be difficult to calibrate the variant calling and / or genotyping methods and, as such, the results may have a lower confidence and reproducibility in low confidence regions. Cohort level sequencing data metrics may be implemented to improve the genotyping rate and / or genotype quality of an individual sample, that may be unable to filter or keep when you have small sample sizes or at a single sample level. Identifying high confidence regions and low confidence regions may allow for calling variants in both high / low confidence regions, which may allow for high-sensitive calls in both high / low confidence regions to improve the genotyping rate and / or quality.

[0080] The training data may include sequencing data related to trios and / or monozygotic twins. The genotyping threshold for trios and / or monozygotic twins may include a consistency threshold, such that the positive labeled training data for trios and / or monozygotic twins includes sequencing data above a genotype consistency and negative labeled training data for trios and / or monozygotic twins includes sequencing data is below a genotype consistency. The genotyping threshold for trios and / or monozygotic twins may be the same as or higher than (e.g., such as a 95.5% or 100% genotyping rate) the genotyping threshold for other types of sequencing data, as the trios and / or monozygotic twins may already have biological constraints that need to be accounted for in the cohort level sequencing data.

[0081] At 406, the computing device may identify one or more feature values related to genotyping rate and / or distribution of quality metrics in the cohort level sequencing data. For example, the computing device may identify features for each sample in the positive and negative labeled training data of the cohort level sequencing data. The computing device may scrape the feature values of each feature for each sample in the positive and negative labeled training data. The computing device may query a database for the feature values, where adatabase is used. The one or more feature values may be a subset of the cohort level sequencing data associated with each sample.

[0082] The one or more feature values may be identified from a predefined set of features and / or feature values. Table 1 provides an example of a list of predefined features and feature values that may be defined for being identified for each allele in the cohort level sequencing data.TABLE 1As shown in Table 1, the features in the feature set may include features related to genotyping rate and / or distribution of quality metrics. The features may include a variant type value. For example, the variant type value may be determined from a stored value for each variant allele in the msVCF file in which the cohort level sequencing data may be stored. The variant type value may indicate a variant type (e.g., single nucleotide variants (SNVs), insertions or deletions (Indels), and / or other variant types) for each variant allele in the cohort level sequencing data. As described herein, different variant types may relate to different genotyping rates. For example, Indel alleles may indicate a relatively lower confidence variant or a less reliable variant than SNVs.

[0083] The features may also, or alternatively, include a value for a percentage of samples in the cohort level sequencing data with a missing genotype (e.g., GT). The percentage of sampleswith missing GT may be calculated from the GT field in the cohort level sequencing data in the msVCF file. As described herein, the GT field may indicate a missing genotype or that a call cannot be made for a sample at a given locus. From RAGT statistics, site statistics may also, or alternatively, be stored in the msVCF file that may include the number of samples with genotype or GT, a number of samples without genotype (but with coverage), and / or the total number of samples. The percentage value for samples with a missing genotype (e.g., GT) may be calculated from the number of samples with a genotype and the total number of samples. The higher the percentage of samples in a region with a missing GT may indicate a relatively lower confidence region or a less reliable region than regions with a higher percentage of samples with a missing genotype.

[0084] The features may also, or alternatively, include a value for a percentage of samples in the cohort level sequencing data with no coverage in the sequence reference materials (e.g, hom- ref data). From the RAGT statistics, site statistics may also, or alternatively, be stored in the msVCF file that may include a number of samples without coverage and / or the total number of samples. The percentage of samples in the cohort level sequencing data with no coverage may be calculated from the number of samples without coverage and the total number of samples. The higher the percentage of samples in a region with no coverage or a depth of coverage below a threshold may indicate a relatively lower confidence region or a less reliable region than regions with a lower percentage of samples with no coverage.

[0085] The features may also, or alternatively, include one or more values for a variant GQ histogram. For example, as shown in Table 1, there may be a predefined number (e.g., 15) of feature values selected and stored from the distribution of values in the variant GQ histogram. Each of the feature values may be within a range of values (e.g., 0-70). The GQ values may be binned to reduce the number of features (e.g., to the predefined number, such as 15) that have to be processed. In one example, instead of processing 71 values (e.g., 0-71), GQ values may be binned by every predefined number of values in the range of values (e.g., every 5 values, starting at value 1). The variant GQ histogram may indicate the genotype quality for variant sites plotted at different values across the calls of the cohort level sequencing data in the msVCF. The variant sites may include sites where a variant has been identified. Lower values in the distribution of values in the variant GQ histogram may indicate a relatively lower confidence region or a less reliable region than regions with higher values.

[0086] The features may also, or alternatively, include one or more values for a reference GQ (RGQ) histogram. The histogram of reference genotype quality scores may include scores for invariant sites called in the cohort level sequencing data in the msVCF. The invariant sites may include sites where no variation has been identified. GATK genotype quality scores at invariant sites (e.g., or RGQ) are calculated differently than those at variant sites and are thus not directly comparable. In one example, as shown in Table 1, there may be a predefined number (e.g., 21) of feature values selected and stored from the distribution of values in the RGQ histogram. Each of the feature values may be within a range of values (e.g., 0-100). The values may be binned to reduce the number of features (e.g., to the predefined number, such as 21) that have to be processed. The RGQ histogram may indicate the genotype quality for invariant sites plotted at different values across the calls of the cohort level sequencing data in the msVCF. Lower values in the distribution of values in the RGQ histogram may indicate a relatively lower confidence region or a less reliable region than regions with higher values.

[0087] Though certain features and / or feature values are provided as examples related to genotyping rate and / or distribution of quality metrics, a subset of these features and / or other features / feature values may be similarly implemented. For example, GQ or other features related to genotyping rate and / or distribution of quality may be added to the feature set. The variant GQ histogram and / or RGQ histogram values may be included in the subset of features, as lower range values of variant GQ histogram and / or RGQ histogram may be used as a greater indicator of genotyping rate and / or distribution of quality metrics than other features. For example, the lower range quality score range of variant GQ histogram and / or RGQ histogram may be informative to determine site importance. The lower range values of variant GQ histogram and / or RGQ histogram may be skewed toward low GQ individuals in the cohort level sequencing data. The percentage of samples with no coverage and / or the missing genotype may be included in a subset of features or other feature set as indicating regions that are difficult to map and / or subject to read length limitations. Analyzing data on short reads may impose limitations, which may be difficult to resolve in post-processing. As such, read length may also be used as a limitation on quality of a region. The features in the feature set may also, or alternatively, include QL, DP, MQR, IC, HWE, allelic depth (AD), localized allelic depth (LAD), sample filter (FT), local normalized, Phred-scaled likelihoods for genotypes (LPL), local GT values (LGT), MQRankSum, allele fractions of the local alt alleles (LAF), diploidfied (DF), allele count (AC),total number of alleles (AN), total number of samples (NS), total number of samples with called genotypes (NS_GT), total number of samples with unknown genotypes (NS_NOGT), total number of samples with no coverage (NS NODATA), ALT allele frequencies (AF), allelic balance among homozygotes (ABHom), allelic balance among heterozygotes (ABHet), P-value for allelic balance among heterozygotes (ABHetP), and / or other features.

[0088] At 408, an ML model may be trained to predict an MLSQ score based on the feature values from the positive and / or negative labeled cohort level sequencing data. The ML model may be trained using the labeling training data from the cohort level sequencing data in the joint high and / or low confidence regions to predict the MLSQ score for each variant allele identified in the cohort level sequencing data. The MLSQ score may be within a range (e.g., between 0 and 1) to indicate a relative level of quality associated with a variant that is being scored. The ML model may be trained by adjusting parameters and / or hyperparameters of the model based on the labeled training data. For example, the ML model may include weights of the one or more feature values in a predefined set of features related to genotyping rate or distribution of quality metrics, which may be adjusted based on the training data. Other parameters and / or hyperparameters may also be adjusted during training, as described herein.

[0089] The training data may train the ML model to identify good regions of the genome with which there is a high-confidence (e.g., above a threshold) of mapping and / or genotyping quality. The training data may train the ML model to identify bad regions of the genome with which there is a low-confidence (e.g., below a threshold) of mapping and / or genotyping quality. The ML model may be trained to predict the MLSQ score for each variant allele. The training may result in parameters (e.g., weight values associated with the one or more feature values in a predefined set of features related to genotyping rate or distribution of quality metrics) and / or hyperparameters defining the trained ML model. The MLSQ score may be trained on parameters related to a genotyping rate and / or distribution of quality to indicate the genotyping rate and / or distribution of quality associated with each variant allele in the cohort level sequencing data. The trained parameters and / or hyperparameters may be stored in memory at 410 for being implemented during production of the ML model.

[0090] FIG. 5 illustrates an example system environment and / or process 500 for training and / or implementing a supervised ML model 509, as described herein. For example, the bioinformatics subsystem 104 shown in FIG. 1, or portions thereof, may be implemented on oneor more computing devices to train and / or implement the supervised ML model 509. The supervised ML model 509 may be included in a trained filtering algorithm for filtering cohort level sequencing data to maintain quality of the cohort level sequencing data based on a predicted MLSQ score. The supervised ML model 509 may be trained and / or implemented during production on one or more computing devices, which may be the same or different computing devices during training and production of the ML model 509. For example, the supervised ML model 509 may be trained on one or more computing devices, such as one or more server devices 102 shown in FIG. 1 and the trained supervised ML model 509, or one or more trained parameters and / or hyperparameters 517 resulting from the training, may be implemented during production on one or more bioinformatics subsystems, such as the bioinformatics subsystem 104 shown in FIG. 1 (e.g., which may include one or more iterative genotypers).

[0091] As shown in FIG. 5, the supervised ML model 509 may include model data and one or more algorithms and / or functions trained to predict an MLSQ score for each variant allele in the cohort level sequencing data. The supervised ML model 509 may include one or more algorithms configured for supervised learning. Supervised learning may be implemented utilizing a supervised ML model 509 that is trained during a training process to determine a predictive model using known outcomes.

[0092] The supervised ML model 509 may be characterized by parameters and / or hyperparameters 517 that may be trained during the training process to configure the structure of the supervised ML model 509. The parameters may include values derived during the training process. The parameters may include a number of layers / nodes, weights(e.g., coefficients), and / or biases. The supervised ML model 509 may also include hyperparameters. The hyperparameters may include values used to control the learning process. The hyperparameters may include a learning rate, stride, a number of epochs, a batch size, a number of layers, a number of nodes in each layer, and / or other hyperparameters. Some may use certain parameters and hyperparameters interchangeably. In one example, the parameters and / or hyperparameters may be provided in a predefined library, such as XGBoost or another gradient boosting library.

[0093] The supervised ML model 509 may be trained during supervised learning by inputting training data 507a, 507b to the supervised ML model 509 and adjusting the parameters and / or hyperparameters 517 toward a known target output 515a, 515b while minimizing a loss or errorin the output generated by the supervised ML model 509. In one example, cross validation may be implemented, such that a subset of the training data may be used as the target output 515a, 515b to train the model and the remaining data may be used as the training data 507a, 507b. The training data 507a, 507b may be generated from pre-processing at 505a based on cohort level sequencing data 503a. Training data 507a, 507b, validation data, and / or testing data may be extracted from a dataset and formatted into the training data 507a, 507b, validation data, and / or test data for training, validation, and / or testing, respectively, the supervised ML model 509 during supervised learning. The dataset may include cohort level sequencing data 503a, which may be stored in an msVCF file, a database, or another storage medium.

[0094] The training data 507a, 507b may be input during a training process into the supervised ML model 509 in one or more formats, such as a tensor format, a vector format, an array format (e.g., including single-dimensional or multi-dimensional arrays) and / or another data format capable of being input into the supervised ML model 509. The training data 507a, 507b may be pre-processed at 505a to be format for being input into the supervised ML model 509. The preprocessing at 505a may include identifying positive labeled training data 507a and negative labeled training data 507b, as described herein. The positive labeled training data 507a may include labeled training data identified at joint high-confidence regions of the genome for mapping and / or variant calling of the cohort level sequencing data 503a. The negative labeled training data 507b may include labeled training data identified at joint low-confidence regions of the genome for mapping and / or variant calling of the cohort level sequencing data 503a. The positive labeled training data 507a may include variant alleles in the cohort-level sequencing data 503a having a genotyping rate above a predefined genotyping threshold. The negative labeled training data 507b may include variant alleles in the cohort-level sequencing data 503a having a genotyping rate below a predefined genotyping threshold. The pre-processing at 505a may include format changes or other types of processing in order to generate the training data 507a, 507b in a format for being input into the supervised ML model 509, as described herein.

[0095] The training data 507a, 507b may include cohort level sequencing data 503a related to trios and / or monozygotic twins. The genotyping threshold for trios and / or monozygotic twins may include a genotyping threshold, such that the positive labeled training data for trios and / or monozygotic twins includes sequencing data above a genotype consistency and negative labeled training data for trios and / or monozygotic twins includes sequencing data is below a genotypeconsistency. The genotyping threshold for trios and / or monozygotic twins may be the same or different than the genotyping threshold for other types of cohort level sequencing data. The positive labeled training data 507a may be trained toward a target output 515a. The negative labeled training data 507b may be trained toward a target output 515b.

[0096] In the configuration illustrated in FIG. 5, the cohort level sequencing data 503a may be pre-processed at 505a during the training procedure by identifying or extracting predefined feature values from the cohort level sequencing data 503a for generating the training data 507a, 507b. For example, the computing device may scrape the feature values for each sample in the cohort level sequencing data 503a from msVCF file, where such a file is used. The computing device may query a database for the feature values, where a database is used. As described herein the feature values may include one or more feature values indicating a variant type (e.g., SNP, Indel, etc.), one or more feature values indicating a percentage of samples with a missing GT, one or more feature values indicating samples with no coverage, one or more feature values of a variant GQ histogram (e.g., 0-70), and / or one or more feature values of an RGQ histogram (e.g., 0-100). The feature values identified in the cohort level sequencing data 503a may be used during the training process to train the ML model 509 for generating a predicted MLSQ score 519. The feature values identified in the cohort level sequencing data 503b may be used during implementation of the supervised ML model 509 to predict the MLSQ score 519 for use in applying the filter 521. The feature values may be identified for each variant allele in the cohort level sequencing data 503a. The feature values may be related to genotyping rate or distribution of quality metrics in the cohort level sequencing data.

[0097] Certain pre-processing steps may be performed at 505a to more efficiently manage the processing of the number of samples in the cohort level sequencing data 503a at the supervised ML model 509. The percentage of samples with missing GT and / or the percentage of samples with no coverage may be normalized into a fraction, rather than a raw count of samples in the cohort level sequencing data. The variant GQ histogram and / or the RGQ histogram may be normalized by frequency distribution prior to being input to the supervised ML model 509 as training data. The pre-processing (e.g., normalization) may be implemented, as the features may be independent of the absolute number of samples, which may vary from thousands to hundreds of thousands of samples.

[0098] The supervised ML model 509 may be trained by adjusting the parameters and / or hyperparameters 517 of the model based on the labeled training data 507a, 507b. For example, the supervised ML model 509 may include weights of the one or more feature values in the set of features related to genotyping rate or distribution of quality metrics, which may be adjusted based on the training data 507a, 507b. Other parameters and / or hyperparameters may also be adjusted during training.

[0099] The supervised ML model 509 may be trained to predict the MLSQ score for each variant allele in the cohort level sequencing data. The training may result in the parameters and / or hyperparameters 517 defining the trained ML model 509. The training data 507a, 507b may be labeled to teach the supervised ML model 509 to learn from the labeled data and to test the accuracy of the supervised ML model 509 for being implemented on unlabeled input data 508 during production / implementation of the supervised ML model 509, or similar supervised learning models utilizing similar parameters and / or hyperparameters 517. The training data 507a, 507b may be used to fit the parameters of the supervised ML model 509 using optimization functions, such as a loss or error function 513. Often the training data 507a, 507b may include pairs of input data and a corresponding target output 515a, 515b to which the parameters 517 may be trained to generate (e.g., within a threshold loss or error) a predicted output 519 during production and / or implementation. The training data 507a, 507b may be input into the supervised ML model 509 to train the supervised ML model 509 to predict the MLSQ score for each variant allele in the cohort level sequencing data. The trained parameters and / or hyperparameters 517may be stored in memory for being implemented during production of the ML model.

[0100] Supervised learning may be implemented for various types of supervised ML model 509, including algorithms that implement, for example, non-linear classification, logistic regression, neural networks (NNs), decision trees, Bayesian logics, random forests, and / or support vector machines (SVMs). NNs and Deep NNs (DNNs) are popular examples of algorithms utilized in ML models that may be trained using supervised learning. An example of the supervised ML model 509 may include a binary or JSON structure, which may follow a standard ML model specification of the XGBoost library that may allow the model to be read and written using different language bindings for training (e.g., Python, etc.) and implementation or prediction (e.g., C++, etc.). A binary model may be encrypted. A JSON model may enablevisualization of parameterizations and other information. Various examples of NNs include: perceptrons, multilayer perceptrons (MLPs), feed-forward NNs, fully-connected NNs, convolutional Neural Networks (CNNs), recurrent NNs (RNNs), long-short term memory (LSTM) NNs, and / or residual NNs (ResNets). A perceptron is a NN that includes a function that multiplies its input by a learned weight coefficient to generate an output value. A feed-forward NN is a NN that receives input at one or more nodes of an input layer and moves information in a direction through one or more hidden layers to one or more nodes of an output layer. In a feedforward NN, one or more nodes of a given layer may be connected to one or more nodes of another layer. A fully connected NN is a NN that includes an input layer, one or more hidden layers, and an output layer. In a fully connected NN, each node in a layer is connected to each node in another layer of the NN. An MLP is a fully connected class of feed-forward NNs. A CNN is a NN having one or more convolutional layers configured to perform a convolution. Various types of NNs may have elements that include one or more CNNs or convolutional layers, such as Generative Adversarial Networks (GANs). A GAN may include a generator submodel and a discriminator sub-model. The generator sub-model may be configured to receive input data and pass true and independently generated data to the discriminator sub-model. The discriminator sub-model may be configured to receive the true and independently generated data from the generator, discriminate the true and independently generated data, and provide feedback to the generator sub-model during training to improve the function of the generator sub-model in independently generating an output based on a received input. The GAN is a popular model for generating data types or data sequences, such as image data, audio data, and / or text, for example. An RNN is a NN that is recurrent in nature, as the nodes include feedback connections and an internal hidden state (e.g., memory) that allows output from nodes in the NN to affect subsequent input to the same nodes. LSTM NNs may be similar to RNNs in that the nodes have feedback connections and an internal hidden state (e.g., memory). However, the LSTM NNs may include additional gates to allow the LSTM NNs to learn longer-term dependencies between sequences of data. A ResNet is a NN that may include skip connections to skip one or more layers of the NN. Some NNs include one or more attention layers or functions to enhance or focus on some portions of the input data, while diminishing or de-emphasizing other portions.

[0101] The supervised ML model 509 may include layers of a similar type (e.g., convolutional layers, feed-forward layers, fully-connected layers, etc.) and / or having a similar or differentconfiguration (e.g., size, number of nodes, etc.) for each layer. The supervised ML model 509 may also, or alternatively, include one or more layers having different types or different subsets of NNs that may be interconnected for training and / or implementation, as described herein.

[0102] The training may be performed by initializing parameters and / or hyperparameters of the supervised ML model 509, generating and / or accessing the training data 507a, 507b, inputting the training data 507a, 507b into the supervised ML model 509, calculating the error or loss from the output of the supervised ML model 509 to a target output 515a, 515b via a loss function 513 (e.g., utilizing gradient descent and / or associated back propagation), and / or updating the parameters and / or hyperparameters 517. The target output 515a, 515b may include a trained label. The trained label (c.g., 0 or 1) in the training data 506a, 507b and / or the testing data may be used to train the model and / or fine-tune the parameters to avoid overfitting. The loss function 513 may be implemented using backpropagati on-based gradient updates and / or gradient descent techniques, such as Stochastic Gradient Descent (SGD), synchronous SGD, asynchronous SGD, batch gradient descent, and / or mini-batch gradient descent. Examples of loss or error functions may include functions for determining a squared-error loss, a mean squared error (MSE) loss, a mean absolute error loss, a mean absolute percentage error loss, a mean squared logarithmic error loss, and / or a cross-entropy loss. The training process may be iterated to update the parameters and / or hyperparameters 517 until an end condition is achieved. The end condition may be achieved when the output of the supervised ML model 509 is within a predefined threshold of the target output 515a, 515b.

[0103] The trained or fitted supervised ML model 509 may receive the validation data as input to evaluate the model fit on the training data set 507a, 507b, while tuning the hyperparameters 517 of the supervised ML model 509. The supervised ML model 509 may receive the test data to evaluate a final model fit on the training data set and to assess the performance of the supervised ML model 509. One or more of the training, validation, and / or testing may be performed during supervised learning for different types of supervised learning models. After the training and / or validation process is complete, the trained supervised ML model 509, or portions thereof, may be stored for being implemented by one or more devices. For example, the trained parameters and / or hyperparameters 517 may be stored, which may include weights and / or coefficients for each of the features in the predefined feature set for predicting an MLSQ score 519 based on received input data 508.

[0104] During implementation and / or production, the cohort level sequencing data 503b may be pre-processed at 505b to generate input data 508 for being input into the supervised ML model 509 in one or more formats, such as a tensor format, a vector format, an array format (e. ., including single-dimensional or multi-dimensional arrays) and / or another data format capable of being input into the supervised ML model 509. The format of the input data 508 may be the same or similar to the format of the training data training data 507a, 507b input during the training procedure. The cohort level sequencing data 503b may be pre-processed at 505b by identifying or extracting the predefined feature values in the feature set from the cohort level sequencing data 503b. For example, the computing device may scrape the feature values for each sample from msVCF file, where such a file is used. The computing device may query a database for the feature values, where a database is used. As described herein the feature values may include one or more feature values indicating a variant type (e.g, SNP, Indel, etc ), one or more feature values indicating a percentage of samples with a missing GT, one or more feature values indicating samples with no coverage, one or more feature values of a variant GQ histogram (e.g., 0-70), and / or one or more feature values of an RGQ histogram (e.g., 0-100).The feature values used in the input data 508 may be the same as or a subset of the feature values that are used in the training data 507a, 507b during the training procedure.

[0105] The feature values identified in the cohort level sequencing data 503b may be used during implementation of the supervised ML model 509 to predict the MLSQ score 519 for use in applying the filter 521. The feature values may be identified for each variant allele. The feature values may be related to genotyping rate or distribution of quality metrics in the cohort level sequencing data.

[0106] The predicted MLSQ score 519 may be generated by the trained supervised ML model 509, or portions thereof, for each variant allele in the cohort level sequencing data 503b. The predicted MLSQ score 519 for each variant allele may be implemented in other downstream algorithms or processes, as may be further described herein. For example, the MLSQ score 519 may be implemented in a filter 521 for filtering the variants identified in the cohort level sequencing data 503b. The trained supervised ML model 509 may be trained such that the number of variants that are filtered out is below a threshold (e.g., 3,000 variants out of 950,000 or 1,000,000, 3,000,000 variants out of 85,000,000 variants, or another threshold number or percentage) to preserve a threshold level of variant calls while reducing the missing genotypingrate and / or inconsistent variant calling. The MLSQ score 519 may be within a relative range of values (e.g., 0-1, in tenths or hundredths). A filter threshold may be configured to filter out variants identified as having an MLSQ score 519 below the filter threshold that fail to meet the quality level indicated by the filter threshold. The variants may be stored with an indication of the MLSQ score 519 and / or an indication of being above or below the filter threshold in an updated variant call file e.g., updated msVCF or database) in memory.

[0107] In one example, the updated variant call file may include the variants with the predicted MLSQ score that meet the genotyping threshold and the variants that have the predicted MLSQ score below the genotyping threshold may be removed from the variant call file. In another example, each of the variant alleles may be included in the updated variant call file with the predicted MLSQ score for each variant allele and / or an indicator that indicates whether the variants meet the filter threshold or are below the filter threshold, such that the updated variant call file includes the original genotyping data including each of the originally called variants for access to the original data. The filter threshold on the predicted MLSQ score may be set to a value and the variant alleles having a predicted MLSQ score below the filter threshold may be assigned an indictor in the file (e.g., a value indicating a low MLSQ score may be assigned).Each of the variants with a predicted MLSQ score below the filter threshold may be stored with a filter flag indicating that the variant is filtered and / or below the filter threshold. The preservation of the original data and the predicted MLSQ scores for each variant allele may allow for a user to update the filtering of certain variant alleles by updating the filter threshold value. The user may search for variants called in the original file to identify potentially filtered variants and adjust the filter threshold and / or use the predicted MLSQ score for other analysis on the particular variant of interest.

[0108] In one example, the supervised ML model 509 and / or the filter 521 may receive a single file (e.g., msVCF file) or database as input. In another example, the cohort level sequencing data may be separated into a number of files or databases that may be received as input separately. As described herein, the cohort level sequencing data in a file or database may include hundreds of thousands of samples to millions of samples, which may be resized to efficiently process the sequencing data. In one example, the cohort level sequencing data in a file or database may include between 300,000 samples to 5,000,000 samples and the sequencing data associated therewith. The files and / or databases may be separated by defined regions of thegenome to efficiently process the data in the files. The separate files and / or databases may be processed in parallel (e.g., depending on the number of threads and / or computing devices capable of processing). The input may be partitioned into files and / or databases based on the number of threads and / or computing devices being used for processing. The input may be partitioned into files according to respective regions of the genome. Each of the files and / or databases in a given partition may be processed by different instances (e.g., cloud instances or other instances on one or more computing devices) that are operating the supervised ML model 509 and / or the filter 521. The one or more computing devices may process each partition sequentially or in parallel, based on the number of threads, computing devices, and / or instances operating the supervised ML model 509 and / or the filter 521. The supervised ML model 509 may be similarly trained and / or operate on different genomic regions based on the same or similar parameters and / or hyperparameters 517. Thus, specific genomic regions may be sent to different computing devices, threads, or processing entities for being processed.

[0109] Updates may be performed to the supervised ML model 509. For example, updates to the ML model 509 may be performed to re-train the entire supervised ML model 509 or the updates may be performed efficiently, so as not to have to re-train the entire supervised ML model 509. For example, the collected feature values from the training data 507a, 507b may be stored in a separate file that may be utilized during an updated training procedure. The MLSQ score 519 may be reconfigured on different features or may be re-trained on a specific feature by merely updating the training process for a particular feature. Updated features and / or feature values may be implemented and the supervised ML model 509 may be retrained on the updated features and / or feature values. Again, because the cohort level sequencing data may include hundreds of thousands of samples to millions of samples and the data associated therewith, the supervised ML model 509 may be retrained on a particular updated feature or set of features and / or the updated feature values associated therewith. For example, the feature list may be updated or may be the same and the feature values may be updated based on updated training data. In another example, the training data may be updated by updating the labels on the positive and / or negative labeled training data. A lookup may be performed on the previous training data to perform updated training on a subset of the training data.

[0110] As the iterative genotyper may be implemented using different (e.g., updated) versions of software, the supervised ML model 509 may be trained on cohort level sequencing data (e.g.,msVCF files, databases, etc.) that is generated from different versions of the iterative genotyper. For example, different versions of the iterative genotyper may generate different forms of cohort sequencing data, which may include different data related to genotyping including genotyping rates and / or distributions of quality. Different versions of the iterative genotyper may generate different feature values for each of the features on which the supervised ML model 509 is trained. The supervised ML model 509 may be updated for different versions of the iterative genotyper by identifying different feature values in the cohort level sequencing data generated by different versions of the iterative genotyper and may re-train the supervised ML model 509 based on the updated feature values.[0U1] Different supervised ML models 509 may be trained for specific cohorts using training data 507a, 507b specific to the cohort. For example, the supervised ML models 509 may be customized for a specific cohort on which secondary and / or tertiary analysis may be performed. The supervised ML model 509 may be updated for improving the predicted MLSQ scores for a particular cohort and / or to add / remove individuals from the cohort.

[0112] FIG. 6 is a flowchart depicting an example procedure 600 for implementing an ML model for predicting an MLSQ score and applying a genotyping filter for assisting with maintaining quality control (e.g., genotyping rate and / or genotyping consistency) in variant calls in cohort level sequencing data. The one or more portions of the procedure 600 may be performed by one or more computing devices. For example, the one or more portions of the procedure 600 may be performed by one or more server devices. One or more portions of the procedure 600 may be stored in memory as computer-readable or machine-readable instructions that may be executed by a processor of the one or more computing devices. One or more portions of the procedure 600 may be performed by one or more subsystems operating on a server device. For example, one or more portions of the procedure may be performed by a bioinformatics subsystem implementing the iterative genotyper, as described herein. Though portions of the procedure 600 may be described herein as being performed by the bioinformatics subsystem implementing the iterative genotyper, the procedure 600, or portions thereof, may be performed by another subsystem operating on a computing device or distributed across multiple computing devices, such as one or more client computing devices and / or one or more server computing devices e.g., in a could computing environment).

[0113] The procedure 600 may begin at 602. As shown in FIG. 6, at 602 the iterative genotyper may receive cohort level sequencing data. The cohort level sequencing data may be received in one or more variant call fdes, in response to a query on a database, or from another storage medium. The cohort level sequencing data in the one or more variant fdes may indicate possible variants identified in the cohort level sequencing data. The variant call file may include one or more gVCF files (e.g., a batch of gVCF files) that may be used to generate an aggregate msVCF file that includes the cohort level sequencing data for performing population-scale variant analysis on sequencing data from different projects, sites, and / or pipelines. The cohort level sequencing data may be received in an aggregate msVCF file that has already been aggregated from gVCF files generated from respective batches of data. The cohort level sequencing data may include sequencing data for bi-allelic and / or multi-allelic sites or regions of the genome. Bi-allelic sites or regions may relate to or affect both alleles of a gene. Multi-allelic sites or regions may include one or more specific locus in a genome that includes three or more observed alleles (e.g., including the reference and two or more variant alleles).

[0114] The cohort level sequencing data may include hundreds of thousands of samples to millions of samples. The cohort level sequencing data may be included in a single file or resized into regions of the genome to effectively reduce the size for efficient handling / processing. If the cohort level sequencing data is separated into different files that include different regions of the genome, each file may be processed in parallel based on a number of available processing threads and / or computing devices implemented in processing.

[0115] At 604, the iterative genotyper may identify one or more feature values related to genotyping rate and / or distribution of quality metrics in the cohort level sequencing data. For example, the iterative genotyper may identify a predefined set of features for each sample in the cohort level sequencing data related to genotyping rate and / or distribution of quality metrics. The iterative genotyper may scrape the feature values of each feature for each sample from msVCF file, where such a file is used. The iterative genotyper may query a database for the feature values, where a database is used. The one or more feature values may be a subset of the cohort level sequencing data.

[0116] The one or more feature values may be predefined feature values defined in a predefined feature set. Table 1 described herein provides an example of a list of predefined features and feature values that may be defined for being identified for each allele in the cohortlevel sequencing data. As shown in Table 1 , the features may include features related to genotyping rate and / or distribution of quality metrics. The features may include one or more variant type values, one or more values for a percentage of samples in the cohort level sequencing data with a missing genotype (e.g., GT), one or more values for a percentage of samples in the cohort level sequencing data with no coverage, one or more values indicating a distribution of values in a variant GQ histogram, and / or one or more values indicating a distribution of values in an RGQ histogram. Though certain features and / or feature values are provided as examples related to genotyping rate and / or distribution of quality metrics, a subset of these features and / or other features / feature values may be similarly implemented. For example, GQ or other features may be added to the feature set that is related to genotyping rate and / or distribution of quality.

[0117] At 606, an ML model may be implemented by the iterative genotyper to predict an MLSQ score based on the feature values in the cohort level sequencing data. The ML model may predict the MLSQ score for each variant allele identified in the cohort level sequencing data. The MLSQ score may indicate a genotyping rate and / or distribution of quality associated with each variant allele in the cohort level sequencing data. The ML model may be previously trained on training data from high and / or low confidence regions of the genome. For example, the ML model may include weights of the one or more feature values in the predefined set of features related to genotyping rate or distribution of quality metrics.

[0118] At 608, the iterative genotyper may store the predicted MLSQ score with the cohort level sequencing data in memory. For example, the iterative genotyper may generate and store an updated variant call file (e.g., updated msVCF or database) with the predicted MLSQ score of each variant allele in memory. At 610, the iterative genotyper may apply a filter threshold to the predicted MLSQ score of each variant allele in the cohort level sequencing data. The filter threshold may be configured to filter out or flag variants with a predicted MLSQ score below the filter threshold, indicating a quality level below the threshold. The variants having a predicted MLSQ score above the filter threshold may be stored with an indication in an updated variant call file (e.g., updated msVCF or database) in memory.

[0119] In one example, the updated variant call file may include the variant alleles with the predicted MLSQ score above the filter threshold and the variant alleles that have the predicted MLSQ score below the filter threshold may be removed from the variant call file. In anotherexample, each of the variant alleles may be included in the updated variant call file with the predicted MLSQ score for each variant allele and / or an indicator that indicates whether the variant allele is above the filter threshold or are below the filter threshold, such that the updated variant call file includes the original genotyping data including each of the originally called variants for access to the original data. The filter threshold on the MLSQ score may be set to a value and the variant alleles below the filter threshold may be assigned an indictor in the file (e.g., a value indicating a low MLSQ score may be assigned). Each of the variants with an MLSQ score below the filter threshold may be stored with a filter flag indicating that the variant is filtered and / or below the filter threshold. The preservation of the original data and MLSQ scores for each variant allele may allow for a user to update the filtering of certain variant alleles by updating the filter threshold value. The user may search for variants called in the original file to identify potentially filtered variants and adjust the filter threshold and / or use the predicted MLSQ score for other analysis on the particular variant of interest.ML Model Training / Production Example

[0120] One particular example training procedure is provided below for training an ML model to predict an MLSQ score on an msVCF comprising cohort level sequencing data. An ML model was trained on a JSON format using XGBoost. The feature set illustrated as an example in Table 1 and the corresponding feature values were used on training data generated from ILLUMINA’s DRAGEN muti-sample VCF variant caller on UK biobank sequencing data.

[0121] Table 2 below provides characteristics of the training data that was used to train the ML model.TABLE 2As shown in Table 2, the ML model was trained on 500,000 positively labeled samples from high confidence regions of the genome and 500,000 negatively labeled samples from low confidence regions of the genome. The high confidence region of the genome included 1,416,948 samples of cohort level sequencing data. The low confidence region of the genome included 1,261,552 of cohort level sequencing data. The positively labeled samples had a genotyping rate above a genotyping threshold at or above 99.5 percent. The negatively labeled samples had a genotyping rate below the genotyping threshold of 99.5 percent.

[0122] As also shown in Table 2, the ML model was trained on 2 million samples of twin and trio sequencing data. The positively labeled twin and trio samples had a genotyping rate and a consistency of 100 percent. The negatively labeled twin and trio samples had a genotyping rate and a consistency of less than 100 percent. The training used 10-fild cross validations.

[0123] The trained ML model was applied to generate a predicted MLSQ score for each variant allele between 0 and 1. We applied a filter threshold to filter out variant alleles with a predicted MLSQ score of .1 or less. The filter may be adjusted based on the dataset, but if the filter is too stringent then high-quality variant calls may be lost. This filter removed about 3,000 variant alleles that were originally called by ILLUMINA’ s DRAGEN muti-sample VCF variant caller. The variant alleles that were filtered were compared to the original dataset of variant alleles originally called and the variant alleles that were filtered had a high rate of missing genotypes from low confidence regions of the genome, while a relatively lower percentage of the variant alleles originally called in the high-confidence regions of the genome remained. The use of the ML model and the MLSQ score for filtering may effectively reduce inconsistencies to improve the quality and / or number of variants that remain unfiltered.

[0124] Although features, elements, and functions are described above in particular combinations, a feature, element, or function is used alone or in any combination with the other features, elements, or functions. Various presently unforeseen or unanticipated alternatives,modifications, variations, or improvements may be subsequently made that are also intended to be encompassed by the following claims.

[0125] The methods described herein are implemented in a computer program, software, or firmware incorporated in a computer-readable medium for execution by a computer or processor. Examples of computer-readable media include electronic signals (transmitted over wired or wireless connections) and computer-readable storage media. Examples of computer-readable storage media include, but are not limited to, a read only memory (ROM), a random-access memory (RAM), removable disks, and optical media such as CD-ROM disks, and digital versatile disks (DVDs).

Claims

CLAIMSWhat is claimed is:

1. A computing system for implementing a trained machine learning model to assist in quality control of variant calls in cohort level sequencing data, the system comprising: a memory; and at least one processor configured to: receive at least one variant call file comprising cohort level sequencing data, wherein the cohort level sequencing data in the at least one variant call file indicates possible variants identified in the cohort level sequencing data; identify, based on a predefined set of features related to genotyping rate or distribution of quality metrics, input data from the at least one variant call file; and predict, via parameters of the trained machine learning model, a machine learning site quality (MLSQ) score for each variant in the at least one variant call file based on the input data from the at least one variant call file, wherein the trained machine learning model is trained based on each feature of the predefined set of features related to genotyping rate or distribution of quality metrics, and wherein the trained machine learning model is trained based on variants with a genotyping rate above a genotyping threshold from at least one high-confidence region in the genome for variant calling and variants with a genotyping rate below the genotyping threshold from at least one low- confidence region in the genome for variant calling; and store the predicted MLSQ score for each variant in the memory.

2. The computing system of claim 1, wherein the at least one processor is configured to: filter the variants identified in the cohort level sequencing data in the at least one variant call file having the MLSQ score below a filter threshold; and store the filtered variants in an updated variant call file in memory.

3. The computing system of claim 2, wherein the at least one processor is configured to store the filtered variants in the updated variant call file with an indication of being filtered based on the MLSQ score being below the filter threshold.

4. The computing system of claim 1, wherein the predefined set of features comprises at least one of a variant type, a number of samples with missing genotype, a number of samples with no coverage, a variant genotyping quality (GQ) histogram, or a reference GQ histogram.

5. The computing system of claim 1, wherein the trained machine learning model is trained based on: a genotype consistency in trios above a consistency threshold in the at least one high- confidence region in the genome; and a genotype consistency in trios below a consistency threshold in the at least one low- confidence region in the genome.

6. The computing system of claim 1, wherein the trained machine learning model is trained based on: a genotype consistency in monozygotic twins above a consistency threshold in the at least one high-confidence region in the genome; and a genotype consistency in monozygotic twins below a consistency threshold in the at least one low-confidence region in the genome.

7. The computing system of claim 1, wherein the variants identified in the cohort level sequencing data comprise each variant allele in the cohort sequencing data.

8. The computing system of claim 1, wherein the at least one processor is further configured to: partition the at least one variant call file comprising the cohort level sequencing data into a plurality of variant call files based on regions of the genome; and process each variant call file of the variant call files sequentially or in parallel based on a number of threads or instances operating the ML model.

9. The computing system of claim 1, wherein the at least one variant call file comprises between 300,000 samples and 5,000,000 samples.

10. The computing system of claim 1, wherein the at least one variant call file comprises the cohort level sequencing data corresponding to a specific genomic region.

11. At least one computer readable storage medium having instructions stored thereon that are configured to, when executed by at least one processor, cause the at least one processor to: receive at least one variant call file comprising cohort level sequencing data, wherein the cohort level sequencing data in the at least one variant call file indicates possible variants identified in the cohort level sequencing data; identify, based on a predefined set of features related to genotyping rate or distribution of quality metrics, input data from the at least one variant call file; and predict, via parameters of a trained machine learning model, a machine learning site quality (MLSQ) score for each variant in the at least one variant call file based on the input data from the at least one variant call file, wherein the trained machine learning model is trained based on each feature of the predefined set of features related to genotyping rate or distribution of quality metrics, and wherein the trained machine learning model is trained based on variants with a genotyping rate above a genotyping threshold from at least one high-confidence region in the genome for variant calling and variants with a genotyping rate below the genotyping threshold from at least one low-confidence region in the genome for variant calling; and store the predicted MLSQ score for each variant in the memory.

12. The at least one computer readable storage medium of claim 11, wherein the instructions are configured to cause the at least one processor to: filter the variants identified in the cohort level sequencing data in the at least one variant call file having the MLSQ score below a filter threshold; and store the filtered variants in an updated variant call file in memory.

13. The at least one computer readable storage medium of claim 12, wherein the instructions are configured to cause the at least one processor to store the filtered variants in the updated variant call file with an indication of being filtered based on the MLSQ score being below the filter threshold.

14. The at least one computer readable storage medium of claim 11, wherein the predefined set of features comprises at least one of a variant type, a number of samples with missing genotype, a number of samples with no coverage, a variant genotyping quality (GQ) histogram, or a reference GQ histogram.

15. The at least one computer readable storage medium of claim 11, wherein the trained machine learning model is trained based on: a genotype consistency in trios above a consistency threshold in the at least one high- confidence region in the genome; and a genotype consistency in trios below a consistency threshold in the at least one low- confidence region in the genome.

16. The at least one computer readable storage medium of claim 11, wherein the trained machine learning model is trained based on: a genotype consistency in monozygotic twins above a consistency threshold in the at least one high-confidence region in the genome; and a genotype consistency in monozygotic twins below a consistency threshold in the at least one low-confidence region in the genome.

17. The at least one computer readable storage medium of claim 11, wherein the variants identified in the cohort level sequencing data comprise each variant allele in the cohort sequencing data.

18. The at least one computer readable storage medium of claim 11, wherein the at least one processor is further configured to:partition the at least one variant call file comprising the cohort level sequencing data into a plurality of variant call files based on regions of the genome; and process each variant call file of the variant call files sequentially or in parallel based on a number of threads or instances operating the ML model.

19. The computing system of claim 1, wherein the at least one variant call file comprises between 300,000 samples and 5,000,000 samples.

20. The at least one computer readable storage medium of claim 11, wherein the at least one variant call file comprises the cohort level sequencing data corresponding to a specific genomic region.

Images