Systems and methods for iterative and scalable population-wide variant analysis
Patent Information
- Application Number
- JP2023579813
- Authority / Receiving Office
- JP · JP
- Patent Type
- Applications
- Current Assignee / Owner
- Priority Date
- 2022-03-31
- Filing Date
- 2022-12-15
- Publication Date
- 2025-12-16
AI Technical Summary
Existing sequencing technologies face challenges in efficiently aggregating and analyzing large-scale genomic variant data due to the size and complexity of genome variant call format (gVCF) files, leading to difficulties in timely delivery of population-wide variant analysis results.
An iterative process is implemented to incrementally aggregate genomic variant data in batches, using cohort and census files to reduce data size, and employ parallel processing and compression techniques to generate multi-sample variant calling files, enabling efficient storage and analysis of large datasets.
This approach reduces processing and storage requirements, allowing for high-throughput analysis of genomic data, facilitating scalable and timely population-wide variant analysis.
Smart Images

Figure 00000000_0000_ABST
Abstract
Description
[Technical field]
[0001] CROSS-REFERENCE TO RELATED APPLICATIONS This application claims priority to U.S. Provisional Patent Application No. 63 / 361,386, filed December 15, 2021, and U.S. Provisional Patent Application No. 63 / 326,227, filed March 31, 2022, which are incorporated by reference in their entireties herein. [Background technology]
[0002] A population-scale genomic experiment may include aggregating and / or merging variant-related data from many (e.g., hundreds of thousands) of samples. Existing sequencing analyses may focus on per-sample analysis. As sequencing throughput continues to increase, timely delivery of population-scale variant analysis results becomes increasingly desirable. Furthermore, in existing sequencing analyses, sample sequencing may be performed such that some data may not be available at any given time.
[0003] Existing sequencing analyses may use genome variant call format (gVCF) files. gVCF files store sequencing information for both variant and non-variant positions. gVCF files may allow for representation of genotype, annotation, and other information across all sites in a genome. A gVCF genotyper may be a population-based analysis tool that jointly analyzes variants from unrelated individuals. Summary of the Invention [Means for solving the problem]
[0004] Described herein are systems, methods, and apparatus for enabling an iterative process of incrementally aggregating available batches of sample data with previously available batches. One or more computing devices can be configured to receive one or more genomic variant call files associated with one or more samples. An example of a genomic variant call file can be a genomic variant call format (gVCF) file. Data in a genomic variant call file can include a list of variants and genomic blocks. The genomic variant call files can be received for batches of samples guided by sequencing devices at different sites.
[0005] Cohort and census files can be generated for each batch of samples. The cohort and census files can include a subset of the fields in the genomic variant call file and can include summary information for the batch of samples for the subset of fields. Census data in multiple census files generated from different batches of samples can be aggregated into a global census file.
[0006] The multi-sample variant call file may be generated based on the global census file, one or more cohort files, and one or more census files. The multi-sample variant call file may be stored in memory to perform a sequencing analysis on the data in the file. For example, one or more computing devices may be implemented to perform a genome-wide sequencing analysis using one or more multi-sample variant call files.
[0007] The genome variant call files may be processed using parallel processing on multiple computation nodes as described herein. Each batch of genome variant call files may be divided into equal sized shards to enable parallel processing of the genome variant call files. Each shard may be processed using one of the multiple computation nodes. The parallel processing may be implemented using multi-threading by region of sequence data. At least two computation nodes may be configured to implement at least two levels of parallelization to process, aggregate, and / or generate data for corresponding regions of sequence data. Each computation node may be configured to process a particular region. Each core may have a particular threading. A variable number of software threads may be implemented. One or more threads may be implemented per CPU core. For example, a single thread may be implemented by each CPU core. The number of threads implemented by each CPU core may be changed in response to user input.
[0008] One or more computing devices may be implemented as described herein to perform allele ordering and genotype reindexing. Each genomic variant call file may be associated with each sample of the plurality of samples that includes a reference alternate genotype (RAGT) statistic. The RAGT statistic may be used to identify a plurality of reference alleles and a plurality of alternate alleles associated with the sample in the genomic variant call file. The instances of each of the plurality of reference alleles and each of the plurality of alternate alleles may be summed for normalization to determine a number of unique alleles. A normalized reference allele may be selected from the plurality of reference alleles. The longest reference allele may be selected as the normalized reference allele. Other reference alleles of the plurality of reference alleles may be normalized by extending them to correspond to the normalized reference allele. The plurality of alternate alleles may be normalized by extending each alternate allele the same amount that each corresponding reference allele was extended. The normalized reference alleles and the normalized alternate alleles may be used to generate a multi-sample variant call file.
[0009] The one or more computing devices may be implemented to perform compression on data stored in the cohort file and / or the census file. A field of the cohort file or census file to be compressed may be identified. In one example, a plurality of reference alleles and a plurality of alternative alleles associated with a plurality of samples may be stored in a RAGT field that may be identified for compression. The one or more computing devices may be configured to determine which of the plurality of samples have a common reference allele and an alternative allele. The samples may be distributed into allele groups. Each allele group may include one or more samples having a common reference allele and an alternative allele. The one or more computing devices may be configured to select a binary value length based on the number of allele groups. The binary value length may be the shortest binary value length that can be used to uniquely identify each of the allele groups. The one or more computing devices may be configured to assign a unique binary value to each of the allele groups using the determined binary value length. The unique binary value for each of the allele groups may be stored in a bitmap that is used to encode the plurality of reference alleles and the plurality of alternative alleles into a bit array.
[0010] The genomic variant data in the cohort file and the census file may be aggregated into an output buffer including a certain number of buffer locations, as described herein. For example, one or more computing devices may be configured to receive a record of genomic variant data and determine whether the genomic variant data for the received record overlaps with one or more other previously stored records in the output buffer. If the genomic variant data for the received record does not overlap with a previously stored record in the buffer, the genomic variant data for the received record may be stored in one or more buffer locations in the buffer. If the genomic variant data for the received record overlaps with genomic variant data of another record, the buffer location including the overlapping portion of the record may be updated to include the genomic variant data for the received record. Any non-overlapping portions of the previously stored records may be copied and stored in sequential buffer locations along with the overlapping portions of the records. Any non-overlapping portions of the received records may be added to sequential buffer locations along with the overlapping portions of the records. [Brief description of the drawings]
[0011] [Figure 1] FIG. 1 shows a schematic diagram of the system environment (or "environment") in which an iterative genomic variant calling format (gVCF) genotyper may operate. [Diagram 2] 1 illustrates blocks of an exemplary computing device. [Diagram 3] FIG. 1 illustrates an exemplary process that may be implemented by one or more computing devices to perform iterative gVCF genotyping. [Figure 4] 1 is a flowchart of an exemplary method for processing a batch of sample data from a sequencing device to generate one or more multi-sample variant call format (msVCF) files. [Figure 5A]FIG. 1 illustrates an exemplary process performed by an iterative gVCF genotyper. [Figure 5B] FIG. 1 illustrates another exemplary process performed by an iterative gVCF genotyper. [Figure 6] 1 is a graphical representation of an exemplary cohort record generated for a batch of samples. [Figure 7] FIG. 1 illustrates an example iterative process that may be implemented by one or more computing devices to generate a multi-sample variant call file after receiving a gVCF file associated with a batch. [Figure 8] FIG. 1 illustrates an exemplary process that may be implemented by one or more computing devices to incorporate parallel processing in variant analysis. [Figure 9] FIG. 1 illustrates an exemplary process incorporating parallel processing across batches. [Figure 10] FIG. 1 illustrates an exemplary genomic data operator design. [Figure 11] FIG. 1 illustrates an exemplary process that may be implemented by one or more computing devices for re-indexing genotypes for an msVCF file. [Figure 12A] FIG. 1 illustrates an example process that may be implemented by one or more computing devices to encode cohort data. [Figure 12B] FIG. 1 illustrates an example process that may be implemented by one or more computing devices to encode census data. [Figure 13] FIG. 1 is a flow diagram illustrating an example method that may be implemented by one or more computing devices to perform compression of cohort and / or census data. [Figure 14] FIG. 1 illustrates an example process that may be implemented by one or more computing devices for serializing and compressing data. [Figure 15]FIG. 2 illustrates another example process that may be implemented by one or more computing devices for serializing and compressing data. [Figure 16] FIG. 1 illustrates an example process that may be implemented by one or more computing devices to add records, copy, and / or update data during buffer aggregation involving overlapping genomic regions. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS
[0012] 1 shows a schematic diagram of a system environment (or "environment") 100 in which an iterative genomic variant calling format (gVCF) genotyper may operate, as described herein. As shown, the environment 100 includes one or more server devices 102 connected to client devices 108 and sequencing devices 114 via a network 112.
[0013] 1, the server device 102, the client device 108, and the sequencing device 114 can communicate with each other via a network 112. The network 112 can include any suitable network with which computing devices can communicate. The network 112 can include a wired and / or wireless communication network. An exemplary wireless communication network can 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 Wi-Fi communication protocol, and / or another wireless communication protocol. In addition to or instead of communicating across the network 112, the server device 102, the client device 108, and / or the sequencing device 114 can communicate, bypass the network 112, and communicate directly with each other.
[0014] As shown by FIG. 1, the sequencing device 114 can include a device for sequencing a biological sample. The biological sample can include human and non-human DNA for determining individual nucleotide bases of a nucleic acid sequence (e.g., sequencing by synthesis). The sequencing device 114 can analyze nucleic acid segments and / or oligonucleotides extracted from the sample to generate nucleotide reads and / or other data utilizing the computer-implemented methods and systems described herein, either directly or indirectly on the sequencing device 114. More specifically, the sequencing device 114 can receive and analyze nucleic acid sequences extracted from the sample in a nucleotide-sample slide. The sequencing device 114 can utilize SBS to sequence the nucleic acid segments into nucleotide reads.
[0015] As further illustrated by FIG. 1, the server device 102 can generate, receive, analyze, store, and / or transmit digital data, such as data for determining nucleotide-base calls or sequencing a nucleic acid polymer. As illustrated in FIG. 1, the sequencing device 114 can transmit (and the server device 102 can receive) call data from the sequencing device 114. The server device 102 can also communicate with the client device 108. In particular, the server device 102 can transmit to the client device 108 data including variant call files or other information indicative of nucleotide-base calls, sequencing metrics, error data, and / or other metrics associated with the nucleotide-base calls, such as call quality, genotypes, and / or genotype quality.
[0016] The server device 102 includes a distributed collection of servers, where the server device 102 may include several server devices distributed across the network 112 and located at the same or different physical locations. Furthermore, the server device 102 may include a content server, an application server, a communication server, a web hosting server, or another type of server.
[0017] As further shown in FIG. 1, the server device 102 can include a sequencing system 104. The sequencing system 104 can analyze base call data, such as sequencing metrics, received from the sequencing device 114 to determine a nucleotide base sequence for a nucleic acid polymer. For example, the sequencing system 104 can receive raw data from the sequencing device 114 and determine a nucleotide base sequence for a nucleic acid segment. The sequencing system 104 can determine a sequence of nucleotide bases in a DNA and / or RNA segment or oligonucleotide. In addition to processing and determining a sequence for a biological sample, the sequencing system 104 can generate a variant call file indicating one or more nucleotide-base calls for one or more genomic coordinates. The sequencing system 104 can include one or more iterative gVCF genotypers capable of performing parsing, file generation, data aggregation, compression, and / or serialization as described herein. The iterative gVCF genotypers can also be distributed across the sequencing system 104, the sequencing application 110, and / or the database 116 and implemented as described herein. Primary analysis may be defined as calling individual nucleotide bases. Secondary analysis may be defined as alignment and / or assembly of DNA or RNA fragments (e.g., one or more nucleotide bases) by mapping reads to a reference genome. Tertiary analysis may be defined as variant identification / calling based on primary and / or secondary analysis. Primary, secondary, and / or tertiary analysis may be performed on-instrument or off-instrument.
[0018] The client device 108 can generate, store, receive, and / or transmit digital data. In particular, the client device 108 can receive sequencing metrics from the sequencing device 114. Additionally, the client device 108 can communicate with the server device 102 to receive variant call files including nucleotide base calls and / or other metrics such as call quality, genotype designation, and genotype quality. The client device 108 can present or display information regarding the nucleotide-base calls in a graphical user interface to a user associated with the client device 108. For example, the client device 108 can present a contribution measure interface that includes a visualization or depiction of various contribution measures associated with or attributable to individual sequencing metrics for a particular nucleotide-base call.
[0019] The client devices 108 shown in FIG. 1 may include various types of client devices. In examples, the client devices 108 may include non-mobile devices such as desktop computers or servers, or other types of client devices. In other examples, the client devices 108 may include mobile devices such as laptops, tablets, mobile phones, or smartphones. Additional details regarding the client devices 108 are discussed below with respect to FIG. 2.
[0020] 1, the client device 108 can include a sequencing application 110. The sequencing application 110 can be a web application or a native application (e.g., a mobile application, a desktop application) stored on and executed on the client device 108. The sequencing application 110 can include instructions that (when executed) cause the client device 108 to receive data from the sequencing device 114 and present data from a variant call file for display on the client device 108. Additionally, the sequencing application 110 can instruct the client device 108 to display a visualization of the contribution measures for the sequencing metrics of the nucleotide-base calls.
[0021] As further shown in FIG. 1, the environment 100 can 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 102, the client device 108, and / or the sequencing device 114 can 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 can also store one or more models, such as a call recalibration machine learning model and / or a call generation model.
[0022] The environment 100 may be included in a local network or a local high-performance computing (HPC) system. For example, the iterative gVCF genotyper described herein may be executed on one or more client devices 108 (e.g., as part of a sequencing application 110), one or more server devices 102 (e.g., as part of a sequencing system 104), and / or one or more sequencing devices 114 within the local network or HPC system. The environment 100 may be included in a cloud computing environment including multiple server devices, such as server devices 102 having distributed software and / or data. For example, the iterative gVCF genotyper described herein may be executed on one or more server devices 102 in the cloud computing environment. The sequencing system 104 may be implemented to execute the gVCF genotyper described herein and may be distributed across server devices 102 with access to a database 116 via a network 112 in a cloud-based computing system.
[0023] 1 illustrates components of environment 100 communicating over network 112, it will be understood that components of environment 100 may also communicate directly with one another, for example, bypassing network 112. For example, client device 108 may communicate directly with sequencing device 114.
[0024] 2 illustrates a block diagram of an exemplary computing device 200. One or more computing devices, such as computing device 200, may implement one or more features of the iterative gVCF genotyper and / or sequencing system 104 described herein. One or more computing devices, such as computing device 200, may operate as a client device, a server device, or a sequencing device, as described herein. As illustrated by FIG. 2, computing device 200 may include a processor 202, a memory 204, a storage device 206, a VO interface 208, and / or a communication interface 210, which may be communicatively coupled by a communication infrastructure 212. It will be understood that computing device 200 may include fewer or more components than those illustrated in FIG. 2.
[0025] The processor 202 may include hardware for executing instructions, such as instructions constituting a computer program. In an example, to execute instructions for dynamically modifying a workflow, the processor 202 may retrieve (or fetch) instructions from an internal register, an internal cache, the memory 204, or a storage device 206, decode and execute these instructions. The memory 204 may be a volatile or non-volatile memory used to store data, metadata, and programs for execution by the processor. The storage device 206 may include storage, such as a hard disk, a flash disk drive, or other digital storage device, for storing data or instructions for performing the methods described herein. The memory 204 may have computer-readable or machine-readable instructions stored therein for performing one or more processes or methods described herein.
[0026] The I / O interface 208 may enable 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 interface 208 may include a mouse, a keypad or keyboard, a touch screen, a camera, an optical scanner, a network interface, a 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., a display driver), 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 represent one or more graphical user interfaces and / or any other graphical content.
[0027] Communications interface 210 may include hardware, software, or both. In any case, communications interface 210 may provide one or more interfaces for communications (such as, for example, packet-based communications) between computing device 200 and one or more other computing devices or networks. By way of example and not limitation, communications interface 210 may include a network interface controller (NIC) or network adapter for communicating with an Ethernet or other wired-based network, or a wireless NIC (WNIC) or wireless adapter for communicating with a wireless network, such as Wi-Fi.
[0028] Additionally, the communication interface 210 can facilitate communication with various types of wired or wireless networks. The communication interface 210 can also facilitate communication using various communication protocols. The communication infrastructure 212 can also include hardware, software, or both that couples the components of the computing device 200 to one another. For example, the communication interface 210 can enable multiple computing devices connected by a particular infrastructure to communicate with one another using one or more networks and / or protocols to perform one or more aspects of the processes described herein. To illustrate, a sequencing process can enable multiple devices (e.g., a client device, a sequencing device, and a server device) to exchange information such as sequencing data and error notifications.
[0029] The computing devices, systems, and portions thereof described herein may be implemented to assist in genomic sequencing, including processing of genomic variant calling and genotyping files, as described herein. 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 a genomic variant call file may be a genomic variant call format (gVCF) file. As described herein, a gVCF genotyper may be implemented by one or more computing devices to perform sequencing analysis on the genomic variant call files. For example, a gVCF genotyper may be implemented in one or more computing devices and may receive sample data from one or more sequencing devices to generate a gVCF file. A gVCF file may be a digital file generated in a publicly available standard text format that includes several predefined fields of summary information related to a sample, such as genomic variant data related to the sample to which the gVCF file corresponds. The summary information in a gVCF file may include genomic variant data for variants and non-variant genomic blocks at a particular genomic coordinate, including a meta-information line, a header line, and a data line, where each data line includes information about a single nucleotide-base call (e.g., a single variant). The genomic variant data in a gVCF file can include one or more nucleotide-base calls (e.g., variant calls) along with other information about the nucleotide-base calls (e.g., variant calls, quality, mapping alignments, and other metrics).
[0030] Because the gVCF genotyper can focus on high performance sequencing analysis per sample, the genomic variant data in each gVCF file can include information related to a single sequencing run, sequencing cycle, or a single sample from multiple sequencing runs on a sequencing device. The gVCF genotyper can obtain summary information from a batch of multiple gVCF files, each corresponding to a single sample from a sequencing device, and analyze the information to identify aggregate genomic variant data and / or other genotype data from the batch of gVCF files. Depending on the size of the gVCF file, storage of the gVCF file may utilize a large amount of memory resources, and analysis of the gVCF file may utilize a large amount of processing resources. In one example, a gVCF file may utilize 48 threads and / or 250 GB of memory for analysis on a local server or other computing device. On an exemplary cloud instance (e.g., AWS c5d, 18xlarge), a gVCF file may utilize 72 threads and / or 144 GB of memory for analysis. As another example of the processing power that may be utilized to analyze a single gVCF file, the gVCF file may utilize one thread and / or 4 GB of memory, although additional resources may be implemented.
[0031] Additionally, it may be desirable to analyze samples acquired by multiple sequencing devices across multiple sites to identify aggregate genomic variant data and / or other genotypic data across the sequencing devices and build a dataset with a larger aggregated dataset of genomic variant data and / or other genotypic data. Indeed, in some exemplary embodiments, the gVCF genotyper feature may be used to implement Genome Analysis Toolkit (GATK) algorithms to aggregate and combine genomic variant data from large cohorts. Due to the size of the gVCF file, aggregating and combining genomic variant data and other genotypic data in a batch of gVCF files from multiple sequencing devices at multiple sites may be difficult to scale (e.g., more than several thousand samples). A batch of gVCF files from multiple sequencing devices may be difficult to store, analyze, and / or aggregate to identify aggregate genomic variant data and / or other genotypic data collected across the sequencing devices. Because the gVCF file is composed of text fields that may occupy a larger amount of storage than other field types, the size of a single gVCF file may be relatively large when compared to other file types that may be used to store and analyze data. For example, a gVCF file with a 30x sequencing depth may be 4-5 GB. As gVCF files are aggregated to perform analysis of data across batches or different sites in which samples have been received, storing and analyzing the data within the gVCF file may become increasingly difficult. The size of the gVCF file may also make it difficult to communicate batches of sample data from different sites and aggregate variant information or other summary information from the sites to identify variant information across sites.
[0032] To provide an example of the size of a standard gVCF file used in sequencing analysis, a number of fields (e.g., including exemplary genomic variant data) utilized in the standard format are described herein. For example, the number of fields in a gVCF file can include a genotype (GT) field, a genotype quality (GQ) field, a minimum genotype quality (GQX) field, a filtered base call depth (DP) field, a filtered base call from input (DPF) field, an allele depth (AD) field, a read depth associated with indels (DPI) field, a mapping quality (MQ) field, a filter (FT) field, a quality (QL) field, a Phred-scaled genotype likelihood (PL) field, as well as a reference allele, one or more alternative alleles + genotype (GT) field, a contig name (CHROM), a record start and end position (POS, END), a reference allele sequence (REF), and / or one or more alternative allele sequences (ALT).
[0033] The GT field may be encoded as allele values separated by a delimiter (e.g., either / or |). The allele value may be zero for a reference allele (e.g., one in the REF field), 1 for the first allele listed in the ALT, 2 for the second allele listed in the ALT, etc. For diploid calls, example allele values may include 0 / 1, 1|0, or 1 / 2, etc. For example, for a haploid call on Y, male non-pseudoautosomal X, or mitochondrial, an allele value of 1 (e.g., only one) may be given, and a triploid call may be 0 / 0 / 1. If a call (e.g., an allele call) cannot be made for a sample at a given locus, a "." may be specified for each missing allele 5 in the GT field (e.g., ". / ." for diploid and haploid genotypes). The separator ° / : may represent a non-phasic genotype. The separator °|: may represent a phasic 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 alternative non-reference alleles that were called in at least one of the samples.
[0034] The GQ field can indicate the conditional genotype quality, coded as the Phred quality-10 log10 probability of an incorrect genotype call (e.g., genotyping error rate), conditional on the site being a variant (an integer). As shown, the GQ field can store the log-transformed probability of error indicating that the genotype call is correct or incorrect.
[0035] The GQX field may indicate the genotype quality assuming a variant position or assuming a non-variant position.
[0036] The filtered base call depth can be used for site genotyping. The DP field can be an integer value.
[0037] The AD field may indicate the allele depth for the ref and alt alleles, e.g., in the order listed. For indels, this value may include reads that support each allele with confidence, e.g., reads where the posterior probability that the read contains the specified allele and other (e.g., all other) crossover indel alleles is 0.999 or greater.
[0038] The DPI field can be obtained from the site preceding the indel.
[0039] The MQ field may indicate the RMS mapping quality. The MQ field may be an integer value. The system may identify or generate a mapping quality score for a nucleobase call at a genomic coordinate, where the MAPQ score represents the -10log10 probability of a read mapping position being incorrect (e.g., the error rate of the read mapping position), rounded to the nearest integer. As shown, the MAPQ score may include a log-transformed probability of error indicating that the read mapping position is correct or incorrect. Additionally or alternatively, the system may determine a soft clipping metric for the sample nucleic acid sequence, for example, by determining the total number of soft-clipped nucleobases that span the genomic coordinate. If any of the fields are missing, it is replaced with a missing value. For example, if FORMAT is GT:GQ:DP:HQ, then 0|0:.:23:23, then 34 indicates that GQ is missing. Trailing fields may be dropped (except for 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 may be defined in the meta information, however software support for such fields is not guaranteed.
[0040] The FT field may contain a sample genotype filter indicating whether 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 were passed. A semicolon separated list of codes may be used in the FT field to indicate one or more filters that failed. A period "." may be used in the FT field to indicate that no filters were applied. These FT field values may be listed in the meta information in the same way as FILTERS. The FT field may be a string where spaces or semicolons are not allowed.
[0041] The QL field may indicate a Phred-scaled quality score (e.g., for an assertion made in the ALT field, the -10log10 probability of a call in the ALT being incorrect (e.g., the error rate of the 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), the Phred-scaled quality score may be calculated as -10log10 probability (variant), and if the ALT is not ".", the Phred-scaled quality score may be calculated as -10log10 probability (no variant). If unknown, a MISSING value must be specified.
[0042] The PL field may indicate a Phred-scaled genotype likelihood rounded to the nearest integer (eg, and otherwise precisely defined as the GL field). The PL field may be an integer value.
[0043] The RAGT (reference alternative genotype) field can be a combination of "reference allele + alternative allele + genotype" from one sample at one genomic region or one genomic location. The RAGT field can be important in normalizing allele frequencies, allele counts, and / or different variant alleles across many samples. The RAGT statistic (e.g., in the RAGT field) can indicate how many samples with what kind of alleles are at each genomic location. The RAGT statistic (e.g., in the RAGT field) can indicate whether a particular region in the genome is difficult to sequence and / or genotype. The RAGT value can include key information about allele frequencies and variants. RAGT can hash the reference allele, alternative allele, and genotype of one sample at one genomic location or accession into one key. The key can be population-genotype specific.
[0044] Described herein is an iterative gVCF genotyper that may be implemented on one or more computing devices (e.g., one or more server devices 102 and / or client devices 108 shown in FIG. 1 ) to iteratively aggregate variant data across a plurality of samples (e.g., hundreds of thousands or more). FIG. 3 illustrates an exemplary process 300 that may be implemented by an iterative gVCF genotyper on one or more computing devices (e.g., such as the server device 102, client device 108, and / or sequencing device 114 shown in FIG. 1 ) to perform iterative gVCF genotyping as described herein. The iterative gVCF genotyper may be implemented by a processor via computer-readable or machine-readable instructions stored in and accessed from a memory. For example, process 300, or portions thereof, may be implemented in computer-executable instructions stored in a memory and executed by a processor in one or more computing devices. Process 300, or portions thereof, may be implemented to efficiently aggregate and / or store variant data. Process 300, or portions thereof, may be performed on a single computing device or may be distributed across multiple computing devices (e.g., multiple servers, sequencing devices, and / or client computing devices). Process 300, or portions thereof, may reduce the amount of processing resources and / or storage space (e.g., memory) used by a computing device during variant analysis. Process 300 may enable high throughput for sequencing analysis of gVCF files.
[0045] The iterative gVCF genotyper can receive gVCF data per sample as input. For example, the iterative gVCF genotyper can receive a gVCF file for each sample. As described in FIG. 3, at 305, one or more batches of sample data can 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 can have a batch size of about 1000 individuals (or samples). The sample size within a batch can range from one to several thousand or other numbers. Each gVCF file 302A, 302B, 302C, 302D in batch 1 can include a summary of genomic variant data associated with 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 genomic variant data associated with a single sample. The gVCF files 304A, 304B, 304C, 304D may be received from a second sequencing device at a second site. Genomic variant data may be used interchangeably with genotypic variant data.
[0046] The genomic variant data may include a NON REF field indicating whether there are any possible alternative alleles at the respective position. The genomic variant data may include a LowQual field indicating whether the record is of low quality. The genomic variant data may include an allele depth (AD) field indicating the allele depth for the ref and alt alleles in the order listed. The genomic variant data may include a DP field indicating an approximate read depth (e.g., reads with MQ=255 or bad mates may be filtered). The genomic variant data may include a GQ field indicating the genotype quality associated with the record. The genomic variant data may include a GT field indicating the genotype associated with the record. The genomic variant data may include a MIN DP value indicating the minimum DP observed within the GVCF block. The genomic variant data may include a PGT value indicating physical phasing haplotype information describing how the alternative alleles are phased with respect to each other. The genomic variant data may include a PID value indicating physical phasing ID information, where each unique ID within a given sample connects records within a phasing group (e.g., not across samples). The genomic variant data can include PL values, which indicate the normalized, Phred-scaled likelihood for a genotype as defined in the VCF standard. The genomic variant data can include SB values, which indicate component statistics per sample, including Fisher's exact test for detecting strand bias.
[0047] The genomic variant data may include a minimum amount of coverage observed at any one site within a block of records. The genomic variant data may include a BaseQRankSum value indicating a Z-score from a Wilcoxon rank sum test of Alt vs. Ref base quality. The genomic variant data may include a ClippingRankSum value indicating a Z-score from a Wilcoxon rank sum test of Alt counts vs. Ref counts of hard-clipped bases. The genomic variant data may include a DP value indicating an approximate read depth (e.g., some reads may have been filtered). The genomic variant data may include a DS value indicating which samples were downsampled. The genomic variant data may include an END value indicating the stopping position of the interval. The genomic variant data may include an ExcessHet value indicating a Phred-scaled p-value for an exact test of excess heterozygosity. The genomic variant data may include an InbreedingCoeff value indicating the inbreeding coefficient as estimated from the genotype likelihoods per sample when compared against the Hardy-Weinberg expectation value. The genomic variant data can include a maximum likelihood expectation allele count (MLEAC) value indicating, for each ALT allele, a maximum likelihood expectation (MLE) for the allele count (e.g., not necessarily the same as AC), in the same order as listed. The genomic variant data can include a maximum likelihood expectation allele frequency (MLEAF), indicating, for each ALT allele, a maximum likelihood expectation (MLE) for the allele frequency (e.g., not necessarily the same as AF), in the same order as listed. The genomic variant data can include an MQ value indicating RMS mapping quality.The genomic variant data may include an MQRankSum value indicating a Z-score from a Wilcoxon rank sum test of Alt vs. Ref read mapping quality. The genomic variant data may include a RAW value indicating raw data for RMS mapping quality. The genomic variant data may include a ReadPosRankSum value indicating a Z-score from a Wilcoxon rank sum test of Alt vs. Ref read position bias.
[0048] An iterative gVCF genotyper can perform aggregation in an "end-to-end mode," for example, when a single batch of samples (e.g., only a single batch) is available. End-to-end mode may be implemented when there is one batch of samples (e.g., the first batch). In a "step-by-step" mode, gVCF files may be aggregated into cohort and census files, and subsequently msVCF files may be generated from the cohort and census files of the same batch. Aggregation may not be performed in an "end-to-end" mode, where the gVCF genotyper can write a multi-sample VCF file without writing cohort and census files. The "end-to-end" process may not be iterative, as cohort and census files may not be generated iteratively.
[0049] The iterative gVCF genotyper may perform aggregation in a "step-by-step mode," for example, when multiple batches of samples are available, such that files are generated and information is aggregated batch by batch. In step-by-step mode, when multiple batches of samples are available, the iterative gVCF genotyper may perform steps to aggregate gVCFs into cohort and census files batch by batch. Once aggregation is performed batch by batch, the census files from the batches may be aggregated into a global census file in a next step. An msVCF may be generated using the global census file, cohort file, and census file for each batch, as described herein.
[0050] 3, an iterative gVCF genotyper can aggregate a batch of gVCF files into a generated cohort file and a generated census file using process 300. At 310, a batch of gVCF files 302A, 302B, 302C, 302D can be converted into a cohort file 312 and a census file 314. At 315, a batch of gVCF files 304A, 304B, 304C, 304D can be converted into a cohort file 311 and a census file 313. The formats of the cohort files 311, 312 and the census files 313, 314 together can enable a user to extract variant data from the publicly available standard gVCF file format and ingest and compress the data for large scale aggregation.
[0051] The cohort files 311, 312 and the census files 313, 314 may each be a smaller size file than the gVCF file from which they are generated. The cohort files 311 and the census files 313 may each include a subset of data from the respective batch of gVCF files 304A, 304B, 304C, 304D from which they are generated. The cohort files 312 and the census files 314 may each include a subset of data from the respective batch of gVCF files 302A, 302B, 302C, 302D from which they are generated. For example, the cohort files 311, 312 may each include at least the AD, GQ, FT, QL, PL, and RAGT fields from the batch of gVCF files from which they are generated. The iterative gVCF genotyper may be configured to store single sample or multi-sample level variant data in the cohort files 311, 312. For example, an iterative gVCF genotyper can generate cohort files 311, 312 for each batch of samples (e.g., batch 1, batch 2, etc.) by aggregating data from the gVCF files associated with each of the samples in the respective batch for a subset of fields in the cohort files 311, 312.
[0052] The cohort files 311, 312 may be generated in a compressed data format used to store gVCF data for multiple samples. An exemplary cohort file for a batch having 12 samples may be represented by:
[0053] [Table 1]
[0054] As shown in the above example, each field in the cohort file may include a summary of metrics or values corresponding to similar fields in each of the aggregated gVCF files (e.g., gVCF files 302A, 302B, 302C, 302D for cohort 312, and gVCF files 304A, 304B, 304C, 304D for cohort 311, etc.). For example, as shown in the "DP" field, an iterative gVCF genotyper may identify each of the unique metrics in the "DP" field for each gVCF file in the batch. Each of the unique metrics or values in the "DP" field in the above example may include metric "32", metric "30", metric "42", metric "38", metric "41", metric "35", metric "37", and metric "21". Each metric may be followed by an array of integers including an integer identifier of the gVCF file, or a sample including the metric in the corresponding field in the gVCF file. For example, gVCF files with identifiers "1" and "7" contain a metric of "30" in the "DP" field, a gVCF file with identifier "2" contains a metric of "42" in the "DP" field, etc. The metrics or values may be stored in a text format in the cohort file. The integer identifiers of the gVCF files or samples corresponding to each metric may be stored as integer arrays in the cohort file to conserve the size of the cohort file and allow for additional compression of the file, as described herein. The cohort files 311, 312 may be used to calculate and store the number of variants and / or metric counts of which positions in the genome are hom-ref or no data.
[0055] The gVCF file format includes redundancy information. In each gVCF file of a sample, the same value of a metric in each field (e.g., FILTER, GT or ALT, NON REF) may be repeated in string format in different records. When aggregated as a whole, a gVCF file (e.g., gVCF files 302A, 302B, 302C, 302D, 304A, 304B, 304C, 304D, etc.) may include multiple repeated values for various metrics. For example, in the example provided above, a gVCF file with identifier "1" and a gVCF file with identifier "7" should each be stored with a text field containing the field value "DP" and a metric of "30". The above example for formatting a cohort file may allow for avoiding redundancy. For example, the formatting of a cohort file may allow for storage of summary information from a batch of gVCF files for a predefined subset of fields in a gVCF file efficient manner. As the number of samples sequenced continues to grow, the formatting of the cohort files described herein can continue to provide savings in storage and processing resources when analyzing sequencing data for genomic variant data and other genotype data. Also, the cohort files 311, 312 are more stable than gVCF files because the fields are fixed. Fixed fields in a gVCF file can be a subset of fields that contain data in the gVCF file that may be rarely changed by the user. The cohort files 311, 312 contain less repetition of field identifiers and values across samples, and therefore are less redundant than the storage of each of the gVCF files in a batch. Instead, one particular area of each of the cohort files 311, 312 can represent data across each of the samples in a batch of gVCF files for efficient storage.In the cohort files 311, 312, the data may be ordered by unique metric values, and for each unique value, a sample identifier value may be stored to avoid repeating the same value across multiple records in multiple samples.
[0056] The iterative gVCF genotyper can further analyze the data in each of the cohort files 311, 312 to generate a corresponding census file 313, 314. The census files may have a smaller file size than the cohort files 311, 312 from which the census files 313, 314 are generated. The census files 313, 314 may be used to store summary statistics for the variants (e.g., each of the variants) and / or reference blocks between samples within the cohort (e.g., batch) data represented in the corresponding cohort files 311, 312. An exemplary census file for a batch having 12 samples may be represented by the following:
[0057] [Table 2]
[0058] When a large number of samples are available, the user can split the samples into multiple batches, each with a similar sample size (e.g., 1000 samples). As shown in the example above, the census file 313, 314 can include a sample count for each unique metric or value in each field in the census file. The sample count can be calculated and stored in the census file in place of the sample identifier for each metric, as shown in the example cohort file. As shown in the example above, each field in the census file can include a summary of the metrics that correspond to similar fields in the corresponding cohort file. For example, as shown in the "DP" field, the iterative gVCF genotyper can calculate the total number of unique gVCF identifiers for each metric in the "DP" field of the cohort file. Each of the unique metrics in the "DP" field in the example census file above can include metric "32", metric "30", metric "42", metric "38", metric "41", metric "35", metric "37", and metric "21". Each metric may be followed by an integer value of the total number of unique gVCF identifiers that have the metric in the corresponding field. For example, the total number of gVCF files in a batch with a metric of "32" is "1", the total number of gVCF files in a batch with a metric of "30" is "2", etc. The metrics may be stored in a text format. A summary of the identifiers of the gVCF files corresponding to each metric may be stored as an integer value to conserve the size of the census file and allow for additional compression of the file, as described herein. The storage and processing of the sample counts in the census files 313, 314 may be smaller than the list of sample identifiers listed in the metric in the cohort files 311, 312, which may remove sample level information (or identifiers) from the census data in the census files that may be used to aggregate genomic variant data and other genotyping data at a higher level.Therefore, the census files 313, 314 may require less memory storage and processing resources to analyze than the cohort files 311, 312.
[0059] The cohort and census files may be in a file format that can be efficiently stored, analyzed, and / or communicated. In one example, the cohort and census files may be in a file format that can be efficiently compressed, indexed, and / or queried by genomic region using genomic and / or compression algorithms. The cohort file may store multiple fields extracted from the gVCF file for each genomic region (positions with variants or blocks without variants). As described further herein, for each field, the cohort file may store multiple values (e.g., metrics), and for each value (e.g., metric), an identifier for each sample that has a particular value of that field. The cohort file stores data per field across multiple samples, rather than per sample with multiple fields, and thus, querying one field across many samples is more efficient. The census file may store fields in a similar structure, but instead of storing an identifier for a sample for a particular value of one field, the census file may store the number of samples, which allows for the continuous accumulation of sample counts for different fields as batches of samples are aggregated. Storing sample counts instead of sample identifiers also reduces census file size, allowing for aggregation of multiple (e.g., thousands) batches of samples. In both cohort and census files, the data may be further encoded using bit encoding or hash encoding to reduce the data footprint, and the encoded binary data may be further compressed into compressed binary blocks using gzip or other public compression algorithms before being stored to disk. The compressed blocks may be indexed to allow random access of genomic regions. The ability to efficiently compress, index, and query files by genomic region may be a feature that can make variant data accessible when aggregating at scale.The use of these compressed files can reduce the disk footprint of the hardware and can reduce local and / or cloud workflows since data transfer can incur significant inbound and outbound costs as well as storage costs. Given the data structures of the files and described herein, the use of these compressed files can allow for more efficient transfer of data and lower costs for remote storage.
[0060] When census files 313, 314 are generated for multiple batches, the iterative gVCF genotyper can aggregate 320 the census files 313, 314 into a global census file 322. The census files 313, 314 can be aggregated from multiple batches (e.g., batch 1 and batch 2) from sequencing devices at different sites to generate the global census file 322. Although the example shown in FIG. 3 has two batches, it will be understood that 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 global census file 322 can include summary statistics for variant sites and non-variant genomic blocks across a single sample and / or multiple samples (e.g., multiple samples on a single sequencer). Each variant record in the global census file 322 can include RAGT statistics and variant statistics from the global census file, RAGT statistics and variant statistics from the batch census files 313, 314 (e.g., for samples included in the msVCF files for each batch), and other metrics (e.g., non-RAGT metrics) from other fields, such as FT, GQ, AD, PL (with sample identifiers) from the batch cohort files 311, 312. In one example, the global census file 322 can include a RAGT field that includes RAGT statistics for different alleles identified from different sample data received from batches of samples at different sites. The RAGT field can be a summary field used to determine each of the variant statistics (e.g., POS, REF, ALT, INFO fields) in summary fields in other files (e.g., in the msVCF files described herein). From the RAGT statistics, site statistics can be generated that are stored in census files (eg, per-batch census files 313, 314 and global census file 322).Site statistics can include the number of samples with genotype, without genotype (but with coverage), and without coverage, the total number of samples, alternative allele counts, and allele order, and / or a mapping between global alternative alleles and ALT alleles per sample as stored in the RAGT.
[0061] Processing of data received in gVCF files (e.g., gVCF files 302A, 302B, 302C, 302D, 304A, 304B, 304C, 304D) may be performed in parallel across multiple samples to generate a global census file 322 across multiple samples. For example, using parallel processing as described herein, cohort file 312 and census file 314 may be generated from batch 1 of gVCF files 302A, 302B, 302C, 302D, and cohort file 311 and census file 313 may be generated from batch 2 of gVCF files 304A, 304B, 304C, 304D. The global census file 322 may include one or more secondary metrics, such as, for example, allele counts, percentage of samples without sequencing coverage, percentage of samples without high confidence genotypes, etc. Generating the global census file 322 can be scaled to aggregate many (e.g., thousands) of batches more efficiently than aggregating gVCF files from a batch (e.g., all of them). For example, instead of having to open each gVCF file and search the file to find a piece of information (e.g., summary statistics or secondary metrics), the global census file 322 can aggregate the summary statistics and secondary metrics in one place.
[0062] When the global census file 322 is generated, the iterative gVCF genotyper may generate 330, 335 respective multi-sample VCF (msVCF) files 331, 332 for each batch. The msVCF files 331, 332 may be generated using the global census file, one of the cohort files 311, 312, and / or one of the census files 313, 314. For example, the msVCF file 332 may be generated using the cohort file 312, the census file 314, and the global census file 322 (e.g., for batch 1). The msVCF file 332 may be generated using the cohort file 311, the census file 313, and the global census file 322 (e.g., for batch 2). The process 300 may end after the generation of the msVCF files 331, 332. The msVCF files 331, 332 can include batch-specific genomic variant data and / or other genotype data, as well as global census data including genomic variant data and / or other genotype data identified from samples acquired at other sites. Examples are the global set of variant sites and variant alleles, genotypes called per sample, read depth for variant alleles, and hom-ref positions. Likelihood and quality scores for alleles.
[0063] In some examples, the msVCF file may include one or more global statistics, such as allele frequency, read depth and / or home-ref location for variant alleles, number of samples with or without genotype, and number of samples without coverage. The msVCF file may include likelihood and / or quality scores for alleles. Similar statistics between samples within a batch are also included. The msVCF file may be generated for each batch of samples to facilitate data read / write efficiency, data transfer, storage, and / or querying. The number of records in each msVCF file may be the same across all batches, allowing for quick search for the same variant across multiple batches.
[0064] Fields in each of the msVCF files 331, 332 may include summaries of genotypes and / or samples. For example, fields in the msVCF files 331, 332 may include allele counts in genotypes, total number of alleles in called genotypes, total number of samples, total number of samples with called genotypes, total number of samples with unknown genotypes, and / or total number of samples with no coverage. The metrics calculated for each of these fields may be generated in multiple versions. For example, the metrics for each of these fields may be calculated once per batch to provide a local version of the metric and once for each of the samples in the cohort file to provide a global version of the metric.
[0065] As will be appreciated, the DNA sequencing process may have a random component. A sequencing run may not yield any information (=coverage from sequencing reads) for a small set of genomic locations. This may be captured by a "no coverage" metric. The number of these locations may be used to distinguish them from other locations. For example, variant information or statistics may include variant data indicating that the location has a DNA mutation, home-ref data indicating that the location does not have a DNA mutation, or no coverage data indicating that the variant information of the location is unknown. A gVCF variant file may include this information for a particular sample. A gVCF file may include, for each location in the genome, variant data indicating that the location has a DNA mutation, home-ref data indicating that the location does not have a DNA mutation, or no coverage data indicating that the variant information of the location is unknown. The census files 313, 314 can store this information for batches of samples, merging the census files for each batch into a global census file 322, and then sending the global census back to each batch during msVCF writing to include these metrics in the msVCF file for each batch.
[0066] msVCF files may be a standardized, public data format that allows users to perform downstream analyses, such as genome-wide association studies (GWAS), imputation and phasing, gene burden analysis, rare variant discovery, population-specific allele frequencies, population substructure analysis, and clinical analyses to estimate pathogenicity of variants and classify pathogenic / benign variants. msVCF files may allow access to and query of aggregated variant data (e.g., within a batch). msVCF files may include variant site-level information both within the current batch and globally across all batches of samples from the sequencing device at each site. Variant site-level information may include allele count, allele frequency, total number of samples, number of samples with no coverage, number of samples with coverage but not confidently genotyped, and number of samples with genotypes. Users may prefer to use msVCF files to access and query aggregated variant data, for example, instead of other file formats. An exemplary msVCF file for a batch may be represented by:
[0067] [Table 3]
[0068] When a global census file (e.g., global census file 322) has already been generated and a batch of samples becomes available from the sequencer, the iterative gVCF genotyper can aggregate the sample data in the batch into a cohort file (e.g., cohort files 311, 312) and a census file (e.g., census files 313, 314) for the batch. The iterative gVCF genotyper can then aggregate the census file from the batch with the global census files from (e.g., all) previous batches to generate, for example, an updated global census file.
[0069] After the global census file has been updated (e.g., new variant sites have been discovered and / or variant statistics have been updated at existing variant sites), the iterative gVCF genotyper can again generate msVCF files for each batch of samples using the cohort file for the batch, the census file for the batch, and the updated global census file. The msVCF files can include the variants and alleles found in each of the samples from each of the batches received from the genotyping devices across the sites. The msVCF files 331, 332 can then include the updated global census data, including batch-specific genomic variant data and / or other genotype data, as well as genomic variant data and / or other genotype data identified from samples obtained at other sites. The iterative gVCF genotyper can provide an iterative population-based analysis option for analyzing samples from unrelated individuals together, e.g., as samples become available for analysis.
[0070] The gVCF file may be subject to block compression, with certain columns used to indicate which chromosome the data belongs to, and one or more columns indicating start / end points. The process 300 may incorporate a hybrid of binary and ASCI or string compression and / or serialization to store the variant data and minimize the amount of storage space and / or processing resources required to generate the msVCF file.
[0071] An iterative gVCF genotyper can use the exemplary process 300 to aggregate sample data from multiple sequencing devices at multiple sites into an existing data file, such as an msVCF file. An iterative gVCF genotyper can use the exemplary process 300 to, for example, incrementally aggregate newly available genomic variant data and other genotype data from a batch of sample data with genomic variant data and other genotype data available in a preceding batch, without having to redo the analysis on the previously available (e.g., and aggregated) batch. For example, an iterative gVCF genotyper can incrementally aggregate genomic variant data and / or other genotype data as it becomes available. An iterative gVCF genotyper can be scalable, e.g., run on multiple computing platforms, such as a cloud platform, a high performance cluster, and / or a single server.
[0072] FIG. 4 is a flow chart of an exemplary method 400 for processing a batch of sample data from a sequencing device to generate one or more msVCF files. The exemplary method 400 may be performed by one or more computing devices. For example, the exemplary method 400 may be performed by one or more processors executing computer-readable or machine-readable instructions stored in memory in one or more computing devices to operate an iterative gVCF genotyper. The exemplary method 400 may be performed in response to input received from a user to software for constructing one or more msVCF files. One or more portions of the exemplary method 400 may be performed locally by a sequencing application 110 executing on a client device 108 and / or remotely by a sequencing system 104 executing on one shown in FIG. 1. One or more portions of the exemplary method 400 may be performed remotely, such as by a sequencing system 104 on a server device 102 shown in FIG. 1, in a cloud-based environment or other remote computing device. Although the exemplary method 400 may be described herein as being performed by a single device, the exemplary method 400 may be distributed across multiple devices. Similarly, although the exemplary method 400 may be described herein as being performed by a single processing entity, the exemplary method 400 may be performed using parallel processing of batches of sample data, as described herein. For example, one or more portions of the exemplary method 400 may be implemented by a virtual machine or compute node having separate resources for performing one or more portions of the exemplary method 400.
[0073] As shown in FIG. 4, method 400 may begin at 402 when a batch of sample data is received in one or more gVCF files. The batch of sample data may be received from a sequencing device capable of processing samples to generate sample data, or another device capable of generating and / or transmitting a gVCF file. The sample data may include genomic variant data and / or other genotype data. For example, the sample data may include data in one or more fields of a gVCF file as described herein. At 404, a processor of the computing device may aggregate cohort data for the batch of samples from the sample data in the gVCF file. The cohort data may include sample data in a subset of the fields in the gVCF file. The processor of the computing device may aggregate census data for the batch of sample data at 406. The census data may include variant summary statistics and reference blocks of the sample data in the batch. The cohort data and census data may be stored in separate file types, databases, or other data structures.
[0074] The processor of the computing device may determine 408 whether there is an existing global census file for aggregating global census data corresponding to sample data from each of the batches being received from the sequencing devices at the multiple sites. If no global census file exists or does not exist for aggregating sample data in the received batches, the processor of the computing device may generate 414 a global census file including census data for the batch of samples. The global census file may include variant summary statistics that are global metrics (e.g., allele counts, allele frequencies, etc.) for the sample data from each of the batches of samples received. If a global census file already exists, the processor of the computing device may update the global census file at 410 by aggregating the census data for the batch with the existing data in the global census file. This may include adding the census data for the batch received to the sample data and / or generating variant summary statistics using the sample data in the batches and adding the variant summary statistics to the global census file.
[0075] At 416, the processor of the computing device can generate an msVCF file for the batch being processed. The msVCF file can be generated for the batch using the cohort data, the census data for the batch, and / or the global census file. Each variant record in the global census file can include RAGT and variant statistics from the global census file, RAGT and variant statistics from the batch census file (e.g., for the samples included in the msVCF file for each batch), and other metrics (e.g., non-RAGT metrics) from other fields, such as FT, GQ, AD, PL (with sample identifier) from the batch cohort file. The RAGT and variant statistics can include variant position, variant allele, sample genotype and / or quality of the variant call, and / or other RAGT and variant statistics. The computing device can store the msVCF file (e.g., locally, remotely, and / or in the cloud). After the msVCF file is generated at 416, the process 400 may determine at 418 whether there are additional batches of sample data available from a sequencer or other computing device. If there are additional batches to be processed, the procedure may return to 404 to aggregate the cohort and census data for the batch and proceed as described herein. If there are no additional batches of sample data available from the sequencer, the method 400 may end. The cohort files, census files, global census files, gVCF files, and / or msVCF files may be stored in memory in a manner that improves or optimizes memory usage. For example, the data may be stored in memory using hash codes and / or hash maps, as further described herein.
[0076] As shown in the method 400 depicted in FIG. 4 and described otherwise herein, a global census file may be generated from multiple batches of sample data received from sequencing devices at different sites, and may be updated as additional batches of sample data are received from sequencing devices at additional sites. FIGS. 5A and 5B show exemplary processes 500, 550 illustrating how a global census file may be generated and updated, and how msVCF files may be generated and / or updated in response to sample data in an updated global census file. FIG. 5A shows an exemplary scenario when there are two batches of gVCF files simultaneously available for processing. FIG. 5B shows a further example when an additional batch of gVCF files is later available. The exemplary process 500 depicted in FIG. 5A may be performed by an iterative gVCF genotyper to initially generate a global census file. The exemplary process 550 depicted in FIG. 5B may be performed by an iterative gVCF genotyper to update the global census file as additional batches of sample data are received. The processes 500, 550 may be performed by one or more computing devices. For example, the processes 500, 550 may be performed on multiple computing devices, each of which performs an independent region or chromosome in the genome. The processes 500, 550 may be performed by one or more processors executing computer-readable or machine-readable instructions stored in memory in one or more computing devices. The processes 500, 550 may be performed in response to input received from a user in software for constructing one or more msVCF files. One or more portions of the processes 500, 550 may be performed locally by a sequencing application 110 executing on a client device 108 shown in FIG. 1. One or more portions of the processes 500, 550 may be performed remotely, such as by a sequencing system 104 on a server device 102 shown in FIG. 1, in a cloud-based environment or other remote computing device.Although the processes 500, 550 may be described herein as being performed by a single device, the processes 500, 550 may be distributed across multiple devices. Similarly, although the processes 500, 550 may be described herein as being performed by a single processing entity, the processes 500, 550 may be performed using parallel processing of batches of sample data, as described herein. For example, one or more portions of the processes 500, 550 may be implemented by a virtual machine or compute node having separate resources for performing one or more portions of the processes 500, 550.
[0077] As shown in FIG. 5A, a processor of a computing device running an iterative gVCF genotyper may receive two or more batches of sample data at 510. For example, the two or more batches of sample data may be available from a sequencing device capable of processing samples to generate sample data or another device capable of generating and / or transmitting gVCF files. Each batch of sample data may include data from a plurality of samples (e.g., 1000 or more samples). For example, the processor may receive at 510 a gVCF file associated with a first batch of sequencing data and a gVCF file associated with a second batch of sequencing data. The gVCF files associated with the first batch of sequencing data and / or the second batch of sequencing data may be split into shards of equal size as described herein. Each of the shards may be processed using one of a plurality of computation nodes.
[0078] A processor running the iterative gVCF genotyper may aggregate (e.g., in parallel) the sample data within each batch to generate a cohort file and a census file for each of the batches at 510. For example, a first cohort file for a first batch (e.g., batch 1) may include a subset of the fields in each of the gVCF files associated with the first batch (e.g., batch 1) of sequencing data. The census file for the first batch of sequencing data may include variant summary statistics (e.g., RAGT statistics) and / or hom-ref blocks for the first batch of sequencing data. A second cohort file for a second batch (e.g., batch 2) of sequencing data may include a subset of the fields in each of the gVCF files associated with the second batch of sequencing data. The census file for the second batch of sequencing data may include variant summary statistics (e.g., RAGT statistics) and / or hom-ref blocks for the second batch of sequencing data. It should be appreciated that each of the cohort files generated for each batch may include the same subset of fields.
[0079] A processor running the iterative gVCF genotyper can aggregate census files from each of the batches to generate a global census file at 520. The global census file can include census data from batches of samples received from sequencing devices at different sites. A processor running the iterative gVCF genotyper can generate msVCF files for the first and second batches of sequencing data at 530. For each of the two or more batches, the iterative gVCF genotyper can use the cohort file for each batch, the census file for each batch, and the global census file to generate an msVCF file for each batch (e.g., in parallel). One or more of the msVCF files can be used to perform genome-wide sequencing analysis. Although process 500 shows an initial global census file generated at 520 using variant summary statistics from two batches (e.g., batch 1 and batch 2), it will be understood that process 500 is not limited to generating an initial global census file using two batches of sample data. Alternatively, process 500 may use multiple batches (e.g., three or more batches) of sample data that are available. Each of the batches of sample data may be processed in parallel and / or used to generate an initial global census file at 520. It should also be appreciated that the census file may be updated by aggregating census files of subsequent batches of sample data with a previously generated global census file.
[0080] 5B, the global census file may be updated in response to additional sample data received in an additional batch of sample data. For example, process 550 may begin after process 500 has ended. Process 550 may receive an additional batch of sample data (e.g., batch 3) after the first two msVCF files have been generated, for example.
[0081] The processor running the iterative gVCF genotyper may aggregate sample data in additional batches (e.g., batch 3) at 560 to generate cohort files and census files for the additional batches. The iterative gVCF genotyper may aggregate data across batches of samples in a global census file at 570 such that the global census file is updated based on sample data in the additional batches of samples to generate an updated global census file. For example, the processor running the iterative gVCF genotyper may generate an updated global census file by incorporating data in the census file from the additional batches into the global census generated in process 500. Being able to update the global census file and not having to re-aggregate census data from previous census files may conserve processing resources (e.g., for generating a global census file with census data from all received batches and / or for generating msVCF files for each batch).
[0082] The processor running the iterative gVCF genotyper can use the cohort files for the additional batches, the census files for the additional batches, and the updated global census file to generate msVCF files for the additional batches, at 580. The iterative gVCF genotyper can also use the cohort files for the existing batches, the census files for the existing batches, and the updated global census file to generate updated msVCF files for each of the previously received batches, at 580. For example, the processor running the iterative gVCF genotyper can generate a new msVCF file for each batch upon generation of the updated global census file (e.g., and upon receipt of additional batches of sequencing data).
[0083] In Figures 5A and 5B, the batches of gVCF files may be processed by two levels of parallelization by genomic region and sample batch. These two levels may be independent of each other. For example, the analysis of each genomic region may be performed in parallel. For each genomic region, different batches may be further parallelized to generate cohort and census files. Once the cohort and census files are generated for each batch, a processor running the iterative gVCF genotyper may aggregate the census files from each of the batches into a global census file. The generation of the global census file may be parallelized by region. Figures 5A and 5B show each batch containing 1000 samples, but it will be understood that each batch may contain more or less than 1000 samples. Although Figure 5B shows a total of three batches, it will be understood that the iterative gVCF genotyper may receive gVCF files for any number of batches.
[0084] FIG. 6 shows a graphical representation 600 of cohort records generated for a batch to show how similar types of data in multiple gVCF files can be aggregated into a cohort file. Similar types of information can also or alternatively be aggregated into a census file. As shown in FIG. 6, a set of input gVCF files 610 can be received by a processor running an iterative gVCF genotyper. The batch in the example shown in the graphical representation 600 includes five samples. As shown in FIG. 6, the gVCF files 610 can include related samples, such as a trio of three related samples (e.g., father, mother, and child). Each sample can be associated with a respective gVCF file. Each of the gVCF files 610 includes hom-ref data, variant data, and does not include coverage data distributed into sections throughout the gVCF file. For example, each of the gVCF files 610 indicates hom-ref data, one or more variants, and does not indicate a coverage area for its respective sample. The hom-ref data, variant data, and no coverage data from the samples in a batch (e.g., a set of input gVCF files) may be aggregated into a cohort file 620.
[0085] The cohort file 620 may be a single file (e.g., a set of input gVCF files) showing hom-ref, variant, and no coverage data for a batch. In the cohort file 620, data of similar types (e.g., hom-ref, variant, and no coverage data) may be aggregated into rows of the cohort file. The cohort file 620 may include blocks (e.g., all hom-ref blocks or all no coverage) that include regions where each of the samples in the cohort has the same type of record in the original gVCF. Although not shown, census records may be generated for the batch and the census records may be similarly aggregated. While FIG. 6 shows an input gVCF file 610 for a batch containing five samples, it will be understood that more than five input gVCF files 610 may be received. For example, the number of input gVCF files 610 received may be equal to the number of samples in the batch.
[0086] Records in a cohort file can be grouped into regions with the same type of records (e.g., hom-ref, no coverage, or variant) to reduce data size and / or reduce processing resources for compression. In contrast, the cohort file for each region can be decomposed into one record per genomic location, which has adjacent records with the same value across samples.
[0087] FIG. 7 illustrates an exemplary iterative process 700 that may be implemented by one or more computing devices (such as, for example, the server device 102, the client device 108, and / or the sequencing device 114 shown in FIG. 1) to generate a multi-sample variant call file upon receipt of a gVCF file associated with a batch. The one or more computing devices performing the exemplary iterative process 700 may be running an iterative gVCF genotyper 760 as described herein. The iterative gVCF genotyper 760 may use the iterative process 700 to aggregate newly available batches / sample data into existing cohorts. The iterative gVCF genotyper 760 may use the iterative process 700 to, for example, incrementally aggregate newly available batches of sample data with previously available batches without having to redo the analysis for previously available (e.g., and aggregated) batches. The iterative process 700, or portions thereof, may be implemented to more efficiently store variant data in census files, cohort files, and / or msVCF files. For example, the iterative process 700, or portions thereof, may be performed to aggregate gVCF files into a cohort file and census files into a global census file.
[0088] The iterative process 700, or portions thereof, may be performed on a single computing device or may be distributed across multiple computing devices (e.g., multiple servers, sequencing devices, and / or client computing devices). The process 700, or portions thereof, may reduce the amount of processing resources and / or storage space (e.g., memory) used by a computing device during variant analysis. The iterative process 700 may be used to process multiple batches of sequencing runs. The iterative process 700 may store the cohort file 720, the census file 730, the global census file 740, and / or the msVCF file 750 locally or remotely.
[0089] The iterative process 700 may begin upon receipt of a gVCF file 710 associated with a batch of sample data by an iterative gVCF genotyper 760. The iterative gVCF genotyper 760 may read the gVCF file 710 and identify a subset of VCF fields in the gVCF file 710 to generate a cohort file 725. The iterative gVCF genotyper 760 may aggregate the gVCF data read from each of the gVCF files in the batch of gVCF files 710 in the batch. The gVCF file may include sample data from samples of a sequencing run. The cohort data 710 may include a subset of the sample data in the gVCF file 710 that is identified in a predetermined field. The cohort data 720 may be aggregated from each of the predetermined fields in the batch of sample data read from the gVCF file 710, and the iterative gVCF genotyper 760 may write the subset of data to the cohort file 725. As described herein, the cohort data may include a summary of unique identifiers for the gVCF files that contain a metric or value for each of the predefined fields in a batch of gVCF files 710.
[0090] The iterative gVCF genotyper 760 can convert the cohort data 720 into census data 730 and write the census data to a census file 730 for the batch of sample data. The census data 730 can include a count of the number of unique gVCF files that contain a common metric or value for each of the predetermined fields in the cohort file 725. The iterative gVCF genotyper 760 can aggregate the census data 730 from the multiple census files 740 to generate a global census file. For each of the multiple batches, the iterative gVCF genotyper 760 can generate and write an msVCF file 750 for the respective batch of samples using the cohort data 720 and / or census data 730 for the respective batch of samples and the data in the global census file, as described herein.
[0091] FIG. 8 illustrates an exemplary process 800 that may be implemented by one or more computing devices (e.g., multiple server devices 102, multiple client devices 108, and / or one or more server devices and one or more client devices 108, etc., shown in FIG. 1) that may utilize parallel processing in variant analysis of sample data received in gVCF files. For example, the exemplary process 800 may incorporate parallel processing for processing each of batches 802, 804, 806 of sample data from sequencing devices at different sites to generate one or more msVCF files. The process 800, or portions thereof, may be performed by an iterative gVCF genotyper distributed across multiple computational nodes to more efficiently store and process variant data in msVCF files. The process 800, or portions thereof, may reduce the amount of processing resources and / or storage space (e.g., memory) used by a computing device during variant analysis.
[0092] The iterative gVCF genotyper can utilize at least one processor and / or at least one memory to perform parallel processing using multiple compute nodes or virtual machines. The iterative gVCF genotyper can allocate processing resources as described herein to enable the parallel processing described. The iterative gVCF genotyper can divide each sample (e.g., each gVCF file) in a batch (e.g., such as batch 802, batch 804, and / or batch 806) into separate shards. As each sample is assigned to a shard, the samples can be freely grouped into defined groups (e.g., subpopulations or case / control groups). Each shard can be of equal or different size. Each shard can be processed by a separate compute node. The iterative gVCF genotyper may attempt to create shards of equal size, but the shards may be of unequal size. For example, a shard may be blocked from spanning a chromosome, which may result in shards ending up on chromosomes of unequal shard size. Each shard may be of equal size, except for the ends of the chromosomes (eg, remnants), which may be less than equal in size.
[0093] Each core processor accessible by a compute node or virtual machine may have a particular threading. Each core may have a particular threading. A variable number of software threads may be implemented. One or more threads may be implemented per CPU core. For example, a single thread may be implemented by each CPU core. The number of threads implemented by each CPU core may be changed in response to user input. Each thread of the iterative gVCF genotyper may be assigned to extract sequencing data from a respective gVCF file or portion of a gVCF file to create a respective cohort and census file for the respective gVCF file or portion of a gVCF file at 810. The iterative gVCF genotyper may split the cohort and / or census file generation of each of the batches 802, 804, 806 into respective shards such that each of the shard-generated cohort and census files is generated by a separate shard. Each of the census files generated by the shards processing a common batch of samples may be aggregated into a batch census file at 820. A batch census file may contain similar fields to the global census file, including region, RAGT statistics, and optionally site statistics. A batch census file may contain census data for a given batch of gVCF files, whereas a global census file may contain census data from multiple batches of samples, and thus the batch census file may differ from the global census file in sample size.
[0094] The census files created from shards processing gVCF files in batch 802 may be aggregated to create a first batch census file for batch 802, the census files created from shards processing gVCF files in batch 804 may be aggregated to create a second batch census file for batch 804, and the census files created from shards processing gVCF files in batch 806 may be aggregated to create a third batch census file for batch 806. Each of the batch census files may be aggregated to generate a global census file.
[0095] Batch census file RAGT statistics may be calculated from the batch cohort files (e.g., by taking sample counts per metric value). The global census file may contain RAGT statistics aggregated from the RAGT statistics in the census files per batch by summing the sample counts per metric value. The global census file may contain variant metrics and hom-ref metrics for each of the samples in the cohort files. The global census file may also contain this information per genomic location. The batch census file and the global census file may contain site statistics fields that contain variant calling data, e.g., normalization of reference and alternative alleles, ordering and remapping of alternative alleles in indexes between common alternative alleles, and alternative alleles in each RAGT case.
[0096] Since each shard can be implemented to generate cohort data for a respective gVCF file or portion of a gVCF file, the shards can each use their respective cohort files in the batch to generate a portion of the multi-sample VCF file. For example, each shard implemented to process a respective gVCF file or portion of a gVCF file in batch 802 to generate a corresponding cohort file can use their respective cohort file and the global census file to create a portion of msVCF842. Each shard implemented to process a respective gVCF file or portion of a gVCF file in batch 804 to generate a corresponding cohort file can use their respective cohort file and the global census file to create a portion of msVCF844. Each shard implemented to process a respective gVCF file or portion of a gVCF file in batch 806 to generate a corresponding cohort file can use their respective cohort file and the global census file to create a portion of msVCF846. Thus, as shown at 830 in process 800, each msVCF 842, 844, 846 may be generated on a per-shard, per-batch basis. When each sample is assigned to a shard, the samples may be freely grouped into defined groups (e.g., subpopulations or case / control groups, etc.). When each sample is assigned to a shard, the shard may handle (e.g., only handle) one sample, which may allow for efficient compression compared to a cohort of multiple samples.
[0097] Each compute node or virtual machine can implement at least two levels of parallelism to process, aggregate, and / or generate data for a corresponding region of the sequence data. Each compute node or virtual machine can process a particular region of data in the gVCF file within a given batch.
[0098] The iterative gVCF genotyper can include runtime option parameters for implementing process 800. The runtime option parameters can set the number of regions to split the genome into to parallelize the process on distributed nodes. Because each core may have a specific threading, the runtime option parameters can also set the number of sub-regions to split a region into to run on different threads on the same node. The runtime option parameters can set data buffer sizes (e.g., for further splitting sub-regions) to ensure that each thread does not utilize an excessive amount of memory at any one time.
[0099] FIG. 9 illustrates an exemplary process 900 incorporating parallel processing across batches of sample data. In an example, an iterative gVCF genotyper can allocate processing resources to enable the parallel processing described. The diagram illustrated in FIG. 9 illustrates each batch being processed by multiple shards. In the exemplary process 900, each shard can be assigned a particular portion of the cohort and census files for each batch. For example, a shard (e.g., shard 1) can process its portion of the cohort and census files for batch 1, and then proceed to its portion of the cohort and census files for batch 2 once processing of its portion of the cohort and census files for batch 1 is complete. The shard can then proceed to its portion of the cohort and census files for batch 3 once processing of its portion of the cohort and census files for batch 2 is complete. Each shard can have a uniform definition, such that the same genome shards from different batches can be aggregated (e.g., without discrepancies in shard boundaries). Each shard can utilize multiple threads to process a portion of each batch. Within each shard, the thread-level partitioning (or sub-regions) may be adjusted differently between batches because the output files (e.g., cohorts or censuses) of a multi-threaded process are (e.g., always) concatenated into one output file (e.g., cohort or census) per shard.
[0100] In the exemplary process 900, each shard may be assigned a particular portion of the global census file to be generated by the shard. For example, a shard (e.g., shard 1) may process its portion of the global census file as the batch cohort and its portion of the census file are completed. The ends of chromosomes may be used to create the ends of one or more shards. For example, one shard may be assigned to mitochondrial (MT) genomic data and / or one shard may be assigned to alternate configurations and Human Leukocyte Antigen (HLA). These additional shards (e.g., alternate contigs including mitochondrial and HLA) may be present in the human reference genome. They may be kept in separate shards because the ploidy may differ from the autosomes and sex chromosomes. The size of these contigs may be relatively small and may not require implementing further parallelization with more shards.
[0101] In the example process 900, each shard can be assigned a particular portion of each msVCF file. For example, a shard (e.g., shard 1) can process its portion of the msVCF file for batch 1 and then proceed to its portion of the msVCF file for batch 2 once processing of its portion of the msVCF file for batch 1 is complete. The shard can generate a shard-specific portion of the msVCF file for each batch using its portion of the batch file and its portion of the global census file that correspond to the msVCF file being generated. The shard can then proceed to its portion of the msVCF file for batch 3 once processing of its portion of the msVCF file for batch 2 is complete.
[0102] 10 illustrates an exemplary genomic data operator design 1000. The exemplary genomic data operator design 1000 may be used to read, process, and / or write genomic data. Region and / or subregion sizes may be configurable depending, for example, on the ram size per thread and / or the type and number of input files and in-memory data structure size. The exemplary genomic operator design 1000 may be used to process a batch of sample data from a sequencing device to generate one or more msVCF files.
[0103] FIG. 10 illustrates parallelization that may be performed on the same compute node (e.g., when the shards are already pinned). Given a region as defined by one genome shard and the number of threads on a given compute node, the analysis may be further parallelized by concurrent threads by dividing the shard into regions per thread. As each thread processes one such region, a portion of the data in the region may be loaded from the file system. For example, each thread may process records in one subregion at a time. The subregion size may be configured to maximize memory usage per thread and / or avoid overloading the total memory of the system when each of the threads is fully loaded. Each thread may process the subregions contained in the region in serial. The threads may process the regions contained in the shard in parallel.
[0104] In generating the global census files and msVCF files as described herein, the global census files and msVCF files themselves may grow in size as global census data is aggregated from different batches of sample data, updated in the global census files, and / or included in the msVCF files. The global census data may store therein representations of alleles represented in fields of the gVCF files where data is aggregated from multiple sites (e.g., RAGT fields). The size of the global census data in the global census files and / or msVCF files may be reduced by reducing overlapping alleles that may be represented differently in the global census data. Reducing the size of these files may allow for reduced memory and processing requirements when storing and analyzing the global census data.
[0105] One method to reduce the representation of overlapping alleles would be to normalize each allele based on the reference allele. This normalization may also allow for more accurate allele counts for possible variance that may be found for a cohort. FIG. 11 shows an example process 1100 that may be implemented by one or more computing devices (such as, for example, the server device 102, the client device 108, and / or the sequencing device 114 shown in FIG. 1) to normalize and reindex genotypes for msVCF files. Although the process 1100 shows a normalization and reindexing technique that may be implemented on msVCF files, similar normalization and reindexing techniques may be implemented for global census files. The example process 1100, or portions thereof, may use RAGT statistics to normalize population-level reference alleles and alternative alleles, order population-level alternative alleles, reindex sample-level genotypes, and store the data in census files and msVCF files. The reference and alternate alleles and GT re-index and ALT index mapping generation in FIG. 11 can be performed during census aggregation and / or gVCF aggregation resulting in site statistics. The site statistics can be stored in a site statistics field of the census file. The site statistics can be inserted (e.g., in a site statistics field) into the msVCF file during generation. Process 1100, or portions thereof, can be performed to efficiently store genomic variant data in the msVCF file and / or the global census file. Process 1100 can be implemented by an iterative gVCF genotyper running on one or more processors of one or more computing devices. Process 1100, or portions thereof, can be performed on a single computing device or distributed across multiple computing devices (e.g., multiple servers, sequencing devices, and / or client computing devices).Process 1100, or portions thereof, may reduce the amount of processing resources and / or storage space (e.g., memory) used by a computing device during variant analysis.
[0106] The exemplary process 1100 may begin when a processor running an iterative gVCF genotyper receives RAGT statistics in gVCF files at 1102. One or more of the gVCF files may have been recently received by the iterative gVCF genotyper. Each of the gVCF files may be associated with a respective sample of the plurality of samples. Additionally or alternatively, one or more of the gVCF files may have been previously received and stored in one or more computing devices and / or databases (such as, for example, database 116 shown in FIG. 1). The RAGT statistics may be associated with the plurality of samples. For example, the RAGT statistics may be stored in one or more computing devices and / or databases (such as, for example, database 116 shown in FIG. 1), and the process 1100 may be used to update genotype re-indexing, reference alleles, and / or alternative allele mappings based on the received (e.g., new) RAGT statistics. A processor running an iterative gVCF genotyper can use RAGT statistics to identify a plurality of reference alleles and a plurality of alternative alleles associated with a plurality of samples.
[0107] The iterative gVCF genotyper can sum 1104 the instances (e.g., occurrences) of each unique allele (e.g., multiple reference alleles and multiple alternative alleles) into an allele count, e.g., based on a RAGT statistic. For example, the one or more computing devices can determine how many instances of each unique allele (e.g., reference alleles and alternative alleles) are present in the RAGT statistic. Although the example shown in FIG. 11 includes eight unique alleles, it will be understood that the exemplary process 1100 is not limited to this number of alleles. Rather, the exemplary process 1100 can be used when the number of unique alleles is greater than or less than eight.
[0108] The iterative gVCF genotyper can order the allele counts, for example, based on length, at 1106. The allele order determined at 1106 can start with the reference allele count and end with the alternative allele count. Stated differently, the reference allele counts can be listed before the alternative allele counts in the allele order. For example, the longest reference allele can be listed first in the allele order, and the remaining reference alleles can be listed in order of decreasing length. Each of the ordered allele counts can be assigned a number between 0 and the total number of unique alleles in the RAGT statistic (e.g., 7 in the example shown in FIG. 11). Zero (0) can be assigned to the longest reference allele, and the remaining alleles can be ordered based on occurrence. For example, the alleles can be listed in descending order by occurrence in the RAGT statistic. It will be appreciated that the alleles are not limited to being ordered in descending order by occurrence, and instead, the alleles can be ordered in other ways, such as, for example, by descending order of length.
[0109] At 1108, the alleles may be normalized, for example, based on the longest reference allele. For example, the iterative gVCF genotyper may select a normalized reference allele. The normalized reference allele may be the longest reference allele. The alleles per sample in the RAGT statistic may be normalized by extending the alleles to the length of the longest reference allele. The alternative alleles of each reference allele may be extended by the same amount (e.g., number of bases) as the corresponding reference allele. For example, when a reference allele is extended by two (e.g., two base pairs or one nucleotide), each alternative allele of the reference allele is also extended by two (e.g., two bases). A normalized representation of each sample may be generated using the normalized reference allele such that each of the multiple alternative alleles is indexed using the normalized reference allele.
[0110] Normalization of alternative alleles may allow alternative alleles with the same reference to be placed in the same row of data. The same reference allele and alternative alleles may be grouped by the same reference allele. To allow normalized reference alleles and normalized alternative alleles to be included in the same row,
[0111] At 1110, the genotypes for each sample are reordered, for example, based on the normalized reference alleles and / or alternative alleles. For example, assume that the original genotype in the RAGT field is 1 / 2, the original reference allele is TACAC, and the original ALT allele is TAC,T. After renormalization, if the common reference allele is TACACACACAC (ACACAC added), the normalized alternative alleles are TACACACAC,TACACAC. According to the respective ordering of the normalized reference alleles and alternative alleles, if the TACACACAC allele has index 3 and TACACAC index 5, the ordered genotype is 3 / 5 (instead of 1 / 2), and the alternative allele mapping (old to new) is 1=>3,2=>5. After allele reordering, at 1112, the alternative allele mapping for each sample can then be generated, for example, based on the normalized reference alleles and / or alternative alleles. A common reference allele may refer to a common (renormalized) reference allele shared among each of the samples at a particular site. Because different samples may have different alternative alleles at the same genomic location, reordering may be performed at 1112. Consistent ordering across locations can be enforced for writing to the msVCF. This reordering of the genotype field may be implemented because the genotype may refer to an index in a list of alternative alleles.
[0112] The normalized reference and alternative alleles, reordered genotypes, and / or alternative mappings may be stored in a site statistics field of a census file that contains variant calling data, e.g., normalization of reference and alternative alleles, ordering and remapping of alternative alleles in indexes between common alternative alleles, and alternative alleles in each RAGT case. An exemplary site statistics field for a batch with 12 samples may be represented by:
[0113] [Table 4A] [Table 4B]
[0114] The site statistics field can contain RAGT statistics for the batch. The site statistics field can contain number counts for each census field value. The number counts can be used for population genomics. The site statistics field can contain the primary census file information used as input to the msVCF file.
[0115] The normalized reference and alleles, reordered genotypes, and alternative mappings can be used to generate msVCF files for multiple samples associated with RAGT statistics. For example, the normalized reference and alleles, reordered genotypes, and alternative mappings can be output to an msVCF file. The msVCF file can include the common reference alleles in a REF column. The generated msVCF files can be stored by the iterative gVCF genotyper on one or more computing devices and / or databases.
[0116] The process 1100 may be repeated when additional gVCF files are received. The RAGT statistics from the additional gVCF files may be added to the previously received RAGT statistics. For example, the iterative gVCF genotyper may receive additional gVCF files associated with one or more additional samples. The iterative gVCF genotyper may identify one or more reference alleles and one or more alternative alleles associated with the one or more additional samples. The one or more reference alleles and one or more alternative alleles may be added to the allele count. The one or more computing devices may determine whether any of the one or more reference alleles is longer than the longest reference allele. When a reference allele (e.g., in an additional gVCF file) is longer than the previous longest reference allele, the iterative gVCF genotyper may select the reference allele as the updated longest reference allele. The multiple reference alleles and multiple alternative alleles may be normalized by extending them to correspond to the length of the updated longest reference allele. The genotype ordering and alternative mapping may be updated based on the information received in the additional gVCF files. When this data may be updated, the global census file may be updated and a separate msVCF file may be generated for each batch of data received.
[0117] Exemplary embodiments are described herein for generating cohort and census files having smaller file sizes than the gVCF files containing the sample data from which the cohort and census files are generated. The size of the cohort and census files may be further reduced using one or more compression and / or serialization techniques described herein. FIG. 12A illustrates an exemplary process 1200 that may be implemented by one or more computing devices (e.g., such as the server device 102, the client device 108, and / or the sequencing device 114 shown in FIG. 1) to encode the cohort and / or census data in a compressed format. The iterative process 1200, or portions thereof, may be performed to efficiently store variant data in the census files, cohort files, and / or msVCF files. The process 1200, or portions thereof, may be implemented by an iterative gVCF genotyper and executed on one or more processors of one or more computing devices. Process 1200, or portions thereof, may be performed on a single computing device or distributed across multiple computing devices (e.g., multiple servers, sequencing devices, and / or client computing devices). Process 1200, or portions thereof, may reduce the amount of processing resources and / or storage space (e.g., memory) used by a computing device during variant analysis.
[0118] The exemplary process 1200 may begin after an iterative gVCF genotyper generates a cohort file and / or census file to be compressed. As described further herein, the fields in the cohort file and the census file may each include a subset of the fields contained in the gVCF file. Each of the fields may include genomic variant data and / or other genotype data associated with the batch of gVCF files. In the example provided in FIG. 12A, the compression process 1200 is provided to compress allele data for the RAGT field in the cohort file 1202. However, similar types of compression may be performed on the cohort data and / or census data for each field of the corresponding cohort file and / or census file, and will be further described elsewhere herein.
[0119] As described in the exemplary process 1200, the iterative gVCF genotyper can encode data in one or more fields of the cohort file such that the data can be compressed into a bit array 1218. As described herein, each field in the cohort file can include a set of unique values or metrics for the field, each followed by an integer array of gVCF files or samples that include the unique values or metrics. As shown in FIG. 12A, the RAGT field in the cohort file 1202 can include a set of unique values or metrics 1204 previously identified in the batch of samples from which the cohort file was generated. The unique values or metrics 1204 can relate to one or more fields in the RAGT statistics, such as, for example, reference alleles, alternative alleles, and / or genotypes for each of the samples associated with the variant data. Each unique value or metric 1204 can be followed by an integer array that identifies one or more samples or gVCF files 1206 that have the preceding unique value or metric 1204.
[0120] The iterative gVCF genotyper may encode the data in the RAGT field using a bitmap 1208 such that the unique values 1204 and corresponding samples or gVCF files may be represented as a series of bits. To generate the bitmap 1208, the iterative gVCF genotyper may identify a total number of unique values or metrics 1204 in the RAGT field. The iterative gVCF genotyper may identify a binary value length of the number of bits 1216 that may be implemented to represent the total number of unique values or metrics 1204. For example, the total number of unique values or metrics 1204 in the RAGT field of the cohort file 1202 is seven. The iterative gVCF genotyper may identify the lowest number of bits that may be implemented to represent seven unique values in the RAGT field of the cohort file 1202. The bits 1216 may be set to a three bit sequence because a three bit sequence is the lowest number of bits that may be implemented to represent a unique value in the RAGT field. It will be appreciated that other binary value lengths for the number of bits 1216 may be implemented to represent other numbers of unique values or metrics.
[0121] For RAGT fields, the iterative gVCF genotyper can identify a plurality of reference alleles and / or a plurality of alternative alleles, for example, using RAGT statistics in a plurality of gVCF files. The iterative gVCF genotyper can determine which of a plurality of samples have a common reference allele and an alternative allele. These samples can be distributed into allele groups. Each allele group can include one or more samples that have a common reference allele and an alternative allele. For example, when two samples are found to have a common (e.g., the same) reference allele and an alternative allele, the two samples can be grouped together into an allele group. Each allele group can be included in a unique value or metric 1204. The allele groups can be aggregated into a list of allele groups represented in the unique value or metric 1204. The iterative gVCF genotyper can determine the number of bits 1216 based on the number of allele groups (e.g., values). For example, the iterative gVCF genotyper may select a binary value length based on the number of allele groups.
[0122] The bitmap 1208 may include unique values or metrics 1204 and corresponding identifiers of samples or gVCF files 1206 having the unique values or metrics 1204. Each of the values or metrics 1204 and corresponding identifiers of the samples or gVCF files 1206 may be listed in the same row in the bitmap 1208 as its corresponding representative bit 1216. In the bitmap 1208, the bits 1216 are incremented with each row in the bitmap 1208, although other implementations will be appreciated. The integer values 1214 identify the integer values of the bits 1216 in each row.
[0123] The iterative gVCF genotyper may generate a bit array 1218 for the multiple gVCF files by aggregating the unique assigned binary values (e.g., unique values or metrics 1204). The bitmap 1208 may be used by the iterative gVCF genotyper to encode data in the RAGT field of the cohort file 120 into the bit array 1218 and to decode the bit array 1218 to generate data in the RAGT field of the cohort file 1202. The iterative gVCF genotyper may encode the bit array 1218 based on the order of associated samples in the multiple gVCF files. For example, the iterative gVCF genotyper may generate a bit array that includes a series of 3-bit sequences that represent the values 1204 of the samples 1206 in numerical order. As shown in FIG. 12A, the bit array 1218 may start with bit 1216 that represents the RAGT value 1204 for sample zero (e.g., "000"). The iterative gVCF genotyper may add a bit 1216 representing the RAGT value 1204 for sample 1 (e.g., "0001"), etc. The iterative gVCF genotyper may store a bit array 1218.
[0124] Bitmaps may differ between cohort and census files. Bitmaps may be generated per record and per metric field, as the bitmap may depend on the number of unique metric values and sample identifiers. In cohort files, the bitmap may encode hashed values per sample. In an example using the GQ field, a value of 10 in sample 1 may be hashed to 0, a value of 20 in sample 2 may be hashed to 1, and a value of 10 in sample 3 may be hashed to 0. The bitmap may encode an array [0,1,0] with the order of sample identifiers having value strings "10,20" (two unique values). In census files, the bitmap may encode sample counts per metric value. In another example using the GQ field, 25 samples with value 10, and 36 samples with value 20, the bitmap may encode an array [25,36] with the order of value strings "10,20".
[0125] In the example shown in FIG. 12A, the process 1200 can compress 211 bytes of input into 9 bytes of output. The exemplary process 1200 illustrates an example of compressing a metric for a field (e.g., RAGT) that has seven unique values that can be compressed using three bits (e.g., 000 to 110) to represent these seven values. The number of bits varies depending on the number of unique values to compress, e.g., 2 unique values=1 bit, 256 unique values=8 bits. Thus, compression can vary from record to record and metric to metric. The data being compressed can differ between cohorts and censuses. The exemplary process 1200 can be used by an iterative gVCF genotyper to compress one or more predefined fields in a gVCF file described herein. For example, the predefined fields can be compressed when storing data in cohort files and census files.
[0126] Although a process 1200 for performing compression on RAGT statistics is shown in FIG. 12A, it will be understood that the process 1200 is not limited to compressing RAGT statistics. For example, the process 1200 may be used to compress data in other fields of a gVCF file. Although the process 1200 is shown and described with reference to FIG. 12A as compressing data in a cohort file, it will be understood that a similar process may be implemented to compress data in a census file. Although the process 1200 is shown and described as being implemented to compress sequencing data, it will be further understood that the process 1200 is not limited to compressing sequencing data. Instead, the process 1200 may be used to compress other types of data.
[0127] FIG. 12B illustrates an exemplary process 1250 that may be implemented by one or more computing devices (such as, for example, the server device 102, the client device 108, and / or the sequencing device 114 shown in FIG. 1) to encode the cohort and / or census data in a compressed format. The iterative process 1250, or portions thereof, may be implemented to efficiently store the variant data in a census file, a cohort file, and / or an msVCF file. The process 1250, or portions thereof, may be implemented by an iterative gVCF genotyper and executed on one or more processors of one or more computing devices. The process 1250, or portions thereof, may be implemented on a single computing device or distributed across multiple computing devices (e.g., multiple servers, sequencing devices, and / or client computing devices). The process 1250, or portions thereof, may reduce the amount of processing resources and / or storage space (e.g., memory) used by a computing device during variant analysis.
[0128] The exemplary process 1250 may begin after an iterative gVCF genotyper generates a cohort file and / or census file to be compressed. As described further herein, the fields in the cohort file and the census file may each include a subset of the fields contained in the gVCF file. Each of the fields may include genomic variant data and / or other genotype data associated with the batch of gVCF files. In the example provided in FIG. 12B, a compression process 1250 is provided to compress allele data for the RAGT field in the census file 1252. However, similar types of compression may be performed on the cohort data and / or census data for each field of the corresponding cohort file and / or census file, and will be further described elsewhere herein.
[0129] As described in the exemplary process 1250, the iterative gVCF genotyper can encode data in one or more fields of the census file such that the data can be compressed into a bit array 1268. As described herein, each field in the census file can include a set of unique values or metrics for the field, each followed by a sample count of the gVCF files or samples that include the unique value or metric. As shown in FIG. 12B, the RAGT field in the census file 1252 can include a set of unique values or metrics 1254 (such as, for example, alleles shown in FIG. 12B) previously identified in the batch of samples from which the census file was generated. The unique values or metrics 1254 can be associated with one or more fields in the RAGT statistics, such as, for example, reference alleles, alternative alleles, and / or genotypes for each of the samples associated with the variant data. Each unique value or metric 1254 can be followed by a sample count 1256 that indicates how many of the one or more samples or gVCF files 1252 have the preceding unique value or metric 1254.
[0130] The iterative gVCF genotyper may encode the data in the RAGT field using a bitmap 1258 such that the unique values 1254 and sample counts 1256 may be represented as a series of bits. To generate the bitmap 1258, the iterative gVCF genotyper may identify a total number of unique values or metrics 1254 in the RAGT field of the census file. The iterative gVCF genotyper may identify a binary value length of the number of bits 1266 that may be implemented to represent the total number of unique values or metrics 1254. The iterative gVCF genotyper may identify a minimum number of bits that may be implemented to represent unique values (e.g., unique alleles) in the RAGT field of the census file 1252. The bits 1266 may be set to a three bit sequence because a three bit sequence is the minimum number of bits that may be implemented to represent unique values in the RAGT field. It should be understood that other binary value lengths for the number of bits 1266 may be implemented to represent other numbers of unique values or metrics.
[0131] For RAGT fields, the iterative gVCF genotyper can use RAGT statistics in multiple gVCF files to identify sample counts for alleles. For example, when two samples are found to have a common (e.g., the same) reference allele and an alternative allele, the sample count for that allele (e.g., allele group) can be listed in sample count 1256. The allele groups can be aggregated into a list of allele groups represented in sample count 1256. The iterative gVCF genotyper can determine the number of bits 1266 based on the number of allele groups (e.g., values). For example, the iterative gVCF genotyper may select the binary value length based on the number of alleles / allele groups.
[0132] The bitmap 1258 may include unique values or metrics 1254 and sample counts of samples or gVCF files 1256 having unique values or metrics 1254. Each of the values or metrics 1254 and sample counts of samples or gVCF files 1256 may be listed in the same row in the bitmap 1258 as its corresponding representative bit 1266. In the bitmap 1258, the bit 1266 is incremented with each row in the bitmap 1258, although other implementations should be understood. The integer value 1264 identifies the integer value of the bit 1266 in each row.
[0133] The iterative gVCF genotyper may generate a bit array 1268 for the multiple gVCF files by aggregating the unique assigned binary values (e.g., unique values or metrics 1254). The bitmap 1258 may be used by the iterative gVCF genotyper to encode data in the RAGT field of the census file 1252 into the bit array 1268 and to decode the bit array 1268 to generate data in the RAGT field of the census file 1252. The iterative gVCF genotyper may encode the bit array 1268 based on the order of associated samples in the multiple gVCF files. For example, the iterative gVCF genotyper may generate a bit array that includes a series of 3-bit sequences that represent sample counts 1256 of alleles / allele groups 1254 in numerical order. 12B, the bit array 1268 can start with bit 1266 representing the sample count 1256 for allele 1 (e.g., TACAC:TAC,T,.:1 / 2). The iterative gVCF genotyper can append bit 1266 representing the sample count 1256 for allele 2 (e.g., TACACAC:T,TACAC,.:12), etc. The iterative gVCF genotyper can store the bit array 1268.
[0134] A bitmap 1258 may be generated per record and per metric field, as the bitmap may depend on the number of unique metric values and sample identifiers. In a census file, the bitmap may encode the sample count. In an example using the GQ field, 25 samples with value 10, and 36 samples with value 20, the bitmap may encode the array [25,36] in the order of the value string "10,20".
[0135] In the example shown in FIG. 12B, the process 1250 can compress 211 bytes of input into 6 bytes of output. The exemplary process 1250 illustrates an example of compressing a metric for a field (e.g., RAGT) that has seven unique values that can be compressed using three bits (e.g., 000 to 110) to represent these seven values. The number of bits varies depending on the number of unique values to compress, e.g., 2 unique values=1 bit, 256 unique values=8 bits. Thus, compression can vary from record to record and metric to metric. The data being compressed can differ between cohorts and censuses. The exemplary process 1250 can be used by an iterative gVCF genotyper to compress one or more predefined fields in a gVCF file described herein. For example, the predefined fields can be compressed when storing data in cohort files and census files.
[0136] FIG. 13 is a flow diagram illustrating an example method 1300 that may be implemented by one or more computing devices (such as, for example, the server device 102, the client device 108, and / or the sequencing device 114 shown in FIG. 1) to perform compression of cohort and / or census data. The method 1300, or portions thereof, may be performed by an iterative gVCF genotyper to efficiently store genomic variant data and / or other genomic data in census files, cohort files, and / or msVCF files. The method 1300, or portions thereof, may be performed on a single computing device or distributed across multiple computing devices (e.g., multiple servers, sequencing devices, and / or client computing devices). The method 1300, or portions thereof, may reduce the amount of processing resources and / or storage space (e.g., memory) used by a computing device during variant analysis.
[0137] As shown in FIG. 13, the iterative gVCF genotyper can identify fields in a cohort file and / or a census file for compression at 1302. The fields in the cohort file and / or the census file can each include a subset of the fields included in the gVCF file. Each of the fields can include genomic variant data and / or other genotype data associated with a batch of gVCF files. The fields can be formatted as text fields. As described herein, each field in the cohort file can include a set of unique values or metrics for the field, each followed by an integer array of gVCF files or samples that include the unique values or metrics. Also, as described herein, each field in the census file can include a set of unique values or metrics for the field, each followed by a sample count for the number of samples with unique values or metrics in the batch of samples represented by the census file.
[0138] At 1304, the iterative gVCF genotyper can generate a bitmap for encoding the cohort or census data for the identified field. To generate the bitmap, the iterative gVCF genotyper can identify a total number of unique values or metrics in the field. The iterative gVCF genotyper can identify a binary value length for the number of bits that can be implemented to represent the total number of unique values or metrics in the field of the cohort or census file. For example, the iterative gVCF genotyper can identify and select the minimum number of bits that can be implemented to represent the total number of unique values or metrics in the field of the cohort or census file in the bitmap.
[0139] For cohort files, the bitmap may include unique values or metrics for a field and corresponding identifiers of samples or gVCF files that have the unique values or metrics. For census files, the bitmap may include unique values or metrics for a field and corresponding sample counts for the number of samples that have the unique values or metrics in the batch of samples represented by the census file. Each value or metric may be listed in the same row in the bitmap as its corresponding representative bit.
[0140] At 1306, the iterative gVCF genotyper may use a bitmap to encode the cohort or census data for the identified fields in a bit array. The bitmap may be used by the iterative gVCF genotyper to encode the data in the fields of the cohort or census file into a bit array. The iterative gVCF genotyper may encode the bit array for the cohort data based on the order of the associated samples in the multiple gVCF files such that the bits representing unique values or metrics are encoded in ascending order by sample or gVCF identifier. The census data in the census file may be ordered based on the metric values for each field. For metric A, if values are v1, v2, v3 and counts c1, c2, c3, then count [c1, c2, c3] is bit compressed and the number of bits depends on the maximum of c1, c2, c3.
[0141] At 1308, the iterative gVCF genotyper may determine whether additional fields in the cohort or census files should be compressed, or whether additional cohort or census files should be compressed. If there are additional fields or files to compress, method 1300 may return to 1304. If there are no additional fields or files to compress, method 1300 may end.
[0142] The compressed cohort data and the compressed census data may be further compressed and serialized as shown in FIG. 14. FIG. 14 shows an example process 1400 that may be implemented by one or more computing devices (such as, for example, the server device 102, the client device 108, and / or the sequencing device 114 shown in FIG. 1) for data serialization and compression. The process 1400, or portions thereof, may be performed by an iterative gVCF genotyper. In the example process 1400, the cohort data and the census data may be serialized into binary data using bit compression, and may be further encoded into ascii data so that the data (along with the region information) may be compressed using the bgzip compression algorithm, enabling genomic region querying, and indexed by the tabix algorithm. bgzip and tabix are public algorithms implemented in HTSlib. The process 1400, or portions thereof, may be performed to efficiently store genomic variant data and / or other genotype data in census files, cohort files, and / or msVCF files. Process 1400, or portions thereof, may be performed on a single computing device or distributed across multiple computing devices (e.g., multiple servers, sequencing devices, and / or client computing devices). Process 1400, or portions thereof, may reduce the amount of processing resources and / or storage space (e.g., memory) used by a computing device during variant analysis.
[0143] The process 1400 can directly encode the gVCF data into bitwise binary data (instead of ascii serialized data) using a binary data compression scheme (such as, for example, the process 1200 shown in FIG. 12). The process 1400 can implement alternative compression algorithms (e.g., Iz4 instead of bgzip) and customized indexing algorithms (e.g., instead of tabix) to further improve compression and / or query efficiency, for example. Additionally or alternatively, the example process 1400 can incorporate one or more other compression methods.
[0144] Cohort data 1410 (e.g., in one or more cohort files) may be received at an iterative gVCF genotyper. Census data 1412 (e.g., in one or more census files and / or a global census file) may be received at an iterative gVCF genotyper. The cohort data and / or census data may be separated by region. Bit compression may be performed on the cohort and census data to generate a compressed file 1414. For example, the cohort and census data may be represented using a bit array (e.g., a binary bit array) (e.g., using a process similar to process 1200 shown in FIG. 12). The cohort or census file headers including the field values may be serialized in the serialized data file 1416. The cohort or census file headers including the field values may be followed by serialized data for each field in the cohort and census files. The headers and / or data values may be serialized using one or more field separators (e.g., key?val&, key?vall, val2,...&..., and / or key?skeyl:svall, skey2:sval2, ...&..., etc.).
[0145] Cohort and census data may be encoded differently (e.g., with sample identifier lists or sample counts) as described herein. Both encoded cohort and census files store separate copies of value lists (e.g., vail l, vall2...). To serialize the metrics stored in each file, we take the order of the file names defined in the cohort and census headers, respectively, to avoid repeating the file names in each record. For example, the data for each of the metrics, such as field 1 (e.g., having metric values v11, v12, v13, and bit array b1), field 2 (e.g., having metric value v21, bit array b2), and field 3 (e.g., having metric values v31, v32, and bit array b3), are serialized together and the serialized data is v11&v12&v13?b1;v21&b2;v31&v32?b3. The characters "&", "?", may be used to separate the fields as other characters may be used for metric values in a VCF file.
[0146] This serialized data for each record may be stored in columns for each row in the cohort and census files. Each file may also include separate columns containing the chromosome name, start position, and end position of the record. These columns may be implemented by the browser extensible data (BED) format. Additional columns for customized data may be included that may be used to store the cohort and census record data. Each row of the cohort and census files may be further compressed. For example, each row of the cohort and census files may be further compressed in blocks (e.g., 64 kilobyte blocks) using the public bgzip algorithm. The compressed blocks may be indexed by genomic region using an indexing algorithm. An exemplary indexing algorithm may be the public tabix algorithm.
[0147] The iterative gVCF genotyper can generate a serialized data file 1418 in BED format with bytes of serialized data. The BED format can be a text file format used to store genomic regions as coordinates and associated annotations. At 1420, the bit compressed and serialized cohort and census hash map data can be further encoded into ascii data so that the data (along with the region information) can be compressed using the bgzip compression algorithm and indexed by the tabix algorithm. The process 1400 (e.g., the bgzip compression algorithm and the tabix algorithm) can enable genomic region querying.
[0148] FIG. 15 illustrates another exemplary process 1500 that may be implemented by one or more computing devices (such as, for example, the server device 102, the client device 108, and / or the sequencing device 114 shown in FIG. 1) for data serialization and compression. The process 1500, or portions thereof, may be performed by an iterative gVCF genotyper and executed by a processor in one or more computing devices. In the exemplary process 1500, the cohort data and census data are serialized into binary data using bit compression, and may be further encoded into ascii data so that the data (along with the region information) may be compressed using the bgzip compression algorithm, which may enable genomic region querying, and indexed by the tabix algorithm. bgzip and tabix are public algorithms implemented in HTSlib. The process 1500, or portions thereof, may be implemented to more efficiently store variant data in census files, cohort files, and / or msVCF files. Process 1500, or portions thereof, may be performed on a single computing device or distributed across multiple computing devices (e.g., multiple servers, sequencing devices, and / or client computing devices). Process 1500, or portions thereof, may reduce the amount of processing resources and / or storage space (e.g., memory) used by a computing device during variant analysis.
[0149] In the exemplary process 1500, cohorts and censuses may be combined into bit arrays (e.g., hex representations). The exemplary process 1500 may include JSON lists and bit / hex representations of numbers. The exemplary process 1500 may support adding and removing metrics. The exemplary process 1500 may support generic metrics (strings, numbers). The exemplary process 1500 may support values including periods, commas, colons, slashes, etc., and may ensure that delimiters are appropriately selected. The exemplary process 1500 may be compatible with BGZF (e.g., BGZF blocks) and tabix (e.g., tabix indexes).
[0150] The iterative gVCF genotyper may receive cohort data 1510 in one or more cohort files. The iterative gVCF genotyper may receive census data 1512 in one or more census files and / or global census files. The cohort data and / or census data may be separated by location or block within the cohort and / or census files. Bit compression may be performed on the cohort and census data at 1514. For example, the cohort and census data may be represented using bit arrays (e.g., binary bit arrays) (e.g., using process 1200 shown in FIG. 12). One or more headers and / or alphanumeric values in the cohort or census file headers may be serialized at 1516. The headers and / or alphanumeric values may be serialized using one or more field separators (e.g., key?val&, key?va11,va12,....&..., and / or key?skey1:sval1,skey2:sval2,...&..., etc.). At 1518, a byte buffer may be generated. The byte buffer may contain the cohort and / or census data. At 1520, the bit-compressed and serialized cohort and census data may be further encoded into ascii data such that the data (along with the region information) may be compressed using the bgzip compression algorithm and indexed by the tabix algorithm. The process 1500 (e.g., the bgzip compression algorithm and the tabix algorithm) may enable genomic region querying.
[0151] In addition to or instead of implementing data compression and other memory storage techniques described herein, data in files described herein can be stored in memory in a manner that preserves memory locations for overlap data. For example, variant data in cohort files and census files can include overlap data. The overlap data can refer to overlaps of genomic regions (e.g., overlapping variants and / or hom-ref blocks). In one example, a first record for chr20:1000-2000 and a second record for chr20:1500-2500 can have an overlap at chr20:1500-2000. When these records are aggregated, a third record can be generated for three regions: chr20:1000-1499, chr20:1500-2000, chr:2001-2500. Cohort or census data in the first and third regions may be copied from the original first record and the original second record, while data from the second region may be updated from both the first record and the second record. The variant data may be received in a genome variant call file that includes records representing genomic regions. The variant data may be stored in a memory location that improves storage of variant data by leveraging overlapping regions of different records.
[0152] FIG. 16 illustrates an exemplary process 1600 that may be implemented by one or more computing devices (such as, for example, the server device 102, the client device 108, and / or the sequencing device 114 shown in FIG. 1) to add, copy, and / or update records utilizing buffer aggregation with overlapping genomic regions. For example, the exemplary process 1600, or portions thereof, may be performed by an iterative gVCF genotyper and executed by a processor in one or more computing devices. The method 1600 may be implemented to aggregate gVCF files into cohort files and / or aggregate census files into a global census file. The process 1600, or portions thereof, may be implemented to efficiently store genomic variant data in census files, cohort files, and / or msVCF files.
[0153] The exemplary process 1600 can begin when an iterative gVCF genotyper receives one or more records. The records can be received from one or more input gVCF files or one or more input census files. Each of the records can include genomic variant data representing a different genomic region.
[0154] Buffer 1602 may be an output buffer that includes a certain number of buffer positions (bps) for storing data in a record in each of the available buffer positions. The iterative gVCF genotyper may receive a first set of records including genomic variant data for storage in buffer 1602. The first set of records may include record 1, record 2, record 3, and record 4. Each record may indicate a number of buffer positions available for storing genomic regions in the genomic variant data. In the example shown in FIG. 16, the genomic variant data of record 1 may occupy three buffer positions, the genomic variant data of record 2 may occupy one buffer position, the genomic variant data of record 3 may occupy two buffer positions, and the genomic variant data of record 4 may occupy five buffer positions.
[0155] The iterative gVCF genotyper can analyze each of the received records (e.g., records 1-4) and the buffer 1602 to determine an operation to be performed on each record to enable storage in the buffer 1602. For example, the iterative gVCF genotyper can analyze the received records and the buffer 1602 to determine whether there are any genomic regions within the received record that overlap with previously received genomic variant data stored in the buffer 1602. If the iterative gVCF genotyper determines that a genomic region within the record or within one or more portions of the record does not overlap with data currently stored in the buffer 1602, the iterative gVCF genotyper can perform an operation to add the record, or those portions that do not overlap with previously stored data, to the buffer (e.g., as indicated by operation "A" in each record). Now, because each of the buffer locations in buffer 1602 is open and does not have any data stored therein, there are no previously stored records that may overlap any of records 1-4 or portions thereof, and the iterative gVCF genotyper can add each of records 1-4 to buffer 1602. The iterative gVCF genotyper can store genomic variant data for each of records 1-4 in buffer 1602 with at least one empty buffer location between non-overlapping adjacent records. The records may be separated by empty base pairs (e.g., genomic locations) when there is no overlapping data. After storing records 1-4 in buffer 1602, the iterative gVCF genotyper can receive a second set of records to be stored in buffer 1602. The second set of records may include record 5, record 6, record 7, and record 8. Again, each record may indicate the number of buffer locations that may be utilized to store genomic regions in the genomic variant data. Because the first batch of records does not have overlapping data, as shown in the range of records 1, 2, 3, and 4 below buffer 1602, the first batch of records may be separated by empty base pairs (e.g., genomic locations), as shown in FIG. 16.If there were duplicates in some positions in the first batch of records, the records may be stored similarly to the duplicates shown in records 5-8.
[0156] The iterative gVCF genotyper may parse each of the received records (e.g., records 5-8) and parse buffer 1602 to determine an action to be performed on each record to enable storage in buffer 1602. If a genomic region in a portion of the record does not overlap with a record previously stored in buffer 1602, the iterative gVCF genotyper may add genomic variant data including the genomic region to buffer 1602 in the respective buffer location. If a genomic region in at least a portion of the genomic variant data in the record completely overlaps with a genomic region of genomic data in a record previously stored in buffer 1602, the overlap buffer location in buffer 1602 may be updated to include the overlap portion of the updated record. If a genomic region in at least a portion of the genomic variant data in the record partially overlaps with a genomic region of genomic data in a record previously stored in buffer 1602, the overlap buffer location in the buffer may be updated and the non-overlapping portion of the record may be copied and stored with the overlap portion to maintain the genomic variant data stored in the previous record.
[0157] For example, as shown in FIG. 6, record 5 may include a portion of genomic variant data that may occupy three buffer locations that overlap with the entirety of previously stored record 1. The portion of genomic variant data in record 5 may include a genomic region that overlaps with a genomic region in the genomic variant data of record 1. As such, the iterative gVCF genotyper may identify that a genomic region in the genomic variant data of record 5 completely overlaps with a genomic region previously stored for record 1 at buffer location 1602b, and buffer location 1602b may be updated with the data in record 5 (e.g., as indicated by action "U" in record). Because record 5 includes an additional portion of genomic variant data that includes a genomic region that does not overlap with other genomic variant data in buffer 1602, the genomic variant data may be added to buffer location 1602a, which may be a contiguous buffer location with updated buffer location 1602b. For example, one or more of the updated buffer locations 1602b may be adjacent to buffer location 1602a. Upon processing record 6, the iterative gVCF genotyper can also identify that the genomic variant data in record 6 completely overlaps with the genomic variant data previously stored for record 2 in buffer location 1602c, and buffer location 1602c can be updated with the data in record 6.
[0158] The iterative gVCF genotyper may process record 7 and identify that a portion of the genomic variant data in record 7 overlaps with a portion of the genomic variant data in previously stored record 3. Because the portion of the previously stored record 3 is maintained, the iterative gVCF genotyper may copy the portion of the genomic variant data for record 3 and update the portion of the genomic variant data for record 3 (e.g., indicated by the action "C,U" in the record). The genomic variant data previously stored in buffer location 1602e for record 3 may be updated with the data in record 7. The non-overlapping portion of that of the previously stored record 3 may be copied and maintained in buffer location 1602d, which is a contiguous buffer location with buffer location 1602e, which contains overlapping genomic variant data for the two records. Because record 7 also contains a non-overlapping portion of a record, the non-overlapping portion of record 7 may be added to buffer location 1602f.
[0159] The iterative gVCF genotyper may process record 8 and identify that the genomic variant data in record 8 overlaps with a portion of the genomic variant data in record 4 that was previously stored. Because the portion of the previously stored record 4 is maintained, the iterative gVCF genotyper may copy the portion of the genomic variant data for record 4 and update the portion of the genomic variant data for record 4 (e.g., indicated by the action "C,U" in the record). The genomic variant data previously stored in buffer location 1602g for record 4 may be updated with the genomic variant data in record 7. The non-overlapping portion of the previously stored record 4 may be copied and maintained in buffer location 1602h, which may be a contiguous buffer location with buffer location 1602g that contains the overlapping genomic variant data for the two records. Although the exemplary process 1600 shown in FIG. 16 uses an output buffer having 20 buffer locations, it will be understood that the output buffer is not limited to having 20 buffer locations. Instead, the output buffer may have any number of buffer locations.
[0160] Although features, elements, and functions are described above in particular combinations, the features, elements, or functions may be used alone or in any combination with other features, elements, or functions. Various presently unforeseen or unanticipated alternatives, modifications, variations, or improvements may be made thereafter, which are also intended to be encompassed by the following claims.
[0161] The methods described herein are implemented in a computer program, software, or firmware embodied in a computer-readable medium for execution by a computer or processor. Examples of computer-readable media include electronic signals (transmitted by wired or wireless connections) and computer-readable storage media. Examples of computer-readable storage media include, but are not limited to, read only memory (ROM), random-access memory (RAM), removable disks, and optical media such as CD-ROM disks and digital versatile disks (DVDs). [Explanation of symbols]
[0162] 100 System Environment 102 Server Device 104 Sequencing System 108 client devices 110 Sequencing Applications 112 Network 114 Sequencing Device 116 Databases 200 computing devices 202 Processor 204 Memory 206 Storage 208 I / O Interface 210 Communication Interface 212 Communications Infrastructure
Claims
1. 1. A computer-implemented method for iterative gVCF genotyping, comprising: receiving a first plurality of genome variant call files associated with a first batch of sequencing data; generating a first cohort file for the first batch by aggregating data from a subset of fields in each of the first plurality of genomic variant call files; generating a first census file containing variant summary statistics and hom-ref blocks for the first batch; receiving a second plurality of genome variant call files associated with a second batch of sequencing data; generating a second cohort file for the second batch by aggregating data from a subset of fields in each of the second plurality of genomic variant call files; generating a second census file containing variant summary statistics and hom-ref blocks for the second batch; generating a global census file by aggregating the first census file and the second census file, the global census file including census data from batches of samples received from sequencing devices at different sites; generating a first multi-sample variant call file for the first batch using the first cohort file, the first census file, and the global census file; generating a second multi-sample variant call file for the second batch using the second cohort file, the second census file, and the global census file.
2. performing a genome-wide sequencing analysis using one or more of the first multi-sample variant call file or the second multi-sample variant call file; The computer-implemented method of claim 1 further comprising:
3. 2. The computer-implemented method of claim 1, wherein the first plurality of genome variant call files associated with the first batch are divided into equal-sized shards, and each shard is processed using one of a plurality of computation nodes.
4. receiving a third plurality of genome variant call files associated with a third batch of sequencing data; generating a third cohort file for the third batch by aggregating data from a subset of fields in each of the third plurality of genomic variant call files; generating a third census file containing variant summary statistics and hom-ref blocks for the third batch; updating the global census file by aggregating the third census file with the global census file; generating a third multi-sample variant call file for the third batch using the third cohort file, the third census file, and the updated global census file; The computer-implemented method of claim 1 further comprising:
5. The computer-implemented method of claim 1 , wherein the computer-implemented method is performed on a local computing system or distributed across a cloud computing system.
6. 1. A system comprising: at least one processor; a computer-readable medium containing instructions; The instructions, when executed by the at least one processor, cause the at least one processor to: receiving one or more genomic variant call files associated with one or more samples; generating one or more cohort files and one or more census files from the one or more genome variant call files; aggregating the one or more census files into a global census file, the global census file including census data from batches of samples received from sequencing devices at different sites; generating at least one multi-sample variant call file based on the global census file, the one or more cohort files, and the one or more census files; storing the multi-sample variant call file in a memory.
7. The system of claim 6 , wherein the samples comprise samples from a sequencing run, a sequencing cycle, or multiple sequencing runs.
8. The system of claim 6 , wherein the instructions cause the at least one processor to perform parallel processing using multiple compute nodes.
9. 9. The system of claim 8, wherein the instructions further cause the at least one processor to perform parallelization and multithreading by regions of array data.
10. 10. The system of claim 9, wherein at least two computing nodes implement at least two levels of parallelism to process, aggregate, or generate data for corresponding regions of the sequence data, each computing node processing a particular region.
11. The system of claim 6 , wherein the cohort files and the census files within a region are bit-compressed and serialized.
12. 1. A system comprising: at least one processor; a computer-readable medium containing instructions; The instructions, when executed by the at least one processor, cause the at least one processor to: receiving a plurality of genomic variant call files, each of the genomic variant call files associated with a respective sample of a plurality of samples; identifying a plurality of reference alleles and a plurality of alternative alleles associated with the plurality of samples using reference alternative genotype (RAGT) statistics in the plurality of genomic variant call files; counting instances of each of the plurality of reference alleles and each of the plurality of alternative alleles; selecting a normalized reference allele from the plurality of reference alleles, wherein the longest reference allele is selected as the normalized reference allele; normalizing the other reference alleles of the plurality of reference alleles by extending the other reference alleles to correspond to the normalized reference allele; normalizing the plurality of alternative alleles by extending each alternative allele by the same amount as its corresponding reference allele was extended; generating a multi-sample variant call file using the normalized reference allele and the normalized alternative allele.
13. 13. The system of claim 12, wherein the instructions further cause the at least one processor to generate a normalized representation of each sample using the normalized reference allele, such that each of the plurality of alternative alleles is indexed using the normalized reference allele.
14. 13. The system of claim 12, wherein the other reference alleles are extended by adding a respective number of bases to correspond to the normalized reference allele.
15. The instructions cause the at least one processor to: receiving additional genomic variant call files associated with the additional samples; and identifying a reference allele and one or more alternative alleles associated with the additional sample; and and updating the normalized representation to include the reference allele and the one or more alternative alleles associated with the additional sample.
16. The instructions configured to cause the at least one processor to update the normalized representation may include instructions to the at least one processor to: determining that the length of the reference allele is shorter than the normalized reference allele; extending the reference allele and the one or more alternative alleles to correspond to the normalized reference allele; and reordering the plurality of reference alleles and the plurality of alternative alleles to include the extended reference allele and the one or more extended alternative alleles.
17. The instructions configured to cause the at least one processor to update the normalized representation may include instructions to the at least one processor to: determining that the length of the reference allele is longer than the normalized reference allele; selecting the reference allele as an updated normalized reference allele; normalizing the plurality of reference alleles and the plurality of alternative alleles by extending them to correspond to the length of the updated normalized reference alleles; and reordering the plurality of reference alleles and the plurality of alternative alleles to include the updated normalized reference allele and the one or more extended alternative alleles.
18. 13. The system of claim 12, wherein the instructions cause the at least one processor to further reorder the genotypes of each sample based on the normalized reference alleles and the normalized representation.
19. 20. The system of claim 18, wherein the instructions further cause the at least one processor to generate a mapping for each of the plurality of alternative alleles based on the normalized reference allele.
20. 20. The system of claim 19, wherein the mapping for each of the plurality of alternative alleles is stored in site information within a census file.