Systems and methods for generating graph references - Patents.com
Patent Information
- Application Number
- JP2023557038
- Authority / Receiving Office
- JP · JP
- Patent Type
- Applications
- Current Assignee / Owner
- Priority Date
- 2021-03-17
- Filing Date
- 2022-03-17
- Publication Date
- 2025-03-26
AI Technical Summary
Traditional techniques for aligning sequence reads against graph reference constructs yield inaccurate and computationally expensive results due to the inclusion of structural and multiple-alignable variants, which introduce ambiguity and increase computational complexity.
A method for generating graph reference constructs involves filtering out structural and multiple-alignable variants in multiple stages, using a first filtering stage to exclude structural variants and a second stage to remove multiple-alignable variants, thereby improving alignment accuracy and reducing computational complexity.
The proposed technique results in more accurate and computationally efficient alignment of sequence reads, enhancing the accuracy of variant calling and reducing computational costs.
Smart Images

Figure 00000000_0000_ABST
Abstract
Description
[Technical field]
[0001] CROSS-REFERENCE TO RELATED APPLICATIONS This application claims the benefit of priority under 35 U.S.C. §119(e) to U.S. Provisional Patent Application No. 63 / 162,400, entitled “SYSTEMS AND METHODS FOR GENERATING GRAPH SEQUENCES,” filed March 17, 2021, the entire contents of which are incorporated herein by reference.
[0002] Reference to sequence listings submitted as text files via EFS-WEB This application contains a Sequence Listing that was submitted via EFS-Web in ASCII format and is incorporated herein by reference in its entirety. The above ASCII copy, created on March 15, 2022, is named S196170030WO00-SEQ-DGR and is 5,033 bytes in size. [Background technology]
[0003] background Advances in sequencing technology, including the development of next-generation sequencing techniques, have made sequencing an important tool used in both research and medicine. Some applications of sequencing technology include aligning sequence reads obtained by sequencing techniques to a reference sequence construct and identifying differences, sometimes referred to as "variants", between the sequence reads and the reference sequence construct. As a result, the identified differences can be used for diagnostic, prognostic, therapeutic, research, and / or other purposes.
[0004] There are different kinds of reference sequence constructs to which sequence reads can be aligned. For example, sequence reads can be aligned to linear reference sequence constructs, such as hg19 and hg38 human reference genomes. As another example, sequence reads can be aligned to a reference sequence construct that describes one or more known variants at one or more respective locations. One example of such a reference sequence construct is a graph-based reference sequence construct (sometimes referred to herein as a "graph reference construct"). A graph reference construct can include a graph (e.g., a directed acyclic graph) in which there can be multiple paths, each of which can represent one or more known variants. Summary of the Invention [Means for solving the problem]
[0005] overview Some embodiments are directed to a method for generating a graph reference construct, the method comprising: obtaining, with at least one computing device, a plurality of variants associated with a reference sequence construct for at least a portion of a genome; and generating a graph reference construct using the plurality of variants and the reference sequence construct, wherein the generating comprises filtering the plurality of variants to obtain a filtered set of variants, the filtered set of variants being a subset of the plurality of variants, and wherein the filtering comprises a plurality of filtering steps, the plurality of filtering steps including a first filtering step and a second filtering step different from the first filtering step and performed after the first filtering step. the first filtering step includes identifying a first subset of variants from the plurality of variants by, at least in part, excluding one or more structural variants from the plurality of variants, the one or more structural variants including a first structural variant, and the second filtering step includes identifying a filtered set of variants from the first subset of variants by, at least in part, excluding one or more multi-alignable variants from the first subset of variants, the filtering step includes generating a graph reference construct using the filtered set of variants and the reference sequence construct, and outputting the generated graph reference construct.
[0006] Some embodiments provide a system comprising at least one computer hardware processor and at least one non-transitory computer readable storage medium storing processor executable instructions, the processor executable instructions, when executed by the at least one computer hardware processor, for performing a process to cause the at least one computer hardware processor to: obtain a plurality of variants associated with a reference sequence construct for at least one portion of a genome; and generate a graph reference construct using the plurality of variants and the reference sequence construct, wherein the generating comprises filtering the plurality of variants to obtain a filtered set of variants, the filtered set of variants being a subset of the plurality of variants, and the filtering ... and generating a graph reference construct using the filtered set of variants and the reference sequence construct. The system further comprises: a plurality of filtering steps, including a second filtering step performed after the first filtering step, the first filtering step including identifying a first subset of variants from the plurality of variants by, at least in part, excluding one or more structural variants from the plurality of variants, the one or more structural variants including a first structural variant, and the second filtering step including identifying a filtered set of variants from the first subset of variants by, at least in part, excluding one or more multi-alignable variants from the first set of variants. The system further comprises:
[0007] Some embodiments include at least one non-transitory computer readable storage medium storing processor executable instructions, the processor executable instructions, when executed by the at least one computer hardware processor, for causing the at least one computer hardware processor to obtain a plurality of variants associated with a reference sequence construct for at least one portion of a genome; and generating a graph reference construct using the plurality of variants and the reference sequence construct, wherein the generating comprises filtering the plurality of variants to obtain a filtered set of variants, the filtered set of variants being a subset of the plurality of variants, and the filtering comprises a first filtering step and a second filtering step distinct from the first filtering step and performed after the first filtering step. and generating a graph reference construct using the filtered set of variants and the reference sequence construct. The present invention provides at least one non-transitory computer-readable storage medium for performing the generating and outputting of the generated graph reference construct, the filtering step including a first filtering step including a filtering step including a filtering step including a first filtering step including identifying a first subset of variants from the plurality of variants by, at least in part, excluding one or more structural variants from the plurality of variants, the one or more structural variants including the first structural variant, and a second filtering step including identifying a filtered set of variants from the first subset of variants by, at least in part, excluding one or more multi-alignable variants from the first set of variants.
[0008] In some embodiments, identifying the first subset of variants from the plurality of variants includes determining whether a first length of the first structural variant exceeds a first specified threshold, and removing the first structural variant from the plurality of variants upon determining that the first length exceeds the first specified threshold.
[0009] In some embodiments, the first structural variant is an insertion event, and determining whether a first length of the first structural variant exceeds a first specified threshold comprises determining whether the first length is at least 5,000 base pairs.
[0010] In some embodiments, the first structural variant is a deletion event, and determining whether a first length of the first structural variant exceeds a first specified threshold comprises determining whether the first length is at least 90,000 base pairs.
[0011] In some embodiments, identifying the first subset of variants from among the plurality of variants comprises aligning the first structural variant to a reference sequence construct.
[0012] In some embodiments, identifying the first subset of variants from the plurality of variants comprises determining whether the reference sequence construct comprises a subsequence, where the subsequence is identical to at least a portion of the first structural variant, and upon determining that the reference sequence construct comprises the subsequence, excluding the first structural variant from the plurality of variants.
[0013] In some embodiments, identifying the first subset of variants from among the plurality of variants comprises aligning the first structural variant to one or more variants of the plurality of variants, where the one or more variants are distinct from the first structural variant.
[0014] In some embodiments, identifying the first subset of variants from the plurality of variants includes determining whether the second structural variant comprises a subsequence, where the subsequence is identical to at least a portion of the first structural variant, and upon determining that the second structural variant comprises the subsequence, removing one of the first structural variant or the second structural variant from the plurality of variants.
[0015] In some embodiments, identifying the first subset of variants from among the plurality of variants comprises aligning the first structural variant to a decoy sequence associated with the reference sequence construct.
[0016] In some embodiments, identifying a first subset of variants from among the plurality of variants includes determining whether a decoy sequence associated with the reference sequence construct comprises a subsequence, which is identical to at least a portion of the first structural variant, and masking the decoy sequence upon determining that the decoy sequence comprises the subsequence.
[0017] In some embodiments, identifying the first subset of variants from the plurality of variants further includes, upon determining that the first length does not exceed a first specified threshold, determining whether the reference sequence construct comprises a first subsequence, where the first subsequence is identical to at least a first portion of the first structural variant, and, upon determining that the reference sequence construct comprises the first subsequence, excluding the first structural variant from the plurality of variants.
[0018] In some embodiments, determining whether the reference sequence construct comprises a first subsequence comprises determining whether the first subsequence has a length greater than a second specified threshold.
[0019] Some embodiments further include, upon determining that the reference sequence construct does not include the first subsequence, determining whether the second structural variant includes a second subsequence, where the second subsequence is identical to at least a second portion of the first structural variant, and, upon determining that the second structural variant includes the second subsequence, excluding either the first structural variant or the second structural variant from the plurality of variants.
[0020] In some embodiments, determining whether the second structural variant comprises a second subsequence comprises determining whether the second subsequence has a length greater than a second specified threshold.
[0021] In some embodiments, the second specified threshold is at least 150 base pairs.
[0022] In some embodiments, removing one of the first structural variant or the second structural variant from the plurality of variants comprises identifying a shortest variant from among the first structural variant and the second structural variant, and removing the shortest variant from the plurality of variants.
[0023] Some embodiments further include, if it is determined that the second structural variant does not include the second subsequence, determining whether a decoy sequence associated with the reference sequence construct includes a third subsequence, where the third subsequence is identical to at least a third portion of the first structural variant, and, if it is determined that the decoy sequence includes the third subsequence, masking the decoy sequence.
[0024] In some embodiments, identifying the filtered set of variants from the first subset of variants includes generating an initial graph reference construct using at least a portion of the first subset of variants.
[0025] In some embodiments, identifying the filtered set of variants from the first subset of variants further includes generating a plurality of graph reads using the initial graph reference construct, each of at least a portion of the plurality of graph reads being associated with a respective path in the initial graph reference construct.
[0026] In some embodiments, the plurality of graph reads includes a first subset of graph reads and a second subset of graph reads, and generating the plurality of graph reads includes generating the first subset of graph reads by traversing an initial graph reference construct over a first interval, and generating the second subset of graph reads by traversing the initial graph reference construct over a second interval, where the first interval and the second interval at least partially overlap.
[0027] In some embodiments, generating the multiple graph reads includes traversing the initial graph reference construct using a moving window with jumps.
[0028] Some embodiments further include aligning at least a portion of the plurality of graph reads to the initial graph reference construct, where the aligning includes, for each of the at least a portion of the plurality of graph reads, determining an alignment quality between the graph read and the graph reference construct, and determining whether the alignment quality exceeds a threshold.
[0029] Some embodiments further include identifying a first group of at least a portion of the plurality of graph reads, wherein each graph read included in the first group of at least a portion of the plurality of graph reads comprises a first combination of one or more variants of the first subset of variants.
[0030] In some embodiments, the first group of at least a portion of the plurality of graph reads includes a first graph read and a second graph read, and the method further includes excluding at least one multi-alignable variant from the filtered set of variants upon determining that neither the first alignment quality determined for the first graph read nor the second alignment quality determined for the second graph read exceeds a specified threshold.
[0031] In some embodiments, at least one multi-alignable variant is included within the first combination of one or more variants.
[0032] In some embodiments, identifying the filtered set of variants from the first subset of variants includes generating an initial graph reference construct using the first subset of variants, traversing the initial graph reference construct to generate a plurality of graph reads, aligning the plurality of graph reads to the initial graph reference construct, determining an alignment quality for each of at least a portion of the plurality of graph reads, and excluding at least a portion of one or more of the first set of variants from the second set of variants based on the alignment quality.
[0033] In some embodiments, one or more of the plurality of graph reads are associated with the same combination of one or more of the first subset of variants. Some embodiments further include determining whether each of the determined alignment qualities for one or more of the plurality of graph reads is below a specified threshold, and excluding at least one variant from the filtered set of variants upon determining that each of the alignment qualities is below a specified threshold.
[0034] In some embodiments, obtaining the plurality of variants includes obtaining a plurality of alternative sequences associated with a reference sequence construct and processing at least a portion of the plurality of alternative sequences, where the processing includes, for a first alternative sequence of the plurality of alternative sequences, aligning the first alternative sequence to the reference sequence construct and obtaining an aligned position, identifying one or more differences between the first alternative sequence and the reference sequence construct at the aligned position, and including at least a portion of the one or more differences in the plurality of variants as a first variant.
[0035] In some embodiments, at least a portion of the plurality of alternative sequences is processed to construct an updated reference sequence construct that does not include the plurality of alternative sequences.
[0036] In some embodiments, the first alternative sequence comprises an inverted sequence patch, and aligning the first alternative sequence to the reference sequence construct and obtaining an aligned position comprises obtaining an alternative aligned position for the inverted sequence patch.
[0037] Some embodiments further include left-normalizing the first variant to the reference sequence construct prior to including the first variant in the plurality of variants.
[0038] In some embodiments, at least a portion of the one or more differences include a first and a second consecutive difference, the first difference being associated with a first subsequence of the first alternative sequence and the second difference being associated with a second subsequence of the reference sequence construct. Some embodiments further include processing the first and second differences prior to including them in the plurality of variants as a first variant, the processing further including determining whether the first subsequence includes one or more regions included in the second subsequence, and removing the one or more regions from both the first and second subsequences upon determining that the first subsequence includes one or more regions included in the second subsequence.
[0039] In some embodiments, the first and second differences comprise an insertion and deletion event, respectively.
[0040] In some embodiments, obtaining the plurality of variants further comprises obtaining a second variant associated with the reference sequence construct, and including the second variant in the plurality of variants.
[0041] Some embodiments further include annotating the second variant with information indicating a source of the second variant.
[0042] In some embodiments, at least a portion of the first variants are associated with a first allele frequency, and at least a portion of the second variants are associated with a second allele frequency, some embodiments further include, for shared variants contained within both at least a portion of the first variants and at least a portion of the second variants, averaging the first and second allele frequencies associated with the shared variants to obtain an average allele frequency.
[0043] BRIEF DESCRIPTION OF THE DRAWINGS Various aspects and embodiments of the disclosure provided herein are described below with reference to the accompanying drawings, which are not intended to be drawn to scale. In the drawings, each identical or nearly identical component shown in various figures is represented by a similar reference numeral. For purposes of clarity, every component may not be labeled in every figure. [Brief description of the drawings]
[0044] [Figure 1] FIG. 1 is a diagram of an exemplary technique for generating a graph reference construct according to some embodiments of the technology described herein (SEQ ID NOs: 1-2). [Figure 2A]FIG. 2A is a flowchart of an example process 200 for generating a graph reference construct, according to some embodiments of the techniques described herein. [Figure 2B] FIG. 2B is a flow chart illustrating an exemplary process 220 for processing variants associated with a reference sequence construct, according to some embodiments of the technology described herein. [Figure 2C] FIG. 2C is a flow chart illustrating an exemplary process 240 for processing structural variants according to some embodiments of the technology described herein. [Figure 2D] FIG. 2D is a flow chart illustrating an example process 260 for identifying a filtered set of variants from among the first subset of variants according to some embodiments of the techniques described herein. [Figure 3A] FIG. 3A is an illustrative example of processing alternative sequences associated with a reference construct according to some embodiments of the technology described herein (SEQ ID NOs: 3-4). [Figure 3B] FIG. 3B is a diagram of an illustrative example of performing the first stage of a multi-stage variant filtering technique according to some embodiments of the technology described herein, the first stage being used to identify a set of structural variants to be excluded from an initial set of variants (SEQ ID NOs: 5-12). [Figure 3C] FIG. 3C is a diagram of an illustrative example of performing the second stage of a multi-stage variant filtering technique according to some embodiments of the technology described herein, the second stage being used to identify a set of multiple alignable variants to be excluded from the initial set of variants (SEQ ID NOs: 13-23). [Figure 4A] FIG. 4A illustrates an example process 400 for generating a graph reference construct, according to some embodiments of the techniques described herein. [Figure 4B]FIG. 4B shows an exemplary process 402 for processing alternative sequences associated with a reference sequence construct, according to some embodiments of the technology described herein. [Figure 4C] FIG. 4C illustrates an exemplary process 422 for identifying a set of structural variants according to some embodiments of the technology described herein. [Figure 4D] FIG. 4D illustrates an example process 424 for identifying a set of multiple alignable variants, according to some embodiments of the techniques described herein. [Diagram 5] FIG. 5 shows a graph illustrating alignment metrics from an experiment to measure the performance of graph reference construction, according to some embodiments of the techniques described herein. [Figure 6] FIG. 6 shows graphs illustrating variant calling metrics from an experiment to measure the performance of a graph reference construct, according to some embodiments of the technology described herein. [Figure 7A] FIG. 7A shows a graph showing cumulative variant number versus allele frequency from an experiment to measure the performance of a graph reference construction, according to some embodiments of the technology described herein. [Figure 7B] FIG. 7B shows a graph showing cumulative variant number versus allele frequency from an experiment to measure the performance of a graph reference construction, according to some embodiments of the technology described herein. [Figure 8A] FIG. 8A shows a graph showing the number of variants from an experiment to measure the performance of a graph reference construction, according to some embodiments of the technology described herein. [Figure 8B] FIG. 8B shows a graph showing the number of variants from an experiment to measure the performance of a graph reference construction, according to some embodiments of the technology described herein. [Figure 8C]FIG. 8C shows a graph showing the number of variants from an experiment to measure the performance of a graph reference construction, according to some embodiments of the technology described herein. [Figure 8D] FIG. 8D shows a graph showing the number of variants from an experiment to measure the performance of a graph reference construction, according to some embodiments of the technology described herein. [Figure 8E] FIG. 8E shows a graph showing the number of variants from an experiment to measure the performance of a graph reference construction, according to some embodiments of the technology described herein. [Figure 8F] FIG. 8F shows a graph showing the number of variants from an experiment to measure the performance of a graph reference construction, according to some embodiments of the technology described herein. [Figure 9] FIG. 9 is a block diagram of an example computer system that may be used to implement some embodiments of the techniques described herein. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS
[0045] Detailed Description Aligning sequence reads against a graph reference structure that describes known genetic variations between people helps to accurately place sequence reads and facilitates the identification of variants based on the alignment results. However, the inventors have recognized and understood that conventional techniques for aligning sequence reads against a graph reference structure can produce inaccurate results and are computationally expensive, and therefore can be improved.
[0046] When a graph reference construct includes all curated variants (e.g., variants selected to represent genetic variation) without considering how those variants may affect the alignment, aligning sequence reads to the graph reference construct may produce inaccurate results. First, the curated variants may include structural variants. Structural variants include insertions of at least a threshold length (e.g., at least 40 bp, at least 50 bp, at least 60 bp, at least 80 bp, at least 100 bp, at least 150 bp, at least 500 bp, at least 1 Kbp, at least 5 Kbp, at least 20 Kbp, at least 50 Kbp, at least 100 Kbp, at least 500 Kbp, etc.), deletions of at least a threshold length (e.g., at least 40 bp, at least 50 bp, at least 60 bp, at least 80 bp, at least 100 bp, at least 150 ... A sequence read may include an inversion of at least 80 bp, at least 100 bp, at least 150 bp, at least 500 bp, at least 1 Kbp, at least 5 Kbp, at least 20 Kbp, at least 50 Kbp, at least 100 Kbp, at least 500 Kbp, etc.), a duplication of at least a threshold length (e.g., at least 40 bp, at least 50 bp, at least 60 bp, at least 80 bp, at least 100 bp, at least 150 bp, at least 500 bp, at least 1 Kbp, at least 5 Kbp, at least 20 Kbp, at least 50 Kbp, at least 100 Kbp, at least 500 Kbp, etc.), and / or any other suitable structural variant. Structural variants may introduce ambiguity into the graph reference construction due to the nature of short read sequencing data. In other words, a sequence read may be misaligned to two or more positions within the graph reference construction if the structural variant includes a subsequence that is (a) identical to other portions of the graph reference, and (b) longer than the sequence read.Secondly, as more variants are incorporated into the graph reference construct, the number of possible paths in the graph increases exponentially, increasing the possibility that the same path will exist in different regions of the graph.As a result, sequence reads may be aligned to multiple regions in the graph reference construct, making them uninformative for variant calling.Such variants may be referred to herein as "multiply-alignable variants".
[0047] In addition, the curated variants may be obtained from multiple different sources, such as multiple variant databases or VCF files. As a result of discrepancies between variant representations in different bioinformatics pipelines, the same variant may be represented differently when obtained from different sources. The addition of such variants may introduce different, but ultimately equivalent, paths into the graph reference, resulting in alignment errors.
[0048] Furthermore, since curated variants may include many variants from many individuals, aligning sequence reads to such graph reference constructions may be computationally expensive.Since known variants in a graph reference may be represented by each path through the graph underlying the graph reference, increasing the number of known variants represented by the graph reference increases the number of paths through the graph that must be evaluated during alignment of sequence reads to the graph reference, which results in an increase in the computational complexity of performing alignment.Furthermore, the added complexity of the graph reference structure may introduce noise during alignment, reducing accuracy.
[0049] Therefore, the inventors have developed a technique for generating a graph reference construct that filters out variants that cause alignment ambiguity (e.g., structural variants and / or multiply alignable variants), resulting in more accurate alignment results, but also reducing the overall computational complexity of such alignment. In some embodiments, the set of variants may be filtered in multiple stages to identify variants to be included in the graph reference construct. For example, different filtering stages may include filtering out different types of variants (e.g., structural variants may be filtered out in one stage, and multiply alignable variants may be filtered out in another stage, e.g., a stage after the stage in which structural variants are filtered). In some embodiments, the identified variants may be used to construct a graph reference construct, for example, by adding nodes and edges representing the filtered set of variants to a linear reference construct.
[0050] Some embodiments provide computer-implemented techniques for generating a graph reference construct (e.g., a directed acyclic graph (DAG)). In some embodiments, the techniques include: (A) obtaining a plurality of variants associated with a reference sequence construct for at least one portion of a genome (e.g., at least one substantial portion, at least one chromosome, at least 10,000 nucleotides, etc.); (B) generating a graph reference construct using the plurality of variants and the reference sequence construct (e.g., an hg19 or hg38 genome reference); and (C) outputting the generated graph reference construct (e.g., storing the graph reference construct in memory such that it may then be used for various applications, including, for example, aligning sequence reads to the graph reference construct). In some embodiments, a technique for generating a graph reference construct includes: (A) filtering a plurality of variants to obtain a filtered set of variants, where the filtered set of variants is a subset of the plurality of variants, and where the filtering includes multiple filtering steps including a first filtering step (e.g., for excluding a first type of variant) and a second filtering step (e.g., for excluding a second type of variant) that is different from the first filtering step and is performed after the first filtering step; and (B) generating a graph reference construct using the filtered set of variants (the filtered set of variants by applying the first and second filtering steps) and a reference sequence construct.
[0051] In some embodiments, the first filtering step comprises identifying a first subset of variants from among the plurality of variants, at least in part, by excluding one or more structural variants (e.g., insertion, deletion, or inversion events) from the plurality of variants. In some embodiments, the second filtering step comprises identifying a filtered set of variants from among the first subset of variants (e.g., variants identified in the first filtering step), at least in part, by excluding one or more multi-alignable variants (e.g., variants that result in multi-mapping sequence reads) from the first subset of variants.
[0052] It should be understood that the techniques described herein may be implemented in any of numerous ways, as the techniques are not limited to any particular manner of implementation. Detailed examples of implementations are provided herein for illustrative purposes only. Furthermore, the techniques disclosed herein may be used individually or in any suitable combination, as aspects of the technology described herein are not limited to the use of any particular technique or combination of techniques.
[0053] Some exemplary aspects of the techniques described herein are described below with reference to FIGS.
[0054] FIG. 1 is a diagram of an example technique 100 for generating a graph reference construct, according to some embodiments of the techniques described herein. In some embodiments, the example technique 100 includes obtaining a plurality of variants 102. A first filtering stage 104 may be used to identify and remove one or more structural variants 106 from the plurality of variants 102, resulting in a first subset of variants 108. A second filtering stage 110 may be used to identify and remove one or more alignable variants 112 from the first subset of variants 108, resulting in a filtered set of variants 114. In some embodiments, the output of the second filtering stage 110 includes the filtered set of variants 114, a discarded set of variants 118 (e.g., removed during the first and second filtering stages), and / or a linear reference sequence construct 116. In some embodiments, the variants included in the filtered set of variants 114 and the linear reference sequence construct 116 are used to construct a graph reference sequence construct.
[0055] In some embodiments, obtaining the plurality of variants 102 includes obtaining variants from one or more sources. In some embodiments, this includes obtaining variants from one or more public variant databases and / or variant call format (VCF) files. For example, the plurality of variants may be obtained from GRCh38 human reference alternative contigs, 1000 Genomes Project Common Variants, Simons Genome Diversity Project Common Variants, Human Genome Structural Variant Consortium (HGSVC), and / or any other suitable variant databases and / or VCF files.
[0056] In some embodiments, the multiple variants 102 are associated with a reference sequence construct. For example, the reference sequence construct may include the GRCh38 genome assembly. In some embodiments, the reference sequence construct is constructed using a primary chromosome, a decoy, and a replacement sequence that represents deviations from the primary assembly. The decoy may include additional common sequences not in the reference. In some embodiments, if the decoy sequence is not included in the reference sequence construct, then sequence reads may erroneously map to regions of the primary chromosome. For example, HS38D1 and EBV decoys may be included in the reference sequence construct.
[0057] In some embodiments, the first filtering stage 104 includes identifying one or more structural variants 106 and excluding them from the plurality of variants to identify a first subset of variants. In some embodiments, the first filtering stage 104 includes evaluating the variants in multiple stages to determine whether including the variant in the graph construction would (a) make the computational cost for sequence alignment too high and / or (b) result in erroneous sequence alignment.
[0058] In some embodiments, including structural variants in a graph reference construct increases the computational complexity of aligning to such a graph reference construct. In some embodiments, the first filtering stage 104 includes filtering out structural variants that are too large. For example, insertions larger than a threshold size (e.g., larger than 1K, 2K, 3K, 5K, 10K, 15K, 20K, 25K, any number of base pairs in the range of 1-25K) may be filtered out from the plurality of variants. As another example, deletions larger than a threshold size (e.g., larger than 50K, 70K, 90K, 100K, 110K, 150K, 200K, 250K, 300K, any number of base pairs in the range of 50K-300K) may be filtered out in the first filtering stage. In some embodiments, the threshold size of different structural variants varies based on the characteristics of the aligner. In some embodiments, excluding these large structural variants from the multiple variants makes the alignment calculation feasible and significantly more computationally efficient, whereas without removing such structural variants, aligning sequence reads to the resulting graph would be computationally expensive or, in some cases, infeasible.
[0059] In some embodiments, (a) a structural variant that includes a subsequence that is identical to another subsequence included in the graph reference construct (e.g., another variant, a linear reference construct, or a decoy sequence) results in an inaccurate or ambiguous alignment. For example, if the length of the sequence read is shorter than such repeated subsequences, the sequence read may be aligned to each of those subsequences or may be misaligned to one of those subsequences. Thus, in some embodiments, the first filtering stage 104 includes determining whether the structural variant includes a subsequence that is identical to a subsequence included in the reference sequence construct, another variant included in the multiple variants, and / or a decoy sequence associated with the reference sequence construct. If a structural variant includes a subsequence that is identical to a subsequence included in the reference sequence construct and the subsequence is determined to have a length that exceeds a specified threshold (e.g., the length of the sequence read), the structural variant may be excluded from the multiple variants. If a structural variant includes a subsequence contained within another variant (e.g., another structural variant), and the subsequence is determined to have a length greater than a specified threshold, the shorter of the two variants may be excluded from the multiple variants. If a structural variant includes a subsequence contained within a decoy sequence, the subsequence is masked within the decoy sequence. In some embodiments, each of these determinations (e.g., regarding the reference sequence construct, other variants, and decoy sequence) may be made, or some of these determinations may be made, or only one of these determinations may be made. Aspects of using a first filtering stage to identify a first subset of variants are described herein, including at least the description of Figures 2C and 3B.
[0060] In some embodiments, the second filtering stage 110 includes identifying and excluding one or more multiply alignable variants 112 from the first subset of variants 108 to obtain a filtered set of variants 114. A "multiply alignable" variant may be a variant that, when incorporated into the graph reference construct, results in two or more identical paths in different, non-contiguous regions of the graph reference construct. For example, incorporating a multiply alignable variant into the graph reference construct may result in a first path in a first region of the graph reference construct that is identical to a second path in a second region of the graph reference construct, where the first path includes at least a portion (e.g., at least some or all) of the multiply alignable variants. Because a multiply alignable variant may result in two or more identical paths in the graph reference construct, a sequence read that aligns to one path in the graph reference construct may also align to at least one other path, the graph reference construct. Hence the name "multiple alignable" - such variants can cause sequence reads to align to multiple regions within the graph reference construct.
[0061] In some embodiments, the second filtering stage 110 includes evaluating whether inclusion of one or more variants in the graph reference construct would result in two or more identical paths in different (e.g., non-contiguous) regions of the graph reference construct (e.g., evaluating whether one or more variants are multiply alignable variants). In some embodiments, aligning sequence reads to a graph reference construct that includes identical paths in different regions (e.g., a graph reference construct that includes multiply alignable variants) may result in multiple mapping reads that are then uninformative for variant calling.
[0062] In some embodiments, the second filtering stage 110 includes generating a plurality of graph reads using the initial graph reference construct that includes the first subset 108 of variants. The graph reads may represent sequences in a particular region of the initial graph reference construct. One or more of the graph reads may then be aligned to the initial graph reference construct, respectively, to determine a mapping quality for each. The resulting mapping quality may indicate the likelihood that the alignment is correct. The mapping quality may then be used to identify multiple alignable variants. For example, when aligning the graph reads results in a low mapping quality (e.g., a mapping quality of 0), this may indicate that the graph reads align to multiple regions in the initial graph reference construct. In some embodiments, the multiple graph reads may represent the same variant or the same combination of variants. In this case, if aligning each of those graph reads results in a low mapping quality, then the shared variant or combination of variants is likely to result in one or more identical paths in the initial graph reference construct. As a result, a second filtering stage 110 may include excluding one or more of the shared variants (e.g., multiple alignable variants) 112 from the first subset of variants 108 to obtain a filtered set of variants 114. Aspects of identifying a filtered set of variants using the second filtering stage are described herein, including at least with respect to Figures 2D and 3C.
[0063] In some embodiments, the linear reference sequence construct 116 includes a linear human genome reference. For example, the linear reference sequence construct 116 may include a hg19 or hg38 human genome reference. In some embodiments, the linear reference sequence construct 116 may be subjected to one or more processing steps. For example, one or more alternative sequences may be removed from the linear reference sequence construct as described herein, including with respect to FIG. 2B. As another example, one or more of the discarded variants 118 (e.g., one or more of the multiple alignable variants 112) may be included as decoy sequences associated with the linear reference sequence construct 116. In some embodiments, the linear reference sequence construct 116 may be output as one or more files (e.g., one or more VCF files).
[0064] In some embodiments, generating the graph reference sequence construct 116 may include converting the linear reference construct 116 into a graph reference by adding nodes and edges that represent genetic variations. For example, the linear reference construct may be converted into a graph reference by adding nodes and edges that represent the filtered set of variants 114. Techniques for adding nodes and edges to a linear reference construct based on a set of variants are described in U.S. Patent Application Publication No. 2015-0057946, published February 26, 2015, entitled "METHODS AND SYSTEMS FOR ALIGNING SEQUENCES," which is incorporated by reference herein in its entirety.
[0065] FIG. 2A is a flowchart of an example process 200 for generating a graph reference construct, according to some embodiments of the techniques described herein.
[0066] In some embodiments, the process 200 begins at operation 202, where a plurality of variants associated with a reference sequence construct for at least one portion of a genome is obtained. In some embodiments, obtaining the plurality of variants includes accessing one or more variant databases and / or VCF files. For example, this may include accessing GRCh38 human reference alternative contigs, 1000 Genomes Project Common Variants, Simons Genome Diversity Project Common Variants, Human Genome Structural Variants Consortium (HGSVC), and / or any other suitable variants from any suitable variant database, data store, file, and / or VCF file. In some embodiments, the variants obtained from different databases and / or files may include variants from different population studies. In some embodiments, the different variant files may include the same variant or set of variants. Techniques for obtaining the plurality of variants are described herein, including at least with respect to FIG. 2B.
[0067] In some embodiments, the variants may be associated with a reference sequence construct, such as the GRCh38 genome assembly. In some embodiments, the reference sequence construct represents at least a portion of the genome. For example, the reference sequence construct may represent at least a substantial proportion of the genome (e.g., 80% of the genome), at least one chromosome, at least 10,000 nucleotides, or substantially the entire genome of a particular organism. In some embodiments, the variants associated with the linear reference construct are defined in relation to the reference sequence construct, as well as a coordinate system. For example, the variants may be represented by an identifier (e.g., a unique alphanumeric, alphabetic, or numeric character) that identifies the location of the variant relative to the reference sequence construct. Techniques for obtaining multiple variants are further described herein, including at least with reference to FIG. 2B.
[0068] After obtaining the plurality of variants, the process 200 proceeds to operation 204, where generating a graph reference construct using the plurality of variants and the reference sequence construct is performed. As described herein, in some embodiments, aligning the sequence reads to the graph reference construct including all the variants obtained in operation 202 may result in inaccurate or ambiguous alignments and may be computationally expensive. Therefore, as shown in FIG. 2A, operation 204 may include filtering the plurality of variants to obtain a filtered set of variants. In some embodiments, filtering the variants includes a first filtering step 206a and a second filtering step 206b.
[0069] In some embodiments, the first filtering stage 206a includes identifying a first subset of variants from the plurality of variants by excluding one or more structural variants from the plurality of variants. For example, the structural variants may include one or more insertions, deletions, inversions, duplications, or translocations of at least 50 bp in length. In some embodiments, identifying the first subset of variants includes identifying one or more structural variants for excluding from the plurality of variants. One example for processing one such structural variant is described herein, including at least the description of process 240 shown in FIG. 2C. In some embodiments, process 240 may be repeated to process a plurality.
[0070] In some embodiments, the second filtering step 206b includes identifying a second subset of variants from the plurality of variants by excluding one or more multi-alignable variants from the plurality of variants. For example, if the first subset of variants is included in a graph reference construct, the second filtering step may include determining whether a path in one region of the graph reference construct is identical to one or more paths in one or more other regions of the graph reference. In some embodiments, if identical paths are identified, the variants (e.g., multi-alignable variants) that give rise to such paths may be excluded from the graph to obtain a unique set of paths in the graph. An exemplary implementation of operation 206b is described herein, including with respect to FIG. 2D.
[0071] In some embodiments, after obtaining the filtered set of variants in operation 206, process 200 proceeds to operation 208, where generating a graph reference construct using the filtered set of variants is performed. In some embodiments, generating the graph reference construct may include adding one or more nodes or edges representing the filtered set of variants to a reference sequence construct.
[0072] In operation 210, the generated graph reference construct may be output. In some embodiments, outputting the graph reference construct may include storing the graph reference construct so that it can be subsequently used for one or more applications (e.g., to align sequence reads to the graph reference construct in any subsequent bioinformatics pipeline). For example, the generated graph reference construct may be stored locally on the computing device used to perform process 200 (e.g., on a non-transitory storage medium coupled to the computing device, or on a portion of the computing device). In some embodiments, the graph reference construct may be stored in one or more external storage media (e.g., a remote database or cloud storage environment, etc.). The stored graph reference construct may then be used, for example, to align sequence reads to the graph reference construct. FIG. 2B is a flow chart illustrating process 220 for processing variants associated with a reference sequence construct, according to some embodiments of the techniques described herein. Process 220 is an example of how operation 202 of process 200 may be implemented.
[0073] As shown, process 220 begins at operation 222 for obtaining a number of alternative sequences associated with a reference sequence construct for at least one portion of a genome. The alternative sequences, or alternative contigs, represent genetic deviations from the reference sequence construct (e.g., primary assembly). Thus, there are nucleotide sequence differences between the alternative sequences and the corresponding portions of the reference sequence construct. In some embodiments, the alternative sequences may be highly divergent (e.g., at least 80% novel) from the corresponding portions of the reference sequence construct. In some embodiments, the alternative sequences may be very similar (e.g., differing by only a few nucleotides) to the corresponding portions of the reference sequence construct.
[0074] In some embodiments, obtaining the alternative sequences in operation 222 includes obtaining one or more files describing an alignment of the alternative sequences to a reference sequence construct. For example, when using the GRCh38 assembly as the reference sequence construct, this may include obtaining one or more files from a general feature format (GFF) file describing an alignment of the alternative sequences to a primary chromosome. In some embodiments, the files describe the alignment of the alternative sequences in any suitable format. For example, the files may describe the alignment of the alternative sequences in a concise idiosyncratic gapped alignment report (CIGAR) format. However, it should be understood that the alternative sequences may be obtained from any suitable source (e.g., a database, a file, etc.) in any suitable format, as aspects of the technology described herein are not limited in this respect.
[0075] In some embodiments, the reference sequence construct includes alternative sequences as part of the primary assembly, and alternative sequences may be removed from the primary assembly, as aspects of the technology described herein include adding at least some of the processed alternative sequences to the reference sequence construct to obtain a graph reference construct.
[0076] As described above, some of the alternative sequences may be very similar to the reference sequence construct. Specifically, some of the alternative sequences may include large subsequences that are identical to subsequences included in the reference sequence construct. As a result, incorporating the alternative sequences into the graph reference construct may cause short sequence reads to erroneously align to multiple identical regions. Thus, process 220 includes techniques to address such concerns. Specifically, operation 224 includes processing at least some of the alternative sequences obtained in operation 222. In some embodiments, processing the alternative sequences includes suboperations 224a, 224b, and 224c.
[0077] As shown in FIG. 2B, sub-operation 224a includes aligning a first alternative sequence to a reference sequence construct and obtaining an alignment position for the first alternative sequence. In some embodiments, the alignment may be performed using any suitable alignment technique, as aspects of the technology described herein are not limited in this respect. For example, in some embodiments, the alignment may be performed using any of the techniques described in U.S. Patent Application Publication No. 2015-0057946, published February 26, 2015, entitled "METHODS AND SYSTEMS FOR ALIGNING SEQUENCES," which is incorporated by reference herein in its entirety. In some embodiments, an alignment position may have been previously obtained for the alternative sequence, making sub-operation 224a optional. For example, as described above, one or more files obtained in operation 222 may describe the alignment of the alternative sequence to the reference sequence construct.
[0078] After aligning the first alternative sequence in sub-operation 224a, process 220 proceeds to sub-operation 224b to identify one or more differences between the first alternative sequence and the reference sequence construct at the aligned positions. In some embodiments, the one or more differences include one or more nucleotide sequence differences. In some embodiments, the one or more differences can be sequence variants, such as substitutions, insertions, deletions, translocations, inversions, or any other suitable type of sequence mutation or variant. For example, the reference sequence construct can include the subsequence "AGGTCA", while the aligned alternative sequence includes the subsequence "AAGTCA". A "G" at the second position of the reference subsequence is substituted for an "A" at the second position of the alternative subsequence. In some embodiments, the one or more differences in sub-operation 224b can be identified using any suitable technique. For example, the technique can include processing one or more files describing the alignment in CIGAR (or any other) format and extracting the differences.
[0079] In some embodiments, the alternative sequence may include an inverted sequence patch. For example, the regions on both sides of the alternative sequence patch may be aligned in a forward direction with the reference sequence construct, while the inverted sequence patch is aligned in a reverse direction with the reference sequence construct. In some embodiments, the technique includes obtaining an alternative alignment for the inverted sequence patch, and then extracting one or more differences from the alternative alignment.
[0080] In some embodiments, the portion of the first alternative sequence that is not identified in operation 224b is excluded from further processing. For example, the portion of the first alternative sequence that is identical to the reference sequence construct may be excluded from further processing. In contrast, in some embodiments, one or more differences identified in operation 224b are included in further processing. This not only reduces the computational complexity (e.g., due to the large size of some alternative sequences before excluding identical portions), but also improves the accuracy of alignment of sequence reads. For example, if identical partial sequences are not excluded from the graph reference construct, sequence reads may be inaccurately aligned to both partial sequences.
[0081] After identifying one or more differences between the reference sequence construct and the first alternative sequence, the exemplary implementation proceeds to sub-operation 224c, which processes at least a portion of the one or more differences to obtain variants. In some embodiments, the differences may include consecutive differences, such as consecutive insertion and deletion events. Sometimes, the consecutive differences may include identical subsequences, which may be identified by aligning the differences to each other. Consider an exemplary insertion event including nucleotides "AGGTCGA" and an exemplary consecutive deletion event including nucleotides "CCGTCGG". After aligning the events to each other, the subsequence "GTCG" is identified as a matching subsequence (e.g., included in both the insertion and deletion events), e.g., using the Needleman-Wunsch algorithm. In some embodiments, sub-operation 224c includes removing the matching subsequence and splitting both differences into smaller mutations. In this example, removing the matching subsequence would result in insertions "AG" and "A", which would result in deletions of "CC" and "G". An example of processing the differences and removing matching subsequences is described herein, including at least with respect to Figure 3A. In some embodiments, processing at least some of the differences in sub-operation 224c further includes left-normalizing the differences to the reference sequence construct.
[0082] In some embodiments, as a result of operation 224c, one or more processed differences may be identified as first variants to be included in the multiple variants. In the example of FIG. 3A, the insertions "AG" and "A" and the deletions "CC" and "G" would be identified as first variants to be included in the multiple variants. In some embodiments, the first variants may be included in one or more input files in any suitable format. For example, the first variants may be included in one or more VCF files. As should be understood from the above, sub-operations 224a, 224b, and 224c may be performed for at least a portion of each of the multiple alternative sequences obtained in operation 222.
[0083] The process 220 then proceeds to operation 226, where a second variant associated with the reference sequence construct is obtained. In some embodiments, the second variant includes any of the variants described with respect to operation 202 of FIG. 2A, except for the alternative sequence obtained in operation 222.
[0084] Process 220 then proceeds to operation 228 where merging variants and obtaining multiple variants is performed. In some embodiments, the obtained multiple variants include multiple variants to be used as part of process 200 starting with operation 204 (as shown in FIG. 2A, FIG. 2B shows an exemplary implementation, the multiple variants output from operation 202 are provided as input to operation 204 and are filtered in operation 204). In some embodiments, merging variants includes processing input files describing the variants and consolidating variant structures for merging. In some embodiments, processing the input files includes splitting multi-allelic variants. In some embodiments, processing the input files may include removing non-standard variant definitions and leaving only fully resolved variants. In some embodiments, processing the input files may include additional filters, such as filtering by allele frequency and selecting second variants to be included. For example, in some embodiments, only variants having an allele frequency of at least a threshold percentage (e.g., at least 2%, at least 5%, at least 10%, at least 15%, etc.) may be included. In some embodiments, processing the input file may further include left-normalizing the variants. In some embodiments, processing the input file may include clearing unused annotations, clearing certain fields (e.g., ID and FILTER fields), and clearing sample information. In some embodiments, processing the input file may include annotating the file with information indicating the allele frequency. In some embodiments, processing the input file may include annotating the variants to indicate the source file (e.g., with an ID assigned to the file).
[0085] After processing the input files, the first and second variants may be merged. In some embodiments, merging the variants includes taking multiple input files and generating a single file (e.g., a VCF file, or a file in any other suitable format) that describes an initial graph reference that includes the first and second variants. In some embodiments, merging the input files may include aggregating annotations when the same variant comes from multiple sources. For example, a new effective allele frequency may be calculated for variants coming from multiple sources (e.g., having differential allele frequencies and different sample sizes). The final allele frequency may be determined by averaging the original allele frequencies weighted by the number of samples used for the corresponding source file.
[0086] In some embodiments, to perform operation 206a of process 200, in which identifying a first subset of variants is performed, one or more structural variants are identified for removal from the plurality of variants to obtain the first subset of variants. Figure 2C is a flowchart of an example process 240 for identifying a structural variant for removal from the plurality of variants. In some embodiments, process 240 can be repeated to identify one or more additional structural variants for removal from the plurality of variants.
[0087] As described herein above, the variants obtained in operation 202 of process 200 may include structural variants that are (a) large in size and / or (b) contain subsequences identical to subsequences contained elsewhere among the reference sequence construct, other variants, or decoy sequences. Thus, process 240 includes filtering such structural variants. An example of filtering structural variants is further described herein, including at least with respect to FIG. 3B.
[0088] In some embodiments, process 240 begins at operation 242, where a determination is made as to whether the length of a first structural variant exceeds a specified threshold. In some embodiments, different types of structural variants may be compared to different thresholds. For example, the length of an insertion may be compared to a first threshold (e.g., 2,500bp, 5,000bp, 7,500bp, 10,000bp, 20,000bp, etc.), while the length of a deletion may be compared to a second, different threshold (e.g., 50,000bp, 75,000bp, 90,000bp, 100,000bp, 150,000bp, 200,000bp, etc.). In other embodiments, different structural variants may be compared to the same threshold.
[0089] Regardless of the threshold, if the length of the first structural variant does exceed the specified threshold, the structural variant is removed from the plurality of variants in operation 254. If the length of the first structural variant does not exceed the threshold, then the exemplary implementation proceeds to operation 244, where a determination is made as to whether the reference sequence construct includes a subsequence that is identical to a portion of the first structural variant.
[0090] In some embodiments, determining whether the reference sequence construct includes a subsequence identical to the first portion of the first structural variant includes aligning the structural variant to the reference sequence construct. The reference sequence construct may be compared to the structural variant at the aligned position to determine whether they include any matching subsequences. If the reference sequence construct includes a subsequence identical to a subsequence included in the structural variant, a length is determined for the matching subsequence. In some embodiments, operation 244 includes determining whether the length of the matching subsequence is greater than a specified threshold. For example, the specified threshold may be similar to the length of the sequence read (e.g., 50bp, 100bp, 150bp, 200bp, 250bp, 300bp, etc.). In some embodiments, the specified threshold may vary based on the length of one or more sequence reads to be aligned. In some embodiments, if the matching subsequence is longer than the sequence read to be aligned to the graph reference construct, the sequence read may be incorrectly aligned to both or either of the subsequences (e.g., included in the structural variant and reference sequence construct) if the structural variant were to be included in the graph reference construct.
[0091] In some embodiments, if the reference sequence construct includes a subsequence that is identical to a portion (e.g., a subsequence) included in the first structural variant, and the length of the subsequence exceeds a specified threshold, then the first structural variant is removed from the plurality of variants in operation 254. If the reference sequence construct does not include a subsequence that is identical to a portion of the first structural variant and has a length that exceeds a specified threshold, then process 240 proceeds to operation 246.
[0092] Operation 246 may include determining whether the second structural variant includes a subsequence that is identical to a portion of the first structural variant. The determination may be made in any suitable manner, including, for example, aligning the first structural variant to one or more other variants. If the second structural variant includes a subsequence that is identical to a subsequence included in the first structural variant, a length is determined for the matching subsequence. For example, the designated threshold may be similar to the length of the sequence read (e.g., 50 bp, 100 bp, 150 bp, 200 bp, 250 bp, 300 bp, etc.). In some embodiments, the designated threshold may vary based on the length of one or more sequence reads to be aligned. In some embodiments, the threshold may be the same threshold used in operation 244. In some embodiments, the threshold may be different from the threshold used in operation 244. In some embodiments, if a matching subsequence is longer than the sequence read to be aligned to the graph reference construct, the sequence read may be inaccurately aligned to both or either subsequence (e.g., included in the first structural variant and second sequence construct) if both the first and second structural variants are included in the graph reference construct.
[0093] In some embodiments, if the second structural variant includes a subsequence that is identical to a portion of the first structural variant, and the length of the subsequence exceeds a specified threshold, then the process 240 proceeds to operation 252. Operation 252 may include determining which structural variants should be excluded. In some embodiments, it may be desirable to exclude shorter ones of the structural variants, since longer structural variants contain more information. Therefore, operation 252 may include determining whether the length of the second structural variant exceeds the length of the first structural variant. If it is determined that the length of the second structural variant exceeds the length of the first structural variant, then in operation 254, the first structural variant is excluded from the plurality of variants. If it is determined that the length of the second structural variant does not exceed the length of the first structural variant, then in operation 256, the second structural variant is excluded from the plurality of variants.
[0094] If it is determined in operation 246 that the second structural variant does not include a subsequence that is identical to a portion of the first structural variant and has a length that exceeds a specified threshold, the process 240 proceeds to operation 248. Operation 248 may include determining whether the decoy sequence includes a subsequence that is identical to a portion of the first structural variant. As described herein, the decoy sequence may include a common sequence that is not included in the reference. However, if one of the common sequences is already represented by the structural variant, then it is not necessary to include that sequence as a decoy. Thus, if the decoy sequence includes a subsequence that is identical to a subsequence included in the first structural variant, in operation 258, that region of the decoy sequence is masked. The process 240 then proceeds to operation 250, where the first structural variant is included in a first partial set of variants.
[0095] 2D is a flow chart illustrating a process 260 for identifying a filtered set of variants from among the first subset of variants, according to some embodiments of the techniques described herein. Process 260 is an example of how operation 206b of process 200 may be implemented.
[0096] As described herein above, the more variants are included in the graph reference construct, the more likely it is that the same path is included in different regions of the graph reference construct. Aligning sequence reads to such a graph reference construct may result in ambiguous and uninformative results due to multiple mapping sequence reads. In some embodiments, the alignment quality indicates the degree of certainty that the alignment is correct. If there are multiple regions in the graph to which the sequence reads are mapped (e.g., multiple mapping), the mapping quality may be low. In some embodiments, a filtering step, such as the exemplary implementation 206b, may be used to exclude some variants that result in multiple mapping sequence reads (e.g., multiple alignable variants) to break the identity of different regions in the graph reference construct. An example of filtering multiple alignable variants is described herein, including at least the description of FIG. 3C.
[0097] In some embodiments, the example implementation 206b starts at 262, where an initial graph reference construct is generated using the reference sequence construct and at least some variants of the first subset of variants identified in operation 250 of process 240. In some embodiments, at least some variants in the first subset of variants may be added to the reference sequence construct using one or more nodes and / or edges to generate the initial graph reference construct. Thus, the initial graph reference construct may include one path representing the reference sequence construct and one or more paths representing variants included in the initial graph reference construct. An "edge combination" may refer to a path in the initial graph reference construct that follows one or more particular edges and thus represents one or more variants associated with those edges (e.g., variants are included as edges, included as nodes following edges, etc.).
[0098] The example implementation 206b then proceeds to operation 264, where it traverses the initial graph reference construct and synthetically generates a number of graph reads of a specified length from the graph reference construct. A graph read may include one or more nucleotides that represent a path in a particular region in the initial graph reference construct. In some embodiments, graph reads are generated for all possible haplotypes in the graph. In some embodiments, traversing the initial graph reference construct to generate graph reads includes traversing the graph reference construct using a moving window with jumps. In some embodiments, operation 264 includes performing sub-operations 264a and 264b.
[0099] Sub-operation 264a, in some embodiments, includes generating a first subset of graph reads by traversing the graph reference construct over a first interval. In some embodiments, the first subset of graph reads may include one reference graph read and one or more non-reference graph reads. The reference graph reads may represent paths through the reference sequence construct, while the non-reference graph reads may represent paths that follow edges (e.g., combinations of edges) in the initial graph reference construct within the interval.
[0100] Sub-operation 264b may include generating a second subset of graph reads by traversing the initial graph reference construct over a second interval that overlaps with the first interval. As described herein above, the second subset of graph reads may include one reference graph read and one or more non-reference graph reads. In some embodiments, because the first and second intervals overlap, the graph reads included in the second subset of graph reads may represent one or more variants represented by the graph reads included in the first subset of graph reads.
[0101] In some embodiments, after generating the multiple graph reads in operation 264, the example implementation 206b proceeds to operation 266, where it aligns the multiple graph reads to the initial graph reference construct and determines an alignment quality for each of at least some of the multiple graph reads. As described herein above, the alignment quality may indicate the likelihood that the graph read is correctly aligned to the initial graph reference construct. Determining the alignment quality for the graph read may include determining whether the graph read maps to more than one region in the initial graph reference construct. In some embodiments, a graph read that maps to more than one region in the initial graph reference construct results in a lower alignment quality than a graph read that maps to only one location in the graph reference construct. This is because a graph read that maps to only one location in the initial graph reference construct is more likely to be mapped to the correct location than a graph read that may map to multiple locations.
[0102] In some embodiments, for a subset of graph reads (e.g., a first subset of graph reads, or a second subset of graph reads), the alignment quality determined for the reference graph read is compared to the alignment quality determined for the non-reference graph read. In some embodiments, if the non-reference graph read has an alignment quality lower than the alignment quality determined for the reference graph read, then a combination of edges represented by the non-reference graph read may be identified for removal from the filtered set of variants. For example, if the alignment quality of the non-reference graph read is 0, while the alignment quality of the reference graph read is greater than 0, then a combination of edges represented by the non-reference graph read may be identified for removal from the filtered set of variants. In some embodiments, if the alignment quality determined for the non-reference graph read is greater than the alignment quality determined for the reference graph read, then a combination of edges represented by the non-reference graph read may be identified for inclusion in the filtered set of variants. Additionally or alternatively, edge combinations represented by non-reference graph reads may be identified for inclusion in the filtered set of variants if the non-reference graph reads have an alignment quality greater than a specified threshold (e.g., at least 10, at least 20, at least 30, at least 40, etc.). However, it should be understood that while edge combinations may be identified for inclusion or exclusion from the filtered set of variants, in some embodiments, edge combinations may not actually be included in or excluded from the filtered set of variants in operation 268 of process 200.
[0103] After determining the alignment quality for at least some of the graph reads, the example implementation 206b proceeds to operation 268, where it removes at least some of the first subset of variants from the filtered set of variants. In some embodiments, in operation 268, non-reference graph reads that include the same combination of edges may be grouped together. For example, because the first and second intervals overlap, a non-reference graph read included in the first subset may represent a combination of edges that are similarly represented by a non-reference graph read included in the second subset. Therefore, those graph reads may be grouped together.
[0104] In some embodiments, if each of the grouped non-reference graph reads is identified for removal from the filtered set of variants in operation 266, this may indicate that the edge combination results in multiple mapping sequence reads (e.g., reads that align to multiple different regions of the graph). Thus, the edge combination may be identified for filtration. In some embodiments, in operation 268, a set of variants represented by the edge combination identified for filtration is removed from the filtered set of variants. For example, a set of variants is identified from the edge combination identified for filtration such that each edge combination has at least one variant removed from the filtered set of variants.
[0105] 3A is a diagram of an illustrative example of processing alternative sequences associated with a reference construct, according to some embodiments of the technology described herein. The example of FIG. 3A serves as an example of performing operation 224 of process 220.
[0106] In some embodiments, the example 300 begins with operation 302, where an alignment of the alternative sequence to the reference sequence construct is performed. As part of the alignment, one or more differences are identified in the aligned positions and are represented by shaded boxes. The example includes annotations identifying regions of match and regions containing structural variants, along with the number of nucleotides contained within the region. In some embodiments, a region annotated with "M" indicates a matching nucleotide. For example, a region annotated with "M3" represents three matching nucleotides, while a region annotated with "M23" represents 23 matching nucleotides. In some embodiments, a matching region may include one or more mismatches. For example, the region "M23" includes two mismatches. First, at position 19, the nucleotide "G" in the reference sequence construct does not match the nucleotide "T" in the alternative sequence. Second, at position 30, the nucleotide "G" in the reference sequence construct does not match the nucleotide "T" in the alternative sequence. In some embodiments, a region annotated with "I" indicates an insertion. For example, region "I5" represents an insertion of five nucleotides. As shown in the alternate sequence, the five shaded boxes represent an insertion of the nucleotides "GACCG". As another example, region "I4" represents an insertion of four nucleotides. As shown in the alternate sequence, the four shaded boxes represent an insertion of the nucleotides "AGTT". In some embodiments, a region annotated with a "D" indicates a deletion. For example, region "D4" represents a deletion of four nucleotides. As shown in the reference sequence construct, the four shaded boxes represent a deletion of the nucleotides "TACC". As another example, region "D3" represents a deletion of three nucleotides. As shown in the reference sequence construct, the three shaded boxes represent a deletion of the nucleotides "AAT".
[0107] In some embodiments, in operation 304, some of the differences identified in operation 302 may be processed. In some embodiments, operation 304 may include splitting complex variants, such as insertion and deletion events, to generate smaller variants. For example, in operation 304, consecutive insertion and deletion events, represented by regions "I5" and "D4", may be processed. As shown, the insertion and deletion events are aligned with respect to each other to determine whether they contain any matching nucleotides. The aligned positions include a matching region "M4" and an insertion region "I1". The matching region includes one mismatch and three matches, as represented by the grey boxes. Thus, complex variants (e.g., insertion and deletion events) can be split into smaller variants. As shown, the mismatched nucleotides in region "M4" may be represented as single nucleotide polymorphisms (SNPs), while the insertions in region "I1" may be represented by single nucleotide insertions. The matching regions are removed, which (a) simplifies the variant and (b) reduces the size of the variant.
[0108] In operation 306, a first variant is obtained that represents a simplified version of the alternative sequence. As shown, the variant is left normalized, i.e., the starting position of the variant is shifted to the left. In some embodiments, the first variant may be annotated to indicate the starting position relative to the reference sequence construct. For example, a first variant annotated with the number "4" indicates that it starts at the fourth position (e.g., the fourth nucleotide) from the left of the reference sequence construct.
[0109] In some embodiments, a VCF file may be output that includes the first variant obtained in operation 306. The VCF file may include the position of the variant and the nucleotides that define the variant relative to the reference sequence construct and the alternative sequence. For example, at position 13, the alternative sequence includes the nucleotide "C", which matches the first nucleotide in the sequence "CAAT" of the reference sequence construct. The nucleotide "AAT" represents a deletion event following the nucleotide at position 13 in the alternative sequence. Thus, the reference sequence includes the nucleotide "AAT" following position 13, but does not include the alternative sequence.
[0110] 3B is a diagram of an illustrative example of performing a first stage of a multi-stage variant filtering technique, according to some embodiments of the techniques described herein, where the first stage is used to identify a set of structural variants to be excluded from an initial set of variants. The example of FIG. 3B serves as an example of performing operation 206a of process 200, as described herein, including at least the description with respect to FIG. 2A, to identify one or more structural variants for exclusion from a plurality of variants.
[0111] In this example, identifying the set of structural variants includes four steps 322, 324, 326, and 328. However, it should be understood that in some embodiments, one or more steps may be omitted. For example, if the decoy sequence is not associated with the reference sequence construct, then step 328 may be omitted. Such omission will not affect the performance of the remaining three steps 322, 324, and 326.
[0112] In some embodiments, the first stage 322 includes aligning the structural variants, such as an insertion, to a reference sequence construct. As shown in FIG. 3B, two structural variants are aligned to the reference sequence construct to determine two alignments, alignment 332 and alignment 334.
[0113] In some embodiments, the first structural variant is compared with the reference sequence construct at the alignment position, and whether the first structural variant comprises a subsequence that is identical to the subsequence contained in the reference sequence construct and has a length greater than a specified threshold. In other words, this may include determining the length of the matching region at the alignment position. For example, when the first structural variant is aligned with the reference sequence construct and alignment 332 is determined, it comprises three matching regions. The first matching region has a length of 8 nucleotides, the second matching region has a length of 42 nucleotides, and the third matching region has a length of 19 nucleotides. Compared with the exemplary threshold of 30 nucleotides, the length of the second matching region (e.g., 42 nucleotides) exceeds the threshold. Therefore, the first structural variant will be excluded from the multiple variants.
[0114] In some embodiments, if the first structural variant was included in the multiple variants and used to generate the graph reference construct rather than being filtered out in operation 322, it may have resulted in ambiguous sequence read alignments. For example, a sequence read having a length less than 42 nucleotides (e.g., 30 nucleotides) may align to both the reference sequence construct and the first structural variant within the matching region. In this case, there would be no way to determine which alignment is correct, and as a result, the alignment would be uninformative.
[0115] As another example, when the second structural variant is aligned to the reference sequence construct and alignment 334 is determined, it contains four matching regions. The first matching region has a length of 8 nucleotides, the second matching region has a length of 20 nucleotides, the third matching region has a length of 18 nucleotides, and the fourth matching region has a length of 19 nucleotides. Since none of the matching regions has a length exceeding the exemplary threshold of 30 nucleotides, the second structural variant is not excluded from the multiple variants.
[0116] In some embodiments, the second stage 324 includes filtering the structural variants based on their size. For example, if the length of a deletion event is greater than a maximum deletion size threshold (e.g., 90,000 bp), then the deletion event may be excluded from the plurality of variants. Similarly, if the length of an insertion event is greater than a maximum insertion size threshold (e.g., 5,000 bp), then the insertion event may be excluded from the plurality of variants. If the length of the insertion or deletion event does not exceed the maximum size threshold, then those structural variants may be included in the plurality of variants or may be subjected to further filtering stages 322, 326, 328. In some embodiments, the structural variants excluded from the plurality of variants may be included as additional decoy sequences.
[0117] In some embodiments, the third stage 326 includes determining whether the two structural variants contain identical subsequences whose length exceeds a specified threshold. As shown in the first alignment 338, there are two matching regions (e.g., two identical subsequences). The first matching region has a length of 8 nucleotides, while the second matching region has a length of 51 nucleotides. Since the length of the second matching region (e.g., 51 nucleotides) exceeds an exemplary threshold of 30 nucleotides, one of the structural variants is excluded from the multiple variants. Since the longer structural variant contains more information, the shorter structural variant is excluded from the multiple variants.
[0118] As another example of step 326, alignment 340 shows an alignment of a different pair of structural variants. As shown, alignment 340 includes three matching regions. The first matching region has a length of 6 nucleotides, the second matching region has a length of 22 nucleotides, and the third matching region has a length of 326 nucleotides. Since none of the matching regions has a length exceeding the exemplary threshold of 30 nucleotides, neither structural variant is excluded from the multiple variants.
[0119] In some embodiments, the filtering step 328 includes aligning the structural variants to the decoy sequence to obtain an aligned position 342. Since the sequence represented by the structural variant will be included in the graph reference construct, there is no reason to additionally include it in the decoy sequence. Furthermore, including the sequence in the decoy sequence will result in the sequence reads aligning to both the decoy sequence and the structural variant representing that sequence. Therefore, the decoy sequence is masked at the aligned position to obtain a masked decoy sequence 344. In some embodiments, the structural variants may not align to the decoy sequence in the filtering step 328. Therefore, the region of the decoy sequence will not be masked.
[0120] In some embodiments, illustrative example 320 is performed as part of operation 206a, and in this example, structural variants that are not excluded from the plurality of variants would have been included as part of the first subset of variants created in process 206a.
[0121] 3C is a diagram of an illustrative example of performing a second stage of a multi-stage variant filtering technique, according to some embodiments of the techniques described herein, where the second stage is used to identify a set of multi-alignable variants to be excluded from the initial set of variants. The example of FIG. 3C serves as an example of performing operation 206b of process 200, which identifies one or more multi-alignable variants for exclusion from the multiple variants.
[0122] In some embodiments, an initial graph reference construct 362 may be generated. In some embodiments, this may include adding a first subset of variants to the reference sequence construct, obtained as a result of using a first filtering stage (e.g., a first filtering stage described herein, including at least the description with respect to Figures 2A and 3B). As shown in the example, the initial graph reference construct includes a variant at position 12, a variant at position 16, and a variant starting at position 36. The variants are represented in the graph using nodes and edges.
[0123] In some embodiments, the first stage 352 includes generating multiple graph reads by traversing the initial graph reference 362 over a specified interval. As shown, a first subset 364 of graph reads is generated for a first interval in the graph. The first subset 364 of graph reads includes one graph read that represents a path through the reference sequence construct represented by a graph read that includes only open squares. The remaining graph reads included in the first subset 364 of graph reads represent paths that include different combinations of edges in the graph. For example, one graph read represents a path that follows along the edge at position 12, including the variant represented at that position. Another graph read represents a path that follows along the edge at position 16. The final graph read represents a path that follows along both edges, the edge at position 12 and the edge at position 16.
[0124] A second subset of graph reads 366 is generated by traversing the initial graph reference construct over a second interval within the initial graph reference construct that overlaps with the first interval. Similarly, the second subset of graph reads 366 includes a reference graph read and three non-reference graph reads that represent three different edge combinations (e.g., an edge at position 12, an edge at position 16, and an edge at positions 12 and 16). As shown, the overlapping intervals result in graph reads that include the same edge combinations.
[0125] Finally, a third subset of graph reads 368 is generated by traversing the initial graph reference construct over a third interval within the initial graph reference construct that overlaps with the second interval. The third interval only includes one combination of edges as represented by the variant included at position 36. Thus, the third subset of graph reads 368 includes one reference graph read and one non-reference graph read.
[0126] In some embodiments, the generated graph reads may be collected as a FASTQ file. As shown in step 354, the FASTQ file and the initial graph reference construct 362 may be used to align the graph reads to the initial graph reference construct using a graph aligner to obtain a BAM file that is used to represent the aligned sequences. In some embodiments, the BAM file may include an alignment quality, or mapping quality ("MQ"), for each of the graph reads. The alignment quality may indicate the likelihood of correct alignment. As described herein above, including the description of FIG. 2D, a graph read with low alignment quality may indicate that the graph read may align to more than one location in the initial graph reference construct.
[0127] As shown in step 356, each graph read is annotated with an alignment quality ("MQ"). If the alignment quality of the non-reference graph read is less than the alignment quality of the reference graph read, the edge combination represented by the non-reference graph read is identified for removal from the filtered set of graph reads (labeled as "bad" in FIG. 3C). If the alignment quality of the non-reference graph read is greater than the alignment quality of the reference graph read and / or is greater than a specified threshold, the edge combination represented by the graph read is identified for inclusion in the filtered set of graph reads (labeled as "good" in FIG. 3C). Otherwise, the edge combination and associated graph read are ignored.
[0128] As shown in the example, a first subset of graph reads includes a reference graph read with an alignment quality of 25. Each of the non-reference graph reads has an alignment quality less than 25 (e.g., 0), so the combination of edges represented by the non-reference graph reads is identified for exclusion from the filtered set of graph reads. A second subset of graph reads 374 includes a reference graph read with an alignment quality of 35. Two of the non-reference graph reads have an alignment quality less than 35, so the combination of edges represented by those graph reads is identified for exclusion from the filtered set of graph reads. One non-reference graph read included in the second subset 374 has an alignment quality of 45, which is greater than the alignment quality of the reference graph read. Thus, the combination of edges represented by that graph read is identified for inclusion in the filtered set of graph reads. Finally, a third subset of graph reads 376 is ignored, because both the reference and non-reference graph reads in this subset have the same mapping quality of 0.
[0129] After classification, the graph reads are grouped by the combination of edges they represent. For example, a first group 378 represents the combination of edges that contain variant "G" at position 16, a second group 380 represents the combination of edges that contain variant "T" at position 12, and a third group 382 represents the combination of edges that contain both variants "T" and "G" at positions 12 and 16, respectively.
[0130] Next, each group 378, 380, 382 is classified based on the classification of the graph reads contained within the group. For example, group 378 includes all graph reads identified for exclusion from the filtered set of variants. This indicates that a combination of edges including variant "G" may result in the same path in different regions within the initial graph reference construct 362, resulting in multiple mapping sequence reads. Thus, group 378 is identified for filtration. Group 380 includes mixed classifications (e.g., graph reads identified for both inclusion and exclusion from the filtered set of variants). Thus, group 380 is not identified for filtration. Finally, group 382 includes all graph reads identified for exclusion from the filtered set of variants. Thus, group 382 is identified for filtration.
[0131] After sorting the groups, a set of variants is excluded from the first subset of variants among the variants included in the groups identified for filtration. Identifying the set of variants may include, in some embodiments, identifying one or more variants that are common to the groups identified for filtration. For example, as shown in FIG. 3C, the variant at position 16 is identified for exclusion because it is included in both groups 378, 382 identified for filtering. Thus, the variant is excluded from the first subset of variants.
[0132] In some embodiments, illustrative example 350 is performed as part of operation 206b, and in this example, variants that are not excluded from the plurality of variants would have been included as part of the filtered ones of the variants created in process 206b.
[0133] In some embodiments, the filtered set of variants is used to generate a graph reference construct at step 384 .
[0134] Further aspects of graph construction Additional aspects of the graph construction techniques described herein are described below with reference to FIGS.
[0135] 4A illustrates an example process 400 for generating a graph reference construct, according to some embodiments of the techniques described herein. In some embodiments, the process 400 is an example implementation of the example technique 100 and process 200 for generating a graph reference construct, described herein, including at least with respect to FIGS. 1 and 2A.
[0136] In some embodiments, the process 400 includes an operation 408 of processing the linear reference construct. In some embodiments, prior to operation 408, the process 400 includes obtaining the linear reference construct 404 and a decoy 406 associated with the linear reference construct 404. For example, as shown in FIG. 4B, the linear reference construct 404 in this example is the GRCh38 genome assembly. In some embodiments, the GRCh38 genome assembly includes primary chromosomes 432, unplaced and unlocalized contigs 434, alternative (ALT) and NOVEL contigs 436, and FIX contigs 438.
[0137] In some embodiments, the ALT and NOVEL contigs 436 represent alternative sequences for certain regions in canonical chromosomes. These regions show high variability within the population, and the ALT and NOVEL contigs 436 are provided as additional sequences to extend the haploid genome. In some embodiments, the ALT and NOVEL contigs 436 are obtained as General Feature Format (GFF) files. The GFF files describe the alignment of the alternative contigs to the canonical regions in the Concise Idiosyncratic Gapped Alignment Report (CIGAR) format. However, it should be understood that the data may be formatted in any other suitable format, as aspects of the technology described herein are not limited in this respect. In some embodiments, the ALT contigs are examples of alternative sequences, as described herein, including at least with respect to Figures 1A-3C.
[0138] In some embodiments, operation 408 includes processing the linear reference construct 404 to remove ALT and NOVEL contigs 436 from the linear reference 404 so that it contains only primary chromosomes and unplaced and unlocated contigs 434. Additionally, or alternatively, in operation 408, decoys 406 may be added to the linear reference construct to obtain a linear reference construct 412. The linear reference construct 412 may be output as a FASTA file, or data in any other suitable format.
[0139] In some embodiments, ALT and NOVEL contigs 436 are removed from the linear reference 404 before mapping them to the primary chromosome 432. Because ALT and NOVEL contigs 436 often contain long sequences that are identical to the linear reference, further processing is performed in operation 408 to (a) decompose the ALT and NOVEL contigs 436 into smaller variants and (b) left-normalize the decomposed variants. The resulting variants 410 may be output as a FASTA file, or data in any other suitable format. Operation 408 is an example of the type of processing that may be performed in operation 224 of process 220, in which processing the alternative sequences and obtaining second variants associated with the reference sequence construct is performed, as described herein, including at least the description with respect to FIG. 2B.
[0140] As an example of decomposing ALT and NOVEL contigs 436, consecutive insertion and deletion events can be combined and the alignment can be simplified (e.g., simplified to a larger number of SNPs). In some embodiments, the mutations are aligned to each other, for example, using the Needleman-Wunsch algorithm, to obtain a more minimal representation of the mutation. If a long sequence of the same match block is identified in the alignment, then the mutation can be split into smaller mutations from the match block.
[0141] In some embodiments, process 400 includes operation 416, where an input 414 is obtained and prepared. In some embodiments, the input 414 includes a variant file (e.g., a VCF file). In some embodiments, the input 414 is obtained from one or more sources. If the input 414 is obtained from different sources, preparing the input in operation 416 includes processing the input 414 and consolidating the variant structures. For example, processing the input 414 may include splitting multi-allelic variants, removing non-standard variant definitions and leaving only fully sequence resolved variants, filtering by allele frequency, left-normalizing the variants, clearing unused annotations, clearing ID and FILTER fields, clearing sample information, annotating with information used to calculate effective allele frequency, and / or annotating the variants to point to their original source file using an ID assigned to the respective VCF file.
[0142] Operations 408 and 414 may be performed in any order. In some embodiments, operations 408 and 414 may be performed simultaneously.
[0143] After obtaining and preparing the input 414 in operation 416 and processing the linear reference construct 404 and the decoys 406 in operation 408, the process 400 proceeds to operation 418, which is merging variants. Operation 418 is an example of the type of operation that may be performed in operation 228 of the process 220, in which merging the first and second variants and obtaining multiple variants is performed, as described herein, including at least the description with respect to FIG. 2B. In some embodiments, merging variants in operation 418 includes processing the multiple input variant files and obtaining a single bi-allelic candidate graph file. For example, operation 418 may include processing the prepared input 414 and the alternative variants 410 and obtaining a single output file that describes the initial graph reference construct. In some embodiments, merging includes combining all variants into a single set. If the same variant comes from multiple sources, then merging in operation 418 includes aggregating annotations associated with the variants and calculating an effective allele frequency for the variant. In some embodiments, the effective allele frequency for a variant is the average of the allele frequencies from all of the source files, weighted by the number of samples used for the corresponding source file.
[0144] Process 400 then proceeds to operation 420, where the variants (e.g., the variants output in operation 418) are filtered. In some embodiments, in operation 420, a multi-stage filtering technique may be used to filter the variants to identify a set 430 of variants to be excluded from the graph reference constructs 426, 428 obtained as an output of process 400. Operation 420 is an example of the type of processing that may be performed in operation 206 of process 200, where filtering of a plurality of variants associated with a reference sequence construct is performed to obtain a filtered set of variants, as described herein, including at least as described with respect to FIG. 2A.
[0145] In some embodiments, the filtering in operation 420 includes a structural variant (SV) filter 422 and a multiple map filter 424. In some embodiments, the SV filter 422 is a type of filter that can be used as part of the first filtering stage 206a of the process 200, described herein, including at least in the description with respect to FIG. 2A. In some embodiments, the SV filter 422 can be used to identify structural variants that should be excluded from the graph reference constructs 426, 428. This can eliminate introducing sequences that would result in duplications in the graph reference constructs. One example of using the SV filter 422 is described herein, including at least in the description with respect to FIG. 4C.
[0146] In some embodiments, the multiple map filter 424 may be used as part of the second filtering stage 206b of the process 200, described herein, including at least in the discussion with respect to Figure 2A. In some embodiments, the multiple map filter 424 may be used to identify multiply alignable variants that, if included in the graph reference construct 426, 428, would cause sequence reads to align to multiple regions of the graph reference construct 426, 428 (e.g., resulting in a multiple mapping problem). In some embodiments, the identified variants are excluded from the graph reference construct 426, 428. One example of using the multiple map filter 424 is described herein, including at least in the discussion with respect to Figure 4D.
[0147] In some embodiments, filtered variants 430 and graph reference constructs 426, 428 are obtained as output of the filtering in operation 420. In some embodiments, filtered variants 430 include those variants identified for removal using the SV filter 422 and the multi-map filter 424. Filtered variants 430 may be output as a VCF file or as data in any other suitable format.
[0148] In some embodiments, the graph reference constructs 426, 428 are obtained by aligning variants excluded from the filtered variants 430 to the linear reference construct 412 output in operation 408. For example, the variants include variants output in operation 418 that were not identified for exclusion in operation 420. In some embodiments, the graph reference constructs are output as a FASTA file 426, as a VCF file 428, and / or as data in any suitable format.
[0149] 4C illustrates an example process 422 for identifying a set of structural variants according to some embodiments of the techniques described herein. In some embodiments, the process 422 includes processing the structural variants represented in an initial graph reference construct 442 (e.g., output in operation 418, where merging is performed).
[0150] In some embodiments, operation 444 includes filtering the structural variants by size. For example, structural variants having a size above a threshold may be excluded from the filtered graph 454 and included in the filtered structural variants 452. In some embodiments, insertions are compared to a different threshold than deletions, or insertions and deletions are compared to the same threshold.
[0151] In some embodiments, operation 446 includes aligning the structural variants that were not filtered out in operation 444 to the linear reference construct 412. The structural variants may be aligned using any suitable alignment technique, such as the Minimap2 method described by Heng Li ("Minimap2: pairwise alignment for nucleotide sequences", Bioinformatics, Vol. 34, Issue 18, 2018, pp. 3094-3100), which is incorporated by reference in its entirety herein. If there is a subsequence (e.g., a match block) in the structural variant that is identical to a non-decoy sequence in the linear reference and is at least the length of the sequence read (e.g., 150 bp), the structural variant is excluded from the filtered graph 454 and included in the filtered structural variants 452. In some embodiments, the threshold is chosen to be the length of the sequence read, because when alignment gaps larger than the sequence read length are present, it becomes difficult for variant callers to reassemble sequence reads and detect structural variants.
[0152] In some embodiments, operation 448 includes aligning the structural variants that were not filtered out in operation 446 to each other. The structural variants may be aligned using any suitable alignment technique, such as, for example, Minimap2. If there are identical subsequences of at least the read length in common, then the smaller of the structural variants are excluded from the filtered graph 454 and included in the filtered structural variants 452.
[0153] In some embodiments, structural variants that are not filtered out in operation 448 are included in the graph reference construct 454. However, because the decoys for the references are captured by additional common sequences that are not in the linear reference, some of those sequences may already be represented by the structural variants included in the graph reference construct 454. Thus, in some embodiments, the structural variants are aligned to the decoy sequences in operation 456. If an alignment is found, those regions in the decoy are masked with a corresponding number of bases. In some embodiments, the masked decoy sequences are concatenated with the linear reference sequence 412 in operation 458 to generate a masked reference 460, which may be output as a FASTA file, or data in any other suitable format.
[0154] 4D illustrates an example process 424 for identifying a set of multiple alignable variants using a second filtering stage, according to some embodiments of the techniques described herein. In some embodiments, process 424 includes processing the variants represented in an initial graph reference construct 462. In some embodiments, initial graph reference construct 462 is the same as graph reference construct 442 (e.g., output in operation 418, where merging is performed). In some embodiments, initial graph reference construct 464 is the same as filtered graph reference construct 454 output from process 422.
[0155] In some embodiments, operation 468 includes simulating reads that traverse all possible paths in the graph reference construct 464. This may include, for example, traversing the graph reference construct 464 at a specified interval for the start position. For a given start position, all possible paths of a specified length are generated as reads for that position. In some embodiments, the generated reads are collected as a FASTQ file, or data in any other suitable format.
[0156] In some embodiments, operation 470 includes aligning the reads to the graph reference 464 using any suitable alignment technique.
[0157] In some embodiments, operation 472 includes filtering the variants based on the alignment. In some embodiments, filtering the variants includes grouping reads at the same starting position. In the group, there will be one read that corresponds only to the linear reference 462, and the rest will follow possible combinations of edges in the graph construct 464.
[0158] After the reads are grouped, in some embodiments, non-reference reads will be identified for exclusion from the filtered graph reference construct 480 (e.g., classified as "bad"). A read is classified as bad if it has a mapping quality of 0 where the reference read had a mapping quality greater than 0. In some embodiments, non-reference reads will be identified for inclusion in the filtered graph reference construct 480 (e.g., classified as "good") if the read has a mapping quality greater than the mapping quality of the reference read or greater than a threshold value (e.g., 20). In some embodiments, non-reference reads will not be ignored if they do not meet the criteria specified above.
[0159] In some embodiments, if there are reads that have different starting locations but follow the same combination of edges, the reads are aggregated. If the aggregated group of reads contains only reads classified as "bad", then the edge combination is identified (e.g., flagged) for filtration.
[0160] In some embodiments, a minimal subset of edges is identified from the combinations of edges identified for filtration in operation 476. For example, a minimal subset of edges may be identified such that each combination of flagged edges has at least one edge in common with the subset.
[0161] In some embodiments, in operation 478 , variants associated with a subset of edges are included in the filtered set of variants 430 and excluded from the filtered graph construction 480 .
[0162] Illustrative Examples Experiments were conducted to evaluate the performance of graph reference construction obtained using the techniques described herein. The experiments show that graph construction can significantly improve both read alignment and variant calling accuracy while being computationally efficient. The results are compared to those obtained using conventional linear non-graph-based techniques, and clearly show that the graph-based approach achieves a significant reduction in read mapping errors, an increase in variant calling sensitivity, and leads to improved joint variant calling without the use of computationally intensive post-processing steps. The conventional technique uses BWA-MEM to align sequence reads to a linear reference, and then uses GATK to identify data differences against the linear reference (variant calling). The conventional technique is referred to herein as "BWA+GATK". BWA-MEM is described by Li H. and Durbin R. ("Fast and accurate short read alignment with Burrows-Wheeler Transform". Bioinformatics, 25:1754-60, 2009). GATK is described by McKenna A et al. ("The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res", 20:1297-303, 2010), which is incorporated by reference in its entirety.
[0163] To demonstrate the capabilities of the techniques described herein, one pan-genome graph and six population-specific graphs were generated according to the techniques described herein, including at least the description of Figures 1-4D. Public databases, such as gnomAD and UK BioBank, were used to construct the pan-genome graph. Similarly, public databases were used to construct the initial population-specific graphs. The initial population-specific graphs were then iteratively extended using construction sets that included Illumina sequencing data of African samples from the 1000 Genomes Project. Pan-African 0 refers to the population-specific graphs obtained using only gnomAD, and Pan-African 5 is the final graph obtained after all five construction sets were added to the graph. In constructing the Pan-African 1 graph, high-quality SVs curated by the Human Genome Structural Variants Consortium (HGSVC) using PacBio HiFi sequencing data for 10 African samples in the 1000 Genomes dataset are also incorporated. Table 1 lists the graph and index memory usage, as well as the total number of variants in each graph. Table 2 shows the contents of each graph.
[0164] [Table 1]
[0165] [Table 2]
[0166] [Table 3]
[0167] We compared pan-genome and population-specific graph references with BWA-MEM for alignment and GATK for variant calling. First, we compared the alignment accuracy as shown in Figure 5. Each panel shows different alignment statistics as violin diagrams. Each violin corresponds to a different graph reference and represents the median and distribution of the statistics across all benchmarking samples. Panel (a) shows the percentage of unmapped reads. BWA maps more reads than any of the graph references. This is due to the lenient alignment approach used by BWA in contrast to the more stringent criteria used by the graph aligners. The percentage of inappropriate reads (classified as either inappropriate orientation for the read pair or insert length outside the expected range) and uninformative reads (MAPQ<20) appears much lower for the graph approach compared to BWA. The ratio of multiple mapping reads is also higher for BWA compared to any of the graph approaches. As can be easily seen from this example, using the graph reference construct generated using the techniques developed by the inventors and described herein results in improved read alignment over conventional techniques. By excluding variants that may introduce ambiguity into the graph reference construct (e.g., by excluding one or more structural variants and / or multiply alignable variants), fewer reads align to multiple locations within the graph reference construct compared to conventional techniques, resulting in more accurate and reliable alignment results.
[0168] A useful metric for measuring the representativeness of population-specific graphs is the alignment error rate, i.e., the base-by-base mismatch rate against the genomic reference. A smaller error rate indicates that the genetic composition of the population is better captured and also the reference bias is reduced. Panel (f) of Figure 5 shows that the error rate decreases consistently from the linear approach to the Pan-African graph. Each extension of the Pan-African graph achieves better error rates, resulting in a reduction of around 50% compared to BWA in the final iteration.
[0169] We also measured the utility of population-specific graphs for variant calling. A graph-aware variant caller that can utilize the information stored in the graph reference was used for variant calling. In Figure 6, the overall performance for single nucleotide polymorphisms (SNPs), insertions and deletions (INDELs), and structural variants (SVs) for all graph references is shown. Panels (a) and (c) show the number of SNPs and INDELs discovered per sample, respectively. It appears that the pan-genome graph provides higher sensitivity compared to the BWA+GATK pipeline. Thus, using the graph reference constructs developed by the inventors and generated using the techniques described herein allows for improved variant calling over conventional techniques.
[0170] Panel (e) of FIG. 6 shows the number of SVs detected by each pipeline (SVs are defined as variants longer than 50 base pairs). Panel (f) also shows the size distribution of SVs for the BWA+GATK, pangenome, pan-African0, and pan-African5 pipelines. It can be seen that the linear approach using BWA+GATK has a significantly lower SV detection rate and is only able to detect short SVs. The pangenome graph provides a significant improvement over the linear approach. This is made possible by the addition of alt contigs in the GRCh38 assembly as alternative paths into the graph reference. Thus, using the graph reference construct generated using the techniques developed by the inventors and described herein allows for more accurate variant calling. By merging variants from different sources and including the merged variants in the graph reference construct, the resulting graph reference construct can be used for more accurate variant calling.
[0171] Using the output of the last iteration as a final graph reference, we compare in more detail the variant calls made by the Pan-African 5 pipeline and those made by the BWA+GATK pipeline. Figure 7 shows the cumulative variant counts for both pipelines against allele frequency. Variants are first classified into SNPs and INDELs (panels A and B, respectively), and then into common (populations detected by both pipelines) and unique (those detected by either pipeline) variant sets. A high concordance is observed between the pipelines, since the majority of variants are detected by both pipelines (solid lines). To distinguish the genotyping effectiveness of these methods, common variants are classified into SNPs and INDELs (panels A and B, respectively), and then into common (populations detected by both pipelines) and unique (those detected by either pipeline) variant sets. GRAF >AF GATK and A.F. GATK >AF GRAFWe further divide the graph into two categories as (dotted lines). The former represents the number of variants detected in the population by both methods but genotyped with higher sensitivity by the graph pipeline (and vice versa for the latter). Among variants observed in the population with high frequency (≧5%), the graph pipeline is able to genotype approximately 120k INDELs and 119k SNPs with higher AF, compared to the same numbers of 106k INDELs and 51k SNPs for GATK. In addition, it is noteworthy that the graph-based approach identifies approximately six times as many unique variants as the linear method.
[0172] To predict the potential clinical significance of variants detected by the graph-based approach and to exclude any bias in variant calling sensitivity towards specific genomic regions or spread in the population, all detected variants were stratified into exonic, intronic, and intergenic regions, as shown in Figure 8. Variants were further divided into three frequency bins as singleton (observed in only one sample), rare (AF<5%, but observed in multiple samples), and common (AF>=5%), and the results were compared with the linear approach BWA+GATK. It can be seen that the use of the Pan-African graph results in the detection of 3-4 times more high- and medium-impact variants in exonic regions for all frequency bins compared to the BWA+GATK pipeline (Panel F). Specifically, there are 429 and 9457 more high- and medium-impact variants detected by the graph pipeline, respectively. It can be seen that the use of the Pan-African graph results in the detection of 3-4 times more high and medium impact variants in exonic regions for all frequency bins compared to the BWA+GATK pipeline (Panel F). Specifically, there are 429 and 9457 more high and medium impact variants detected by the graph pipeline, respectively. As is evident from this example, the use of the graph reference constructs developed by the inventors and generated using the techniques described herein allows for increased sensitivity in variant calling over conventional techniques. The increased sensitivity allows for the detection of more high and medium impact variants that can be used to predict the clinical significance of the detected variants.
[0173] Further implementation details FIG. 9 illustrates an exemplary implementation of a computer system 900 that may be used in connection with any of the embodiments of the technology described herein (e.g., the processes described with reference to FIGS. 2A-D and 4A-4D). The computer system 900 includes one or more computer hardware processors 910 and one or more articles of manufacture that include non-transitory computer-readable storage media (e.g., memory 920 and one or more non-volatile storage media 930). The processor 910 may control the writing of data to and reading of data from the memory 920 and the non-volatile storage device 930 in any suitable manner, as aspects of the technology described herein are not limited in this respect. To perform any of the functionality described herein, the processor 910 may execute one or more processor-executable instructions stored in one or more non-transitory computer-readable storage media (e.g., memory 920), which may act as a non-transitory computer-readable storage medium that stores the processor-executable instructions for execution by the processor 910.
[0174] Computing device 900 may also include a network input / output (I / O) interface 940 through which the computing device may communicate with other computing devices (e.g., over a network), and may also include one or more user I / O interfaces 950 through which the computing device may provide output to a user and receive input from a user. The user I / O interfaces may include devices such as a keyboard, a mouse, a microphone, a display device (e.g., a monitor or touch screen), speakers, a camera, and / or various other types of I / O devices.
[0175] The above-described embodiments may be implemented in any of numerous ways. For example, the embodiments may be implemented using hardware, software, or a combination thereof. When implemented in software, the software code may be executed on any suitable computer hardware processor (e.g., one or more microprocessors, one or more graphic processing units (GPUs)) or collection of computer hardware processors, whether provided in a single computing device or distributed among multiple computing devices. Additionally or alternatively, the embodiments may be implemented using one or more application specific integrated circuits (ASICs) and / or one or more field programmable gate arrays (FPGAs). Thus, the embodiments may be implemented using any suitable computing device (e.g., one or more computer hardware processors, one or more ASICs, and / or one or more FPGAs).
[0176] In this regard, it should be appreciated that an implementation of the embodiments described herein comprises at least one non-transitory computer-readable storage medium (e.g., RAM, ROM, EEPROM, flash memory, or other memory technology, CD-ROM, digital versatile disk (DVD), or other optical disk storage, magnetic cassette, magnetic tape, magnetic disk storage, or other magnetic storage device, or other tangible non-transitory computer-readable storage medium) encoded with a computer program (e.g., a plurality of executable instructions) that, when executed on one or more computer hardware processors, performs the above-described functions of one or more embodiments. The computer-readable medium may be portable such that a program stored thereon may be loaded onto any computing device to implement aspects of the techniques described herein. In addition, it should be appreciated that reference to a computer program that, when executed, performs any of the above-described functions is not limited to an application program running on a host computer. Rather, the terms computer program and software are used in a general sense herein to refer to any type of computer code (e.g., application software, firmware, microcode, or any other form of computer instructions) that can be utilized to program one or more processors to implement aspects of the techniques described herein.
[0177] The above description of implementations provides illustration and description, but is not intended to be exhaustive or to limit the implementations to the precise form disclosed. Modifications and variations are possible in light of the above teachings or may be acquired from practice of the implementations. In other implementations, the methods depicted in these figures may include fewer operations, different operations, operations in a different order, and / or additional operations. Additionally, non-dependent blocks may be performed in parallel.
[0178] The terms "program" or "software" are used in a general sense herein to refer to any type of computer code or set of computer-executable instructions that can be utilized to program a computer or other processor to perform various aspects as described above. In addition, according to one aspect, it should be understood that one or more computer programs that, when executed, perform the methods of the present disclosure need not reside on a single computer or processor, but may be distributed in a modular manner among a number of different computers or processors to perform various aspects of the present disclosure.
[0179] Computer-executable instructions may be in many forms, such as program modules, executed by one or more computers or other devices. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform particular tasks or implement particular abstract data types. Typically the functionality of the program modules may be combined or distributed as desired in various embodiments.
[0180] Also, the data structures may be stored in the computer-readable medium in any suitable format. For ease of illustration, the data structures may be depicted as having fields that are related through locations within the data structure. Such relationships may similarly be achieved by allocating storage for the fields with locations within the computer-readable medium that convey the relationships between the fields. However, any suitable mechanism may be used to establish relationships between information within fields of the data structures, such as through the use of pointers, tags, or other mechanisms that establish relationships between data elements.
[0181] When implemented in software, the software code may be executed on any suitable processor or collection of processors, whether provided within a single computer or distributed among multiple computers.
[0182] Further, it should be understood that a computer may be embodied in any of a number of forms, such as, by way of non-limiting example, a rack-mounted computer, a desktop computer, a laptop computer, or a tablet computer. In addition, a computer may be incorporated within devices not generally considered computers, but having suitable processing capabilities, including a Personal Digital Assistant (PDA), a smartphone, a tablet, or any other suitable portable or fixed electronic device.
[0183] All definitions, as defined and used herein, should be understood to take precedence over dictionary definitions, definitions in documents incorporated by reference, and / or ordinary meanings of the defined terms.
[0184] The indefinite articles "a" and "an," as used in the specification and claims, unless clearly indicated to the contrary, should be understood to mean "at least one."
[0185] The phrase "and / or," as used in the specification and claims, should be understood to mean "either or both" of the elements so conjoined, i.e., elements that are conjunctive in some cases and disjunctive in other cases. Multiple elements listed with "and / or" should be construed in the same manner, i.e., "one or more" of the elements so conjoined. Other elements other than the elements specifically identified by the "and / or" clause may optionally be present, whether related or unrelated to those elements specifically identified. Thus, as a non-limiting example, a reference to "A and / or B," when used in conjunction with open-ended language such as "comprising," may refer in one embodiment to only A (optionally including elements other than B), in another embodiment to only B (optionally including elements other than A), in yet another embodiment to both A and B (optionally including other elements), etc.
[0186] As used in this specification and claims, the phrase "at least one" referring to a list of one or more elements should be understood to mean at least one element selected from any one or more of the elements in the list of elements, but not necessarily including at least one of each and every element specifically listed in the list of elements, and not excluding any combination of elements in the list of elements. This definition also allows that elements other than those specifically identified in the list of elements to which the phrase "at least one" refers may optionally be present, whether related or unrelated to those elements specifically identified. Thus, as a non-limiting example, "at least one of A and B" (or, equivalently, "at least one of A or B" or, equivalently, "at least one of A and / or B") can refer to, in one embodiment, at least one A, optionally including two or more, where B is not present (and optionally including elements other than B), in another embodiment, at least one B, optionally including two or more, where A is not present (and optionally including elements other than A), in yet another embodiment, at least one A, optionally including two or more, and at least one B, optionally including two or more (and optionally including other elements), etc.
[0187] In the claims and in the above specification, all transitional phrases such as "comprising," "including," "carrying," "having," "containing," "involving," "holding," "composed of," and the like, are to be understood to be open-ended, i.e., to mean including but not limited to. Only the transitional phrases "consisting of" and "consisting essentially of," respectively, are intended to be closed or semi-closed transitional phrases.
[0188] The terms "approximately," "substantially," and "about" can be used to mean, in some embodiments, within ±20% of a target value, in some embodiments, within ±10% of a target value, in some embodiments, within ±5% of a target value, and in some embodiments, within ±2% of a target value. The terms "approximately," "substantially," and "about" can include the target value.
Claims
1. 1. A method for aligning sequence reads to a graph reference sequence construct, the method comprising: Obtaining a plurality of variants associated with a linear reference sequence construct for at least one portion of a genome; generating the graph reference sequence construct using the plurality of variants and the linear reference sequence construct, said generating comprising: filtering the plurality of variants to obtain a filtered set of variants, the filtered set of variants being a subset of the plurality of variants, the filtering comprising a first filtering step and a second filtering step different from the first filtering step and performed after the first filtering step; the first filtering step includes identifying a first subset of variants from the plurality of variants, at least in part, by excluding one or more structural variants from the plurality of variants, the one or more structural variants including a first structural variant; the second filtering step includes identifying the filtered set of variants from among the first subset of variants, at least in part, by excluding one or more multiply alignable variants from the first subset of variants. Filtering, generating the graph reference sequence construct using the filtered set of variants and the linear reference sequence construct; generating aligning the sequence reads to the generated graph reference sequence construct; generating a graph reference sequence construct using a method comprising using at least one computing device executing the method.
2. Identifying the first subset of variants from among the plurality of variants comprises: determining whether a first length of the first structural variant exceeds a first specified threshold; and removing the first structural variant from the plurality of variants upon determining that the first length exceeds the first specified threshold; Including, Optionally, the first structural variant is an insertion event and determining whether the first length of the first structural variant exceeds the first specified threshold comprises determining whether the first length is at least 5,000 base pairs; or the first structural variant is a deletion event, and determining whether the first length of the first structural variant exceeds the first specified threshold comprises determining whether the first length is at least 90,000 base pairs. The method of claim 1.
3. Identifying the first subset of variants from among the plurality of variants comprises: (i) aligning the first structural variant to the linear reference sequence construct, determining whether the linear reference sequence construct contains a subsequence of length exceeding a specified threshold, wherein the subsequence is identical to at least a portion of the first structural variant; and upon determining that the linear reference sequence construct comprises the subsequence, excluding the first structural variant from the plurality of variants; (ii) aligning the first structural variant to one or more variants of the plurality of variants, the one or more variants being distinct from the first structural variant; determining whether a second structural variant comprises a subsequence, said subsequence being identical to at least a portion of said first structural variant; and upon determining that the second structural variant comprises the subsequence, removing one of the first structural variant or the second structural variant from the plurality of variants; and (iii) aligning the first structural variant to a decoy sequence associated with the reference sequence construct, the decoy sequence comprising a consensus sequence not included in the linear reference sequence construct; determining whether a decoy sequence associated with the linear reference sequence construct comprises a subsequence of length exceeding a specified threshold, the subsequence being identical to at least a portion of the first structural variant; and masking the decoy sequence upon determining that the decoy sequence comprises the subsequence; The method of claim 1 or 2, comprising one or more of the following:
4. identifying the first subset of variants from among the plurality of variants includes, upon determining that the first length does not exceed the first specified threshold, determining whether the linear reference sequence construct comprises a first subsequence of length greater than a second specified threshold, optionally wherein the second specified threshold is at least 150 base pairs, and wherein the first subsequence is identical to at least a first portion of the first structural variant; and upon determining that the linear reference sequence construct comprises the first subsequence, excluding the first structural variant from the plurality of variants; The method of claim 2 , further comprising:
5. determining whether a second structural variant comprises a second subsequence of length greater than a second specified threshold, optionally wherein the second specified threshold is at least 150 base pairs, and wherein the second subsequence is identical to at least a second portion of the first structural variant, if the linear reference sequence construct does not comprise the first subsequence; upon determining that the second structural variant comprises the second subsequence, removing one of the first structural variant or the second structural variant from the plurality of variants; Optionally, excluding one of the first structural variant or the second structural variant from the plurality of variants comprises: identifying a shortest variant from among the first structural variant and the second structural variant; and removing the shortest variant from the plurality of variants; The method of claim 4 , comprising:
6. If it is determined that the second structural variant does not include the second partial sequence, determining whether a decoy sequence associated with the linear reference sequence construct includes a third partial sequence having a length greater than a specified threshold, the third partial sequence being identical to at least a third portion of the first structural variant, and the decoy sequence includes a common sequence not included in the linear reference sequence construct; and masking the decoy sequence upon determining that the decoy sequence comprises the third subsequence; The method of claim 5 further comprising:
7. Identifying the filtered set of variants from among the first subset of variants includes: generating an initial graph reference sequence build using at least a portion of said first subset of variants; and generating a plurality of graph reads using the initial graph reference construct, wherein at least a portion of the plurality of graph reads are each associated with a respective path within the initial graph reference sequence construct; The method of any one of claims 1 to 6, further comprising:
8. The plurality of graph reads traverse an initial graph reference construct using a moving window with jumps, and / or the plurality of graph reads includes a first subset of graph reads and a second subset of graph reads, and generating the plurality of graph reads includes: generating the first subset of graph reads by traversing the initial graph reference sequence construction over a first interval; and generating the second subset of graph reads by traversing the initial graph reference sequence construction over a second interval, wherein the first interval and the second interval at least partially overlap; The method of claim 7, comprising:
9. The method further comprises aligning at least a portion of the plurality of graph reads to the initial graph reference sequence construct, wherein the aligning comprises, for each of the at least a portion of the plurality of graph reads, determining an alignment quality between the graph reads and the graph reference sequence construct; and determining whether the alignment quality exceeds a threshold; The method of claim 7 or 8, comprising:
10. identifying a first group of the at least some of the plurality of graph reads, each graph read included in the first group of the at least some of the plurality of graph reads comprising a first combination of one or more variants of the first subset of variants; Optionally, the first group of the at least some of the plurality of graph leads comprises a first graph lead and a second graph lead; and upon determining that neither the first alignment quality determined for the first graph read nor the second alignment quality determined for the second graph read exceeds the specified threshold, excluding at least one multiply alignable variant from the filtered set of variants, wherein the at least one multiply alignable variant is included within the first combination of the one or more variants.
10. The method of claim 9.
11. Identifying the filtered set of variants from among the first subset of variants includes: generating an initial graph reference sequence construct using said first subset of variants; traversing said initial graph reference sequence construction to generate a plurality of graph reads; aligning the plurality of graph reads to the initial graph reference sequence construct and determining an alignment quality for each of at least a portion of the plurality of graph reads, the alignment quality comprising determining whether the graph read maps to two or more regions within the initial graph reference sequence construct; and excluding at least a portion of the one or more of the first set of variants from the second set of variants based on the alignment quality, optionally wherein one or more of the plurality of graph reads are associated with the same combination of one or more of the first subset of variants; Including, determining whether each of the alignment qualities determined for the one or more of the plurality of graph reads is below a specified threshold; removing at least one variant from the filtered set of variants upon determining that each of the alignment qualities is below the specified threshold; The method of any one of claims 1 to 10, further comprising:
12. Obtaining the plurality of variants includes: Obtaining a plurality of alternative sequences associated with said reference sequence construct; and processing at least a portion of the plurality of alternative sequences, the processing comprising, for a first alternative sequence of the plurality of alternative sequences: aligning said first alternative sequence to said linear reference sequence construct and obtaining an aligned position; identifying one or more differences between the first alternative sequence and the reference sequence construct at the aligned position; and including at least a portion of the one or more differences as a first variant within the plurality of variants; Including, Optionally, the first alternative sequence comprises an inverted sequence patch, and aligning the first alternative sequence to the linear reference sequence construct and obtaining the aligned position comprises obtaining an alternative aligned position for the inverted sequence patch; and / or the method further comprises left-normalizing the first variant to the reference sequence construct prior to including the first variant in the plurality of variants. The method according to any one of claims 1 to 11.
13. the at least some of the one or more differences include consecutive first and second differences, the first difference being associated with a first subsequence of the first alternative sequence and the second difference being associated with a second subsequence of the linear reference sequence construct; and processing the first and second differences prior to including them in the plurality of variants as a first variant, the processing further comprising: determining whether the first subsequence contains one or more regions contained within the second subsequence; and upon determining that the first subsequence contains the one or more regions contained within the second subsequence, removing the one or more regions from both the first and second subsequences; and optionally, the first and second differences comprise an insertion and a deletion event, respectively. The method of claim 12.
14. Obtaining the plurality of variants includes: Obtaining a second variant associated with the linear reference sequence construct; 14. The method of claim 12 or 13, further comprising: including the second variant within the plurality of variants; and optionally annotating the second variant with information indicating a source of the second variant; and / or for a shared variant contained within both the at least some of the first variants and the at least some of the second variants, averaging a first allele frequency associated with at least some of the first variants that include the shared variant and a second allele frequency associated with at least some of the second variants that include the shared variant to obtain an average allele frequency.
15. 1. A system comprising: at least one computer hardware processor; at least one non-transitory computer-readable storage medium storing processor-executable instructions; wherein the processor-executable instructions, when executed by the at least one computer hardware processor, cause the at least one computer hardware processor to perform the method recited in any one of claims 1 to 14.
16. At least one non-transitory computer readable storage medium storing processor executable instructions that, when executed by at least one computer hardware processor, cause the at least one computer hardware processor to perform the method recited in any one of claims 1 to 14.