Methods and systems for identifying recombinant mutants

Abstract

Disclosed herein are systems, devices, and methods for identifying recombinant variants (e.g., gene conversion variants) of genes, such as the GBA gene and the CYP21A2 gene, the copy number of the recombinant variants, and the genetic variant status (e.g., carrier, compound heterozygous, or homozygous).

Description

[Technical field]

[0001] Related Applications This application claims priority under 35 U.S.C. §119(e) to U.S. Provisional Application No. 63 / 197,936, filed June 7, 2021. The contents of the related application are incorporated herein by reference in their entireties.

[0002] The present disclosure relates broadly to the field of determining genetic variants, and more particularly to determining recombinant variants. [Background technology]

[0003] Segmental duplications are hotspots for structural variants and genetic variants (e.g., gene conversion). Segmental duplications can occur for genes with highly homologous gene family members or pseudogenes. High sequence similarity of genes and homologous gene family members or pseudogenes can result in poor read alignment and variant calling. There is a need to informatically identify variants of genes with highly homologous gene family members or pseudogenes. Summary of the Invention [Means for solving the problem]

[0004] Disclosed herein is a method for determining GBA status. In some embodiments, the method for determining GBA status is under the control of a processor (such as a hardware processor or a virtual processor) and includes receiving a first plurality of sequence reads generated from a sample obtained from a subject. The method may include aligning the first plurality of sequence reads to a reference genome sequence to obtain a second plurality of sequence reads aligned to the GBA gene or the GBAP1 gene in the reference genome sequence. The method may include determining a number of sequence reads of the second plurality of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference genome sequence. The method may include determining a normalized number of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference genome sequence. The method may include determining the total copy number of the GBA gene and the GBAP1 gene using a mixture of Gaussians including a plurality of Gaussians, each of which represents a different integer copy number, given the normalized number of sequence reads aligned to the region between the GBA gene and the GBAP1 gene. The method may include phasing one or more haplotypes derived from the GBA gene or GBAP1 gene in a region of the GBA gene or a corresponding region of the GBAP1 gene that includes a plurality of GBA / GBAP1 discriminatory bases using sequence reads of a second plurality of sequence reads aligned to a region including a plurality of GBA / GBAP1 discriminatory bases. The method may include determining a copy number of each of the one or more haplotypes using a total copy number of the GBA gene and the GBAP1 gene and a number of sequence reads of the second plurality of sequence reads that each include one or more of the plurality of GBA / GBAP1 discriminatory bases that support the haplotype. The method may include determining a GBA status of the subject using one or more haplotypes derived from the GBA gene or GBAP1 gene in a region of the GBA gene or a corresponding region of the GBAP1 gene and / or a copy number of each of the one or more haplotypes.In some embodiments, the method includes generating a user interface (UI) that includes UI elements that represent or include the GBA state.

[0005] In some embodiments, the unique region between the GBA gene and the GBAP1 gene in the reference genome sequence comprises a unique region of about 10 kilobases in length. The unique region between the GBA gene and the GBAP1 gene in the reference genome sequence may comprise chr1:155220429-155230539 of hg38 or the corresponding region of the reference human genome sequence.

[0006] In some embodiments, determining the normalized number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference genome sequence comprises determining the normalized number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference genome sequence using (1a) a depth of the sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene, (1b) a length of the unique region, (2a) a depth of the sequence reads of a first plurality of sequence reads aligned to each of a plurality of regions in the reference genome sequence other than the locus comprising the GBA gene and the GBAP1 gene, and (2b) a length of each of the plurality of regions of the reference genome other than the locus comprising the GBA gene and the GBAP1 gene.

[0007] In some embodiments, the method comprises determining a normalized, corrected number of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference genome sequence from the normalized number of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference genome sequence. Determining the normalized, corrected number of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference genome sequence may comprise determining a normalized, GC content-corrected number of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference genome sequence from the normalized number of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference genome sequence. Determining the normalized GC content corrected number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference genome sequence may include determining the normalized GC content corrected number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference genome sequence from the normalized number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference genome sequence using (1) the GC content of the unique region between the GBA gene and the GBAP1 gene, and optionally (2) the GC content of each of one or more regions in the reference genome sequence other than the locus including the GBA gene and the GBAP1 gene. Determining the total copy number of the GBA gene and the GBAP1 gene may include determining the total copy number of the GBA gene and the GBAP1 gene using a Gaussian mixture model given the normalized corrected number of sequence reads aligned to the region between the GBA gene and the GBAP1 gene. In some embodiments, the method includes generating a user interface (UI) including a UI element representing or including a CYP21A2 status.

[0008] In some embodiments, determining the total copy number of the GBA and GBAP1 genes comprises determining the copy number of the region between the GBA and GBAP1 genes using a Gaussian mixture model given a normalized number of sequence reads aligned to the region between the GBA and GBAP1 genes. The total copy number of the GBA and GBAP1 genes is the copy number of the region between the GBA and GBAP1 genes plus 2.

[0009] In some embodiments, determining the total copy number of the GBA gene and the GBAP1 gene comprises determining the total copy number of the GBA gene and the GBAP1 gene using a Gaussian mixture model and a predetermined posterior probability threshold given the normalized number of sequence reads aligned to the region between the GBA gene and the GBAP1 gene. The predetermined posterior probability threshold can be 0.95.

[0010] In some embodiments, the Gaussian mixture model includes a one-dimensional Gaussian mixture model. The plurality of Gaussians of the Gaussian mixture model can represent integer copy numbers from 0 to 10. The plurality of Gaussians of the Gaussian mixture model can include 5 Gaussians. The average of each of the plurality of Gaussians can be the integer copy number represented by the Gaussian.

[0011] In some embodiments, phasing one or more haplotypes derived from the GBA or GBAP1 gene comprises analyzing linkage information between the GBA / GBAP1 discriminating bases of the plurality of GBA / GBAP1 discriminating bases using sequence reads of the second plurality of sequence reads aligned to a region including or corresponding to the plurality of GBA / GBAP1 discriminating bases. Phasing one or more haplotypes derived from the GBA or GBAP1 gene may comprise phasing one or more haplotypes derived from the GBA or GBAP1 gene using sequence reads of the second plurality of sequence reads aligned to two or more of the plurality of GBA / GBAP1 discriminating bases, respectively.

[0012] In some embodiments, the sequence reads of the second plurality of sequence reads are aligned to a region of the GBA gene that contains a plurality of GBA / GBAP1 discriminatory bases or a corresponding region of the GBAP1 gene with an alignment quality score of 0 or greater.

[0013] In some embodiments, the region of the GBA gene or the corresponding region of the GBAP1 gene that includes a plurality of GBA / GBAP1 discriminator bases is about 1.1 kilobases in length. The region of the GBA gene or the corresponding region of the GBAP1 gene that includes a plurality of GBA / GBAP1 discriminator bases may include exons 9-11 of the GBA gene or the GBAP1 gene, respectively. The region of the GBA gene or the corresponding region of the GBAP1 gene that includes a plurality of GBA / GBAP1 discriminator bases may include p.L483P, p.D448H, c.1263del, RecNciI, RecTL, and c.1263del+RecTL. The plurality of GBA / GBAP1 discriminator bases may include 10 GBA / GBAP1 discriminator bases.

[0014] In some embodiments, the one or more haplotypes comprise a wild-type GBA haplotype, a wild-type GBAP1 haplotype, and / or a GBA / GBAP1 hybrid haplotype. The GBA / GBAP1 hybrid haplotype can comprise a GBA mutant haplotype or a GBAP1 mutant haplotype.

[0015] In some embodiments, determining the copy number of each of the one or more haplotypes comprises determining that the likelihood of one copy of the wildtype GBA haplotype is higher than the likelihood of two copies of the wildtype GBA haplotype given a number of sequence reads of the second plurality of sequence reads each including one or more of the plurality of GBA / GBAP1 discriminating bases supporting a wildtype GBA haplotype. Determining the copy number of each of the one or more haplotypes may include determining that the copy number of the wildtype GBA haplotype is 1. Determining that the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype includes determining whether or not, in each of one or more pairs (or all pairs) of consecutive GBA / GBAP1 discriminatory bases of a plurality of GBA / GBAP1 discriminatory bases, where a first haplotype of the one or more haplotypes includes a GBA base in the consecutive GBA / GBAP1 discriminatory bases and a second haplotype of the one or more haplotypes includes a GBA base and a GBAP1 base, or a GBAP1 base and a GBA base, in the consecutive GBA / GBAP1 discriminatory bases, (1) a second haplotype that includes a GBA base in the consecutive GBA / GBAP1 discriminatory bases, and a second haplotype that includes a GBA base in the consecutive GBA / GBAP1 discriminatory bases, and The method can include determining that the likelihood of one copy of a wild-type GBA haplotype is higher than the likelihood of two copies of a wild-type GBA haplotype given the number of sequence reads of the two plurality of sequence reads, (2) the number of sequence reads of the second plurality of sequence reads that each include a GBA base and a GBAP1 base in consecutive GBA / GBAP1 discriminating bases, (3) the number of sequence reads of the second plurality of sequence reads that each include a GBAP1 base and a GBA base, or a GBA base and a GBAP1 base, in consecutive GBA / GBAP1 discriminating bases, and / or (4) the number of sequence reads of the second plurality of sequence reads that each include a GBAP1 base in consecutive GBA / GBAP1 discriminating bases. The likelihood of one copy of a wild-type GBA haplotype can include a sum (e.g., a weighted or unweighted average) of the likelihoods of one copy of a wild-type GBA haplotype determined for each of one or more pairs of consecutive GBA / GBAP1 discriminating bases.The likelihood of two copies of the wild-type GBA haplotype may include the sum (e.g., weighted or unweighted average) of the likelihoods of two copies of the wild-type GBA haplotype determined for each of one or more pairs of consecutive GBA / GBAP1 discriminatory bases.

[0016] In some embodiments, the copy number of the wild type GBA haplotype is 1. Determining the GBA status of the subject may include determining that the subject is a carrier of a GBA variant haplotype. In some embodiments, the one or more haplotypes include four haplotypes. The total copy number of the GBA and GBAP1 genes can be four. The copy number of each of the four haplotypes can be 1. Determining the GBA status of the subject may include determining that the subject is a carrier of a GBA variant haplotype.

[0017] In some embodiments, the one or more haplotypes include two or more GBA mutant haplotypes. None of the two or more GBA mutant haplotypes can include a GBA base at each of the plurality of GBA / GBAP1 discriminating bases. None of the two or more GBA mutant haplotypes can include a GBA base at all of the plurality of GBA / GBAP1 discriminating bases. Determining the GBA status of the subject can include determining that the subject is a compound heterozygote for a GBA mutant haplotype.

[0018] In some embodiments, the method includes determining, using sequence reads of a second plurality of sequence reads each including a base in the GBA / GBAP1 discriminating base that is not a GBA base, a copy number of the GBA base in each of one or more of the plurality of GBA / GBAP1 discriminating bases is 0. The base in the GBA / GBAP1 discriminating base that is not a GBA base is a GBAP1 base. Determining the GBA status may include determining that the subject is homozygous for each of one or more of the plurality of GBA / GBAP1 discriminating bases.

[0019] In some embodiments, the first plurality of sequence reads comprises sequence reads that are each about 100 base pairs to about 1000 base pairs in length. In some embodiments, the first plurality of sequence reads comprises paired-end sequence reads and / or single-end sequence reads. In some embodiments, the first plurality of sequence reads are generated by whole genome sequencing (WGS). The WGS can be clinical WGS (cWGS). In some embodiments, the sample comprises cells, cell-free DNA, cell-free fetal DNA, amniotic fluid, a blood sample, a biopsy sample, or a combination thereof.

[0020] Disclosed herein includes a method for determining CYP21A2 status. In some embodiments, the method for determining CYP21A2 status is under the control of a processor (such as a hardware processor or a virtual processor) and includes receiving a first plurality of sequence reads generated from a sample obtained from a subject. The method may include aligning the first plurality of sequence reads to a reference genome sequence to obtain a second plurality of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference genome sequence. The method may include determining the number of sequence reads of the second plurality of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference genome sequence. The method may include determining a normalized number of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference genome sequence. The method may include determining the total copy number of the CYP21A2 gene and the CYP21A1P pseudogene using a Gaussian mixture model including a plurality of Gaussians, each Gaussian representing a different integer copy number, given a normalized number of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene. The method may include phasing one or more haplotypes derived from the CYP21A2 gene or the CYP21A1P pseudogene in a region of the CYP21A2 gene or a corresponding region of the CYP21A1P pseudogene that includes the plurality of CYP21A2 / CYP21A1P discriminatory bases, using sequence reads of a second plurality of sequence reads aligned to a region including the plurality of CYP21A2 / CYP21A1P discriminatory bases or a corresponding region of the CYP21A2 pseudogene. The method may include determining the copy number of each of the one or more haplotypes using the total copy number of the CYP21A2 gene and the CYP21A1P gene, and the number of sequence reads of a second plurality of sequence reads, each of which includes one or more of the plurality of CYP21A2 / CYP21A1P discriminatory bases supporting the haplotype.The method may include determining the subject's CYP21A2 status using one or more haplotypes derived from the CYP21A2 gene or a CYP21A1P pseudogene in a region of the CYP21A2 gene or a corresponding region of the CYP21A1P pseudogene, and / or the copy number of each of the one or more haplotypes.

[0021] In some embodiments, determining the normalized number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference genome sequence includes determining the normalized number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference genome sequence using (1a) the depth of the sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene, (1b) the length of the unique region, (2a) the depth of the sequence reads of a first plurality of sequence reads aligned to each of a plurality of regions in the reference genome sequence other than the locus comprising the CYP21A2 gene and the CYP21A1P pseudogene, and (2b) the length of each of the plurality of regions of the reference genome other than the locus comprising the CYP21A2 gene and the CYP21A1P pseudogene.

[0022] In some embodiments, the method comprises determining a normalized, corrected number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference genome sequence from a normalized number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference genome sequence. Determining the normalized, corrected number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference genome sequence may comprise determining a normalized, GC content corrected number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference genome sequence from the normalized number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference genome sequence. Determining the normalized GC content-corrected number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference genome sequence may include determining the normalized GC content-corrected number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference genome sequence from the normalized number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference genome sequence using (1) the GC content of the CYP21A2 gene or CYP21A1P pseudogene, and optionally (2) the GC content of each of one or more regions in the reference genome sequence other than the locus including the CYP21A2 gene and the CYP21A1P pseudogene. Determining the total copy number of the CYP21A2 gene and the CYP21A1P pseudogene may include determining the total copy number of the CYP21A2 gene and the CYP21A1P pseudogene using a Gaussian mixture model given the normalized, corrected number of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene.

[0023] In some embodiments, determining the total copy number of the CYP21A2 gene and the CYP21A1P pseudogene comprises determining the total copy number of the CYP21A2 gene and the CYP21A1P pseudogene using a Gaussian mixture model and a predetermined posterior probability threshold given the normalized number of sequence reads aligned to the CYP21A2 gene and the CYP21A1P pseudogene. The predetermined posterior probability threshold can be 0.95.

[0024] In some embodiments, the Gaussian mixture model includes a one-dimensional Gaussian mixture model. The plurality of Gaussians of the Gaussian mixture model can represent integer copy numbers from 0 to 10. The plurality of Gaussians of the Gaussian mixture model can include 5 Gaussians. The average of each of the plurality of Gaussians can be the integer copy number represented by the Gaussian.

[0025] In some embodiments, phasing one or more haplotypes derived from the CYP21A2 gene or CYP21A1P pseudogene comprises analyzing linkage information between the CYP21A2 / CYP21A1P discriminating bases of the plurality of CYP21A2 / CYP21A1P discriminating bases using sequence reads of the second plurality of sequence reads aligned to a region including or corresponding to the plurality of CYP21A2 / CYP21A1P discriminating bases. Phasing one or more haplotypes derived from the CYP21A2 gene or CYP21A1P pseudogene comprises phasing one or more haplotypes derived from the CYP21A2 gene or CYP21A1P pseudogene using sequence reads of the second plurality of sequence reads aligned to two or more of the plurality of CYP21A2 / CYP21A1P discriminating bases, respectively.

[0026] In some embodiments, the sequence reads of the second plurality of sequence reads are aligned to a region of the CYP21A2 gene that contains a plurality of CYP21A2 / CYP21A1P discriminator bases or a corresponding region of the CYP21A1P pseudogene with an alignment quality score of 0 or greater.

[0027] In some embodiments, the plurality of CYP21A2 / CYP21A1P discriminant bases comprises 14 CYP21A2 / CYP21A1P discriminant bases. The 14 CYP21A2 / CYP21A1P discriminant bases may comprise 9 CYP21A2 / CYP21A1P recombinant variants. The 14 CYP21A2 / CYP21A1P discriminant bases may be chr6:32039081 / 32006353, 32039128 / 32006400, 32039132 / 32006404, 32039143 / 32006407, 32039426 / 32006690, 32039548 / 32006812, 32039802 / 3200 32040421 / 32007686, and 32040535 / 32007800, or these bases of a reference human genome sequence.

[0028] In some embodiments, the one or more haplotypes comprise a wild-type CYP21A2 haplotype, a wild-type CYP21A1P, and / or a CYP21A2 / CYP21A1P hybrid haplotype. The CYP21A2 / CYP21A1P hybrid haplotype can comprise a CYP21A2 variant haplotype or a CYP21A1P variant haplotype.

[0029] In some embodiments, determining the copy number of each of the one or more haplotypes comprises determining that the likelihood of one copy of the wild-type CYP21A2 haplotype is higher than the likelihood of two copies of the wild-type CYP21A2 haplotype given the number of sequence reads of the second plurality of sequence reads that each include one or more of the plurality of CYP21A2 / CYP21A1P discriminatory bases supporting a wild-type CYP21A2 haplotype. Determining the copy number of each of the one or more haplotypes may include determining that the copy number of the wild-type CYP21A2 haplotype is 1. In some embodiments, determining that the likelihood of one copy of the wild-type CYP21A2 haplotype is higher than the likelihood of two copies of the wild-type CYP21A2 haplotype comprises: (1) determining that the likelihood of one copy of the wild-type CYP21A2 haplotype is higher than the likelihood of two copies of the wild-type CYP21A2 haplotype in each of one or more pairs of consecutive CYP21A2 / CYP21A1P discriminatory bases of a plurality of CYP21A2 / CYP21A1P discriminatory bases, wherein a first haplotype of the one or more haplotypes comprises a CYP21A2 base in the consecutive CYP21A2 / CYP21A1P discriminatory bases and a second haplotype of the one or more haplotypes comprises a CYP21A2 base and a CYP21A1P base, or a CYP21A1P base and a CYP21A2 base, in the consecutive CYP21A2 / CYP21A1P discriminatory bases; determining that the likelihood of one copy of the wild-type CYP21A2 haplotype is higher than the likelihood of two copies of the wild-type CYP21A2 haplotype given (1) the number of sequence reads of the second plurality of sequence reads that each comprise a CYP21A2 base at the 1P discriminant base, (2) the number of sequence reads of the second plurality of sequence reads that each comprise a CYP21A2 base and a CYP21A1P base at consecutive CYP21A2 / CYP21A1P discriminant bases, (3) the number of sequence reads of the second plurality of sequence reads that each comprise a CYP21A1P base and a CYP21A2 base at consecutive CYP21A2 / CYP21A1P discriminant bases, and / or (4) the number of sequence reads of the second plurality of sequence reads that each comprise a CYP21A1P base at consecutive CYP21A2 / CYP21A1P discriminant bases.The likelihood of one copy of the wild-type CYP21A2 haplotype can include the sum (e.g., weighted or unweighted average) of the likelihoods of one copy of the wild-type CYP21A2 haplotype determined for each of one or more pairs of consecutive CYP21A2 / CYP21A1P discriminatory bases. The likelihood of two copies of the wild-type CYP21A2 haplotype can include the sum (e.g., weighted or unweighted average) of the likelihoods of two copies of the wild-type CYP21A2 haplotype determined for each of one or more pairs of consecutive CYP21A2 / CYP21A1P discriminatory bases.

[0030] In some embodiments, the copy number of the wild-type CYP21A2 haplotype is 1. The method may include determining that the subject is a carrier of a CYP21A2 variant haplotype. In some embodiments, the one or more haplotypes include four haplotypes. The total copy number of the CYP21A2 and CYP21A1P genes can be 4. The copy number of each of the four haplotypes can be 1. Determining the subject's CYP21A2 status may include determining that the subject is a carrier of a CYP21A2 variant haplotype.

[0031] In some embodiments, the one or more haplotypes include two or more haplotypes. None of the two or more haplotypes may include a CYP21A2 base at each of the plurality of CYP21A2 / CYP21A1P discriminatory bases. Each of the two or more haplotypes may not include a CYP21A2 base at all of the plurality of CYP21A2 / CYP21A1P discriminatory bases. The method may include determining that the subject is a compound heterozygote of a CYP21A2 variant haplotype.

[0032] In some embodiments, the one or more haplotypes include only one haplotype. The only one haplotype may not include a CYP21A2 base at each of the multiple CYP21A2 / CYP21A1P discriminatory bases. The method may include determining that the subject is homozygous for a CYP21A2 variant haplotype.

[0033] In some embodiments, the first plurality of sequence reads comprises sequence reads that are each about 100 base pairs to about 1000 base pairs in length. In some embodiments, the first plurality of sequence reads comprises paired-end sequence reads and / or single-end sequence reads. In some embodiments, the first plurality of sequence reads are generated by whole genome sequencing (WGS). The WGS can be clinical WGS (cWGS). In some embodiments, the sample comprises cells, cell-free DNA, cell-free fetal DNA, amniotic fluid, a blood sample, a biopsy sample, or a combination thereof.

[0034] Disclosed herein is a system (e.g., a computing system) for determining genetic variants. In some embodiments, the system for determining genetic variants comprises a non-transitory memory configured to store executable instructions and a first plurality of sequence reads generated from a sample obtained from a subject. The system can comprise a processor, such as a hardware processor or a virtual processor, in communication with the non-transitory memory. The processor can be programmed by the executable instructions to align the first plurality of sequence reads to a reference sequence to obtain a second plurality of sequence reads aligned to genes or gene paralogs, or regions therebetween, in the reference sequence. The processor can be programmed by the executable instructions to determine the total copy number of genes and gene paralogs using a mixture of Gaussians including multiple Gaussians, each Gaussian representing a different integer copy number, given the number of sequence reads aligned to genes or gene paralogs, or regions therebetween. The processor can be programmed with the executable instructions to phase one or more haplotypes derived from a gene (including recombinant variants of a gene) or gene paralog, or a region of a gene or a corresponding region of a gene paralog, comprising a plurality of gene / gene paralog discriminating bases, using sequence reads of a second plurality of sequence reads aligned to a region comprising the plurality of gene / gene paralog discriminating bases or a corresponding region. The processor can be programmed with the executable instructions to determine a copy number of each of the one or more haplotypes using the total copy number of the gene and gene paralog and the number of sequence reads of the second plurality of sequence reads each comprising one or more of the plurality of gene / gene paralog discriminating bases supporting the haplotype.

[0035] In some embodiments, the genetically engineered variants include reciprocal recombinant variants. The genetically engineered variants can include non-reciprocal recombinant variants. In some embodiments, the reference sequence includes a reference genome sequence.

[0036] In some embodiments, the processor is programmed with executable instructions to execute: determining the number of sequence reads aligned to a gene or gene paralog, or a region between them. In some embodiments, the number of sequence reads aligned to a gene or gene paralog, or a region between them, comprises a normalized and / or GC-corrected number of sequence reads aligned to a gene or gene paralog, or a region between them. In some embodiments, the gene paralog is a gene. The gene paralog can be a pseudogene. In some embodiments, the gene and the gene paralog have at least 90% sequence identity.

[0037] In some embodiments, the gene is the GBA gene and the gene paralog is the GBAP1 gene. In some embodiments, the gene is the CYP21A2 gene and the gene paralog is the CYP21A1P pseudogene. In some embodiments, the gene is ABCC6, ABCD1, ACTB, ACTG1, ACTN4, ADAMTSL2, ADIPOR1, AFG3L2, AGK, ALG1, ALMS1, ANKRD11, ANOS1, AP4S1, ARMC4, ARSE, ASNS, ATAD3A, B3GAT3, BCAP31, BDP1, BMPR1A, BRAF, BRCA1, C2, CACNA1C, CALM1, CD46, CEP290, CFH, CFH, CFH, CHEK2, CISD2, CLCNKA , CLCNKB, CORO1A, COX10, CP, CRYBB2, CSF2RA, CUBN, CUBN, CYCS, CYP11B1, CYP21A2, DCLRE1C, DHFR, DICER1, DIS3L2, DNAH11, DNAH11, DN M1, DSE, DUOX2, EGLN1, ELK1, ELMO2, ERCC6, ESPN, EYS, F8, FANCD2, FANCD2, FAR1, FHL1, FLG, FLNC, FOXD4, FXN, GBA, GH1, GJA1, GK, GLUD1, GLUD1, GOSR2, GUSB, HBA1, HBA2, HNRNPA1, HPS1, HSPD1, HYDIN, IDS, IFT122, IGLL1, KANSL1, KCTD1, KIF1C, KRAS, KRT14, KRT16, KRT17, K RT6A, KRT6B, KRT6C, LEFTY2, LRP5, LRP5, MAT2A, MID1, MOCS1, MSN, MSX2, MYO5B, NCF1, NEB, NECAP1, NEFH, NF1, NF1, NF1, NOTCH2, NXF5, O CLN, OTOA, PARN, PBX1, PIGA, PIGN, PIK3CA, PIK3CD, PKD1, PKP2, PMS2, PMS2, PMS2, PNPT1, POLH, PRODH, PRODH, PROS1, PRPS1, PRSS1, PTE N, RAD21, RBM8A, RBPJ, RDX, RMND1, RNF216, RNF216, RPL15, SALL1, SBDS, SDHA, SHOX, SLC25A15, SLC25A15, SLC33A1, SLC6A8, SMN1, SMN2,SOX2, SPTLC1, SRD5A3, SRP72, STAT5B, STRC, SYT14, TARDBP, TBL1XR1, TBX20, TIMM8A, TPM3, TPMT, TRAPPC2, TRIP11, TTN, TU BA1A, TUBB2A, TUBB2B, TUBB3, TUBB4A, TUBG1, TYR, UBA5, UBE3A, UNC93B1, USP18, VPS35, VWF, WRN, XIAP, ZEB2, or ZNF341. ,

[0038] In some embodiments, the processor is programmed with the executable instructions to determine a genetic variant status of the subject using one or more haplotypes derived from a gene or a gene paralog, or a region of the gene or a corresponding region of the gene paralog, and / or the copy number of each of the one or more haplotypes. In some embodiments, the processor can be programmed with the executable instructions to generate a user interface (UI) that includes UI elements that represent or include the genetic variant status.

[0039] In some embodiments, determining the normalized number of sequence reads aligned to genes or gene paralogs in the reference sequence comprises determining the normalized number of sequence reads aligned to genes or gene paralogs in the reference sequence using (1a) the depth of the sequence reads aligned to the genes or gene paralogs, (1b) the length of the unique region, (2a) the depth of the sequence reads of a first plurality of sequence reads aligned to each of a plurality of regions of the reference sequence other than the loci comprising the genes and gene paralogs, and (2b) the length of each of the plurality of regions of the reference sequence other than the loci comprising the genes and gene paralogs.

[0040] In some embodiments, the processor is programmed with the executable instructions to execute: determining a normalized, corrected number of sequence reads aligned to genes or gene paralogs in the reference sequence from the normalized number of sequence reads aligned to genes or gene paralogs in the reference sequence. Determining the normalized, corrected number of sequence reads aligned to genes or gene paralogs in the reference sequence may include determining a normalized, GC content corrected number of sequence reads aligned to genes or gene paralogs in the reference sequence from the normalized number of sequence reads aligned to genes or gene paralogs in the reference sequence. Determining the normalized GC content corrected number of sequence reads aligned to genes or gene paralogs in the reference sequence may include determining the normalized GC content corrected number of sequence reads aligned to genes or gene paralogs in the reference sequence from the normalized number of sequence reads aligned to genes or gene paralogs in the reference sequence using (1) the GC content of the genes or gene paralogs, and optionally (2) the GC content of each of one or more regions of the reference sequence other than the locus containing the genes and gene paralogs. Determining the total copy number of genes and gene paralogs may include determining the total copy number of genes and gene paralogs using a Gaussian mixture model given the normalized corrected number of sequence reads aligned to the genes or gene paralogs.

[0041] In some embodiments, determining the total copy number of genes and gene paralogs comprises determining the copy number of the region between the genes or gene paralogs using a Gaussian mixture model given the normalized number of sequence reads aligned to the genes or gene paralogs. The total copy number of genes and gene paralogs can be the copy number of the region between the genes or gene paralogs plus 2.

[0042] In some embodiments, determining the total copy number of genes and gene paralogs comprises determining the total copy number of genes and gene paralogs using a Gaussian mixture model and a predefined posterior probability threshold given the normalized number of sequence reads aligned to the genes and gene paralogs. The predefined posterior probability threshold can be 0.95.

[0043] In some embodiments, the Gaussian mixture model includes a one-dimensional Gaussian mixture model. The plurality of Gaussians of the Gaussian mixture model can represent integer copy numbers 0 to 10. The plurality of Gaussians of the Gaussian mixture model can include 5 Gaussians. The average of each of the plurality of Gaussians can be the integer copy number represented by the Gaussian.

[0044] In some embodiments, phasing one or more haplotypes derived from a gene or gene paralog comprises analyzing linkage information between gene / gene paralog discriminating bases of the plurality of gene / gene paralog discriminating bases using sequence reads of the second plurality of sequence reads aligned to a region containing the plurality of gene / gene paralog discriminating bases or a corresponding region. Phasing one or more haplotypes derived from a gene or gene paralog may comprise phasing one or more haplotypes derived from a gene or gene paralog using sequence reads of the second plurality of sequence reads aligned to two or more of the plurality of gene / gene paralog discriminating bases, respectively. In some embodiments, sequence reads of the second plurality of sequence reads are aligned to a region of a gene containing the plurality of gene / gene paralog discriminating bases or a corresponding region of a gene paralog having an alignment quality score of 0 or greater.

[0045] In some embodiments, the one or more haplotypes comprise a wildtype gene haplotype, a wildtype gene paralog, and / or a gene / gene paralog hybrid haplotype. A gene / gene paralog hybrid haplotype can comprise a gene variant haplotype or a gene paralog variant haplotype.

[0046] In some embodiments, determining the copy number of each of the one or more haplotypes may include determining that the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype given a number of sequence reads of the second plurality of sequence reads each including one or more of the plurality of gene / gene paralog discriminating bases supporting the wildtype gene haplotype. Determining the copy number of each of the one or more haplotypes may include determining that the copy number of the wildtype gene haplotype is 1. Determining that the likelihood of one copy of a wild-type gene haplotype is higher than the likelihood of two copies of a wild-type gene haplotype includes, for each of one or more pairs (or all pairs) of consecutive gene / gene paralog distinguishing bases of a plurality of gene / gene paralog distinguishing bases, where a first haplotype of the one or more haplotypes includes a gene base at consecutive gene / gene paralog distinguishing bases and a second haplotype of the one or more haplotypes includes a gene base and a gene paralog base, or a gene paralog base and a gene base at consecutive gene / gene paralog distinguishing bases, (1) determining whether the likelihood of one copy of a wild-type gene haplotype is higher than the likelihood of two copies of a wild-type gene haplotype. determining that the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype given (1) a number of sequence reads of the second plurality of sequence reads each comprising a gene base and a gene paralog base at consecutive gene / gene paralog discriminating bases, (2) a number of sequence reads of the second plurality of sequence reads each comprising a gene base and a gene paralog base at consecutive gene / gene paralog discriminating bases, (3) a number of sequence reads of the second plurality of sequence reads each comprising a gene base at consecutive gene / gene paralog discriminating bases, and / or (4) a number of sequence reads of the second plurality of sequence reads each comprising a gene paralog base at consecutive gene / gene paralog discriminating bases. The likelihood of one copy of the wildtype gene haplotype may include a sum (e.g., a weighted or unweighted average) of the likelihoods of one copy of the wildtype gene haplotype determined for each of the one or more pairs of consecutive gene / gene paralog discriminating bases.The likelihood of two copies of a wild-type gene haplotype may include the sum (e.g., weighted or unweighted average) of the likelihoods of two copies of a wild-type gene haplotype determined for each of one or more pairs of consecutive gene / gene paralog discriminating bases.

[0047] In some embodiments, the copy number of the wild type gene haplotype is 1. The processor can be programmed with the executable instructions to execute determining that the subject is a carrier of a gene variant haplotype. In some embodiments, the one or more haplotypes include four haplotypes. The total copy number of the gene and gene paralogs can be four. The copy number of each of the four haplotypes can be 1. Determining the gene variant status of the subject can include determining that the subject is a carrier of a gene variant haplotype.

[0048] In some embodiments, the one or more haplotypes include two or more haplotypes. None of the two or more haplotypes can include a gene base at each of the multiple genes / gene paralogues distinguishing bases. Each of the two or more haplotypes can not include a gene base at all of the multiple genes / gene paralogues distinguishing bases. The processor can be programmed with executable instructions to execute: determining that the subject is a compound heterozygote of a gene variant haplotype.

[0049] In some embodiments, the one or more haplotypes comprise only one haplotype. Only one haplotype may not comprise a gene base at each of a plurality of genes / gene paralogues discriminating bases. The processor can be programmed by executable instructions to execute determining that the subject is homozygous for a gene variant haplotype.

[0050] In some embodiments, the first plurality of sequence reads comprises sequence reads that are each about 100 base pairs to about 1000 base pairs in length. In some embodiments, the first plurality of sequence reads comprises paired-end sequence reads and / or single-end sequence reads. In some embodiments, the first plurality of sequence reads are generated by whole genome sequencing (WGS). The WGS can be clinical WGS (cWGS). In some embodiments, the sample comprises cells, cell-free DNA, cell-free fetal DNA, amniotic fluid, a blood sample, a biopsy sample, or a combination thereof.

[0051] The details of one or more implementations of the subject matter described herein are set forth in the accompanying drawings and the description below. Other features, aspects, and advantages will become apparent from the specification, drawings, and claims. Neither this summary nor the following detailed description is intended to define or limit the scope of the inventive subject matter. [Brief description of the drawings]

[0052] [Figure 1] 1 shows an analysis of linkage information between a set of reliable base differences or sites between genes and gene analogs provided by reads and read pairs. [Figure 2A1] 1 shows a non-limiting exemplary detection of a challenging GBA variant by targeted copy number calling and haplotype phasing. [Figure 2A2-1] 1 shows a non-limiting exemplary detection of a challenging GBA variant by targeted copy number calling and haplotype phasing. [Figure 2A2-2] 1 shows a non-limiting exemplary detection of a challenging GBA variant by targeted copy number calling and haplotype phasing. [Figure 2B1] 1 shows a non-limiting exemplary detection of a challenging GBA variant by targeted copy number calling and haplotype phasing. [Figure 2B2] 1 shows a non-limiting exemplary detection of a challenging GBA variant by targeted copy number calling and haplotype phasing. [Figure 2C1] 1 shows a non-limiting exemplary detection of a challenging GBA variant by targeted copy number calling and haplotype phasing. [Figure 2C2] 1 shows a non-limiting exemplary detection of a challenging GBA variant by targeted copy number calling and haplotype phasing. [Figure 3A] 1 shows non-limiting exemplary detection of challenge CYP21A2 variants. [Figure 3B] 1 shows non-limiting exemplary detection of challenge CYP21A2 variants. [Figure 4] FIG. 1 is a flow diagram showing an exemplary method for determining or identifying one or more GBA mutations or GBA mutation status (eg, carrier, compound heterozygous, or homozygous). [Diagram 5] FIG. 1 is a flow diagram showing an exemplary method for determining or identifying one or more CYP21A2 variants or variant status (eg, carrier, compound heterozygote, or homozygote). [Figure 6] FIG. 1 is a flow diagram showing an exemplary method for determining or identifying one or more genetic variants or genetic variant status (e.g., carrier, compound heterozygous, or homozygous). [Figure 7] FIG. 1 is a block diagram of an exemplary computing system configured to determine or identify one or more genetic variants (e.g., GBA variants, CYP21A2 variants) or genetic variant status (e.g., carrier, compound heterozygote, or homozygote).

[0053] Throughout the drawings, reference numbers may be reused to indicate correspondence between referenced elements. The drawings are provided to illustrate example embodiments described herein and are not intended to limit the scope of the present disclosure. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0054] In the following detailed description, reference is made to the accompanying drawings, which form a part of this specification. In the drawings, like symbols typically identify like components unless the context dictates otherwise. The exemplary embodiments described in the detailed description, drawings, and claims are not intended to be limiting. Other embodiments may be utilized, and other changes may be made, without departing from the spirit or scope of the subject matter presented herein. It is readily understood that aspects of the present disclosure, as generally described herein and illustrated in the drawings, can be arranged, substituted, combined, separated, and designed in a wide variety of different configurations, all of which are expressly contemplated herein and make part of the disclosure herein.

[0055] All patents, published patent applications, other publications, and sequences from GenBank and other databases referenced herein are hereby incorporated by reference in their entirety with respect to the relevant art.

[0056] Disclosed herein is a method for determining GBA status. In some embodiments, the method for determining GBA status is under the control of a processor (such as a hardware processor or a virtual processor) and includes receiving a first plurality of sequence reads generated from a sample obtained from a subject. The method may include aligning the first plurality of sequence reads to a reference genome sequence to obtain a second plurality of sequence reads aligned to the GBA gene or the GBAP1 gene in the reference genome sequence. The method may include determining a number of sequence reads of the second plurality of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference genome sequence. The method may include determining a normalized number of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference genome sequence. The method may include determining the total copy number of the GBA gene and the GBAP1 gene using a mixture of Gaussians including a plurality of Gaussians, each of which represents a different integer copy number, given the normalized number of sequence reads aligned to the region between the GBA gene and the GBAP1 gene. The method may include phasing one or more haplotypes derived from the GBA gene or GBAP1 gene in a region of the GBA gene or a corresponding region of the GBAP1 gene that includes a plurality of GBA / GBAP1 discriminatory bases using sequence reads of a second plurality of sequence reads aligned to a region including a plurality of GBA / GBAP1 discriminatory bases. The method may include determining a copy number of each of the one or more haplotypes using a total copy number of the GBA gene and the GBAP1 gene and a number of sequence reads of the second plurality of sequence reads that each include one or more of the plurality of GBA / GBAP1 discriminatory bases that support the haplotype. The method may include determining a GBA status of the subject using one or more haplotypes derived from the GBA gene or GBAP1 gene in a region of the GBA gene or a corresponding region of the GBAP1 gene and / or a copy number of each of the one or more haplotypes.

[0057] Disclosed herein includes a method for determining CYP21A2 status. In some embodiments, the method for determining CYP21A2 status is under the control of a processor (such as a hardware processor or a virtual processor) and includes receiving a first plurality of sequence reads generated from a sample obtained from a subject. The method may include aligning the first plurality of sequence reads to a reference genome sequence to obtain a second plurality of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference genome sequence. The method may include determining the number of sequence reads of the second plurality of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference genome sequence. The method may include determining a normalized number of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference genome sequence. The method may include determining the total copy number of the CYP21A2 gene and the CYP21A1P pseudogene using a Gaussian mixture model including a plurality of Gaussians, each Gaussian representing a different integer copy number, given a normalized number of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene. The method may include phasing one or more haplotypes derived from the CYP21A2 gene or the CYP21A1P pseudogene in a region of the CYP21A2 gene or a corresponding region of the CYP21A1P pseudogene that includes the plurality of CYP21A2 / CYP21A1P discriminatory bases, using sequence reads of a second plurality of sequence reads aligned to a region including the plurality of CYP21A2 / CYP21A1P discriminatory bases or a corresponding region of the CYP21A2 pseudogene. The method may include determining the copy number of each of the one or more haplotypes using the total copy number of the CYP21A2 gene and the CYP21A1P gene, and the number of sequence reads of a second plurality of sequence reads, each of which includes one or more of the plurality of CYP21A2 / CYP21A1P discriminatory bases supporting the haplotype.The method may include determining the subject's CYP21A2 status using one or more haplotypes derived from the CYP21A2 gene or a CYP21A1P pseudogene in a region of the CYP21A2 gene or a corresponding region of the CYP21A1P pseudogene, and / or the copy number of each of the one or more haplotypes.

[0058] Disclosed herein is a system (e.g., a computing system) for determining genetic variants. In some embodiments, the system for determining genetic variants comprises a non-transitory memory configured to store executable instructions and a first plurality of sequence reads generated from a sample obtained from a subject. The system can comprise a processor, such as a hardware processor or a virtual processor, in communication with the non-transitory memory. The processor can be programmed by the executable instructions to align the first plurality of sequence reads to a reference sequence to obtain a second plurality of sequence reads aligned to genes or gene paralogs, or regions therebetween, in the reference sequence. The processor can be programmed by the executable instructions to determine the total copy number of genes and gene paralogs using a mixture of Gaussians including multiple Gaussians, each Gaussian representing a different integer copy number, given the number of sequence reads aligned to genes or gene paralogs, or regions therebetween. The processor can be programmed with the executable instructions to phase one or more haplotypes derived from a gene (including recombinant variants of a gene) or gene paralog, or a region of a gene or a corresponding region of a gene paralog, comprising a plurality of gene / gene paralog discriminating bases, using sequence reads of a second plurality of sequence reads aligned to a region comprising the plurality of gene / gene paralog discriminating bases or a corresponding region. The processor can be programmed with the executable instructions to determine a copy number of each of the one or more haplotypes using the total copy number of the gene and gene paralog and the number of sequence reads of the second plurality of sequence reads that each comprise one or more of the plurality of gene / gene paralog discriminating bases that support the haplotype.

[0059] Detecting variants in short read data Segmental duplications are hotspots of structural variants (e.g., with deletions or duplications) with genetic variants (e.g., gene conversions). Genetic variants can arise when the sequence of a gene is copied to its paralog or vice versa. A paralog of a gene can be a gene or a pseudogene. Segmental duplications can occur for genes with highly homologous gene family members or pseudogenes. Many clinically relevant genes have highly homologous gene family members or pseudogenes and can be affected by segmental duplications. Such clinically relevant genes include genes important in rare diseases, cancer, immunology and pharmacogenetics.

[0060] Analyzing sequence reads (e.g., read alignment and variant calling) of genes that have undergone segmental duplications can be informatively challenging. Such analysis may require combinatorial evaluation of different variants, including single nucleotide polymorphisms (SNPs), insertions and deletions (indels), copy number variations (CNVs), and structural variations (SVs). High sequence similarity of genes and homologous gene family members or pseudogenes can result in poor sequence read alignment and variant calling. Standard secondary analysis pipelines may not yield results or may yield unreliable results. For example, genes and gene paralogs may differ by only a few bases (e.g., 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30 or more bases), making alignment and variant calling difficult. When genes and gene paralogs differ by a few bases, sequence reads that do not contain such bases can be aligned to genes or paralogs with the same or similar alignment scores (e.g., percentage of mismatches). As a result, it is difficult to align reads to genes and paralogs, resulting in low alignment quality (e.g., MapQ quality). As a result, variant calling using such low-quality read alignments can be difficult.

[0061] Targeted CNV calling of genes and gene paralogs can be useful. The total copy number of genes and paralogs can be determined by read counting and normalization using all reads. Population depth distribution can be modeled to call copy number. Genes and gene paralogs can be distinguished and their copy number can be determined using sequence read counts at fixed base differences. Gene fusions can be determined based on, for example, changes in gene copy number at fixed base differences. Targeted calling of SNPs and indels can be performed. Star alleles can be called based on all called variants and assigned to haplotypes. One, some, most, or all of these processes can be performed sequentially by one or more computing systems, or in parallel by multiple computing systems.

[0062] The methods disclosed herein can be used to detect short genetic modifications, such as gene conversions, where the sequence of a gene is mutated to become identical to the sequence of another gene (e.g., a paralog of a gene). Each mutated sequence can be as small as a single base. The methods can allow for the detection of single base gene conversions by haplotype phasing. Genetic modification variants (reciprocal or non-reciprocal) of the GBA gene can result in Gaucher disease, which has an incidence of 1 in 50,000 to 100,000. Heterozygotes are associated with Parkinson's disease. Genetic modification variants of the CYP21A2 gene can result in 21-hydroxylase deficiency congenital adrenal hyperplasia (21-OHD CAH), which has an incidence of 1:10,000 to 1:16,000.

[0063] A gene of interest and a paralog (e.g., pseudogene) of the gene of interest are referred to herein as gene A and gene B. The presence of sequence reads (e.g., short sequence reads) of gene B in the reference genome can make alignment and variant calling difficult. Often, a small stretch of sequence (or a single base) of gene A can be mutated, for example, by recombination, such as gene conversion, to appear identical to the corresponding sequence in gene B. This type of recombination or gene conversion variant can be extremely difficult to detect, since variants containing reads of gene A may align to gene B instead of gene A.

[0064] Calling of carrier samples or compound heterozygotes. Based on a set of reliable base differences or sites between gene A and gene B, haplotype phasing can be performed. Based on a set of reliable base differences or sites between gene A and gene B, a caller for calling or identifying recombination or gene conversion variants can analyze linkage information between these difference sites provided by the reads and read pairs. Linkage information can be analyzed, for example, by read backphasing. One read or read pair covering site n and site m can indicate whether the haplotype from which the read or read pair is derived has a gene A or B base at site n and a gene A or B base at site m. The caller can then phase all haplotypes derived from either gene A or gene B in the region with a set of reliable base differences, and identify gene A and gene B haplotypes as well as one or more hybrid haplotypes (a mixture of gene A and gene B bases on the same haplotype). With reference to FIG. 1 , a single read pair covering site 1 and site 4 can indicate that the haplotype from which the read pair is derived (haplotype x) has a gene A base at site 1 and a gene A base at site 4. A single read pair covering site 3 and site 5 can indicate that the haplotype from which the read pair is derived (haplotype y) has a gene A base at site 3 and a gene B base at site 5. A single read covering site 4 and site 5 can indicate that the haplotype from which the read is derived (haplotype y) has a gene A base at site 4 and a gene B base at site 5. The caller can phase all haplotypes from either gene A or gene B in regions with a set of five credible base differences and identify gene A haplotypes (haplotype 1) and gene B haplotypes (haplotype 2) as well as hybrid haplotypes (haplotypes 3 and 4). The number of haplotypes, the number of sites, and sequence read coverage of the sites (e.g., sites 1 and 4, or sites 3 and 5) are shown in FIG. 1 for illustrative purposes only and are not intended to be limiting.

[0065] To assess the relative abundance of different haplotypes, the caller can call the CN of each haplotype using the total copy number (CN) of gene A and gene B and the haplotype supporting read count at the discriminating base. The total CN of gene A and gene B can be determined using a Gaussian mixture model. Determining CN using a Gaussian mixture model by SMN callers and CYP2D6 callers is described in PCT Publication No. 2021 / 045947, entitled METHODS AND SYSTEMS FOR DIAGNOSING FROM WHOLE GENOME SEQUENCING DATA, the contents of which are incorporated herein by reference in their entirety. The caller can compare two scenarios: one copy of wild-type gene A haplotype versus two copies of wild-type gene A haplotype. The caller can determine which scenario is more likely given the number of supporting reads in the data. If the caller calls only one copy of the wild-type gene A haplotype, this indicates that the individual is a carrier of the disease-causing variant. If an individual is a carrier of more than one variant haplotype, and no haplotypes carrying gene A bases exist at all variant sites of interest, the caller calls the sample as a compound heterozygote with no copies of the wild-type gene A.

[0066] Calling samples homozygous for variants. Based on the list of gene conversion variants to be detected, the caller can call the copy number (CN) of the gene A base. The number of reads supporting the gene A base or gene B base, and the total CN of gene A and gene B can be used to determine the most likely combination of CN of gene A base and CN of gene B base. If the CN of the gene A base is called 0, this indicates that the individual has no copies of wild type gene A (a haplotype carrying a gene A base at the mutation site of interest) and is homozygous for the gene conversion variant.

[0067] GBA mutants Detection of challenging GBA variants in short-read WGS data The sequence homology between GBA and GBAP1 results in poorer mapping quality (Figure 2A1) and less accurate variant calling by standard secondary analysis pipelines. In addition, recombination variants in which GBA bases are mutated to corresponding bases in GBAP1 are difficult to detect because variant reads align to GBAP1. Therefore, a novel WGS-based bioinformatics method for calling GBA variants, called Gauchian, was developed using 2405 WGS datasets from the 1000 Genomes Project (1kGP).

[0068] Gauchian analysis begins with determining copy number changes. Reciprocal recombination across the homologous region results in copy number gain (duplication) or loss (fusion) of the 20.6 kb region between the two genes. Because the breakpoints can vary in location, Gauchian uses sequencing depth in the unique region between the two genes to detect copy number variation (CNV) (Figure 2A1, Figure 2A2, Figure 2B1, and Figure 2B2). Of the 2504 1 kGP samples, 108 samples were found to have reciprocal recombination (15 deletions and 93 duplications). Duplications can result in GBA-GBAP1 fusions, but always leave two intact copies of GBA, whereas fusions can result in GBA variants if the deletion breakpoint is contained within the GBA gene (GBA-GBAP1 fusions). After determining the copy number changes, Gauchian used the base differences between the homologous regions of GBA and GBAP1 shown in Table 1 to distinguish the two genes and identify the breakpoints of the CNVs (Figure 2A2). The major homologous region, exons 9-11, is where pathogenic deletions are most likely to occur. Thus, Gauchian site-phased both the GBA and GBAP1 haplotypes across the exon 9-11 homology region to further elucidate the breakpoints (Figures 2C1 and 2C2). For the 1kGP deletion sample, Gauchian identified breakpoints that did not alter the GBA gene in all but one sample. The breakpoints were contained in a region that is mostly identical between GBA and GBAP1, extending from the 3'UTR beyond the gene (0 mapping quality region in Figure 2A1). Such CNVs are known and are likely benign because they leave the GBA gene, or at least the coding region, intact. In one deletion sample, Gauchian identified a breakpoint in exons 9-11 that created a RecNciI fusion (Figures 2C1 and 2C2).

[0069] [Table 1]

[0070] In addition to CNVs resulting from reciprocal recombination, Gauchian performed targeted calling of pathogenic GBA variants, including simple small variants as well as gene conversions, and challenged SNVs in exons 9-11 in which GBA bases are mutated to corresponding bases in GBAP1, also likely arising via small gene conversion events. These include p.L483P, p.D448H, c.1263del (55 bp deletion), RecNciI (containing 3 SNVs p.L483P, A495P and Val499=), RecTL (containing RecNciI and p.D448H) and c.1263del+RecTL (containing RecNciI, p.D448H and c.1263del) (Figures 2C1 and 2C2). The high homology and frequent gene conversions between GBA and GBAP1 make exons 9-11 a very challenging region for standard secondary analysis pipelines. For example, three positions in the GBAP1 reference sequence in hg38 erroneously contain GBA bases (Figure 2C1), so GBAp.L483P reads would easily align to GBAP1, causing a false-negative call. In addition, there are GBAP1 haplotypes in the population that are partially converted to GBA, and those converted bases would lead GBAP1 reads to align to GBA, causing false-positive GBA variant calls at nearby positions (Figure 2C1, purple shading / bottom two rows). Given all reads that align to either GBA or GBAP1, Gauchian is able to phase haplotypes across the entire homology region and thus accurately call these variants. This also enabled the identification of larger gene conversion events such as RecTL or c.1263del+RecTL, which standard pipelines would miss because the variant reads align to GBAP1. Gauchian detected 42 samples with simple non-recombinant minor mutations in 1 kGP, as well as 5 samples with p.L483P, 2 samples with c.1263del, and 2 samples with c.1263del+RecTL conversion.

[0071] Figure 2A1. Median mapping quality (red line) across 2504 1kGP samples plotted for each position in the GBA / GBAP1 region (hg38). The median filter is applied to a 50 bp window. The 11 exons of GBA are shown as orange boxes. GBAP1 and MTX1 exons are shown as green and purple boxes, respectively. The 4 kb major homologous region (98.1% sequence similarity, exons 9-11) between GBA and GBAP1 is shaded in pink in Figure 2A1 (corresponding to the green boxed region in Figure 2A2), highlighting the region of low mapping quality. The light blue box indicates the 10 kb unique region between the two genes where copy number calling is performed in Gauchian. Figure 2B1. Distribution of normalized depth in the 10 kb CN calling region in 2504 1kGP samples shows peaks at 1 (deletion), 2, and 3-8 (duplication). This number plus 2 was the total copies of both GBA and GBAP1 combined. The 10 kb unique region is a proxy for the 20.6 kb region between GBA and GBAP1 that would be lost or gained due to reciprocal recombination (the breakpoints can vary). Normally, a diploid sample has a total of 4 copies of GBA+GBAP1 (2 copies of each) and 2 copies of the 10 kb unique region. A deletion between GBA and GBAP1 results in a loss of one copy of the 10 kb region (now CN 1) and a loss of one copy of GBA+GBAP1 (now CN 3). Similarly, a duplication results in a gain of one copy of the 10 kb region (now CN3 here) and a gain of one copy of GBA+GBAP1 (now CN5 here). Thus, the CN (GBA+GBAP1) is 2 more than the CN of the 10 kb region.

[0072] Figure 2B2 shows the race-dependent distribution of CN. Figure 2C1. Recombinant haplotypes in the exon 9-11 homology region are differentiated by the GBA / GBAP1 discriminatory base (x-axis). The reference genome sequence is shaded in yellow. There is an error in the hg38 reference, where the first three sites of GBAP1 show GBA bases, which may result in alignment errors. GBA recombinant haplotypes, including those in which one or several adjacent sites are mutated to GBAP1 bases, resulting from either gene conversion or fusion deletion, are shown on a white background. The grey bases indicate that the bases can be either GBA or GBAP1, depending on the breakpoint location of the fusion / conversion. Two exemplary GBAP1 haplotypes that are partially converted to GBA, causing false-positive GBA variant calls, are shaded in purple. For the first example, specifically for hg38, where the first three sites are wrong for the GBAP1 reference, the reverse L483P variant on GBAP1 directs the aligner to align GBAP1 reads to GBA, causing a nearby A495P FP call. In the second example, the reverse-c.1263del variant inserts 55bp into GBAP1, causing GBAP1 reads to align to GBA, causing a nearby D448H FP call.

[0073] Detection of all classes of GBA variants is possible using ONT long-read sequencing. Cross-validation To validate Gauchian, ONT sequencing (Oxford Nanopore Technologies) was performed on 14 samples in which Gauchian detected reciprocal recombinations (11 with copy number gains and 3 with copy number losses), 2 samples in which Gauchian detected non-reciprocal recombination c.1263del+RecTL, 9 samples carrying SNVs (8 carrying p.L483P and 1 A456P) and 12 GBA negative controls. These samples were selected from 1kGP and AMP-PD and included 2 samples in which GATK missed p.L483P and 2 samples in which GATK miscalled the p.A456P variant. In all cases, including the cases in which GATK results differed, the ONT and Gauchian results were concordant.

[0074] Furthermore, because Gauchian showed evidence of multicopy CNV, we used digital PCR to precisely quantify the copy number of the 20.6 kb region involved in recombination in the four samples in which Gauchian detected copy number gain (increased copy numbers: 1, 3, 5, and 6, respectively). The dPCR results were as expected, and the CNs were consistent with those detected by Gauchian and ONT (Table 2).

[0075] [Table 2]

[0076] Prevalence of GBA recombining and non-recombining variants in healthy, PD and LBD populations Having validated Gauchian, we applied it to Parkinson's disease (PD) and Lewy body dementia (LBD) cohorts from AMP-PD to provide the first large-scale analysis of GBA recombination variants and to estimate the prevalence of different GBA mutations in healthy, PD, and LBD populations.

[0077] For CNVs with breakpoints that do not alter the GBA gene (copy number gains and non-pathogenic fusion alleles), there was no enrichment in PD or LBD cases versus controls (Table 3). Among Caucasians, CNVs were found in 18 of 2234 PD cases (0.81%) (10 duplications and 8 deletions) and 14 of 1214 controls (1.15%) (7 duplications and 7 deletions) (p-value=0.35, Fisher's exact test). Among Africans, CNVs (duplications) were found in 2 of 25 PD cases and 3 of 33 controls (p-value=1). Among Caucasians, CNVs were found in 34 of 2598 LBD cases (1.31%) (21 duplications and 13 deletions) and 24 of 1941 controls (1.24%) (11 duplications and 13 deletions) (p-value=0.89). Across both the 1kGP and PD cohorts, CNVs (especially duplications) were overall more than nine times more frequent in Africans than in Caucasians (1kGP: 11.6% vs. 0.6%; PD: 8.6% vs. 0.9%). These results are consistent with recent evidence in African genomes showing unexplored structural variants and a larger, still largely unexplored, genetic diversity in Africans. Interestingly, 3 of 10 PD cases with duplications harbor a second pathogenic GBA variant, and 4 of 21 LBD cases with duplications harbor a second pathogenic GBA variant. Although duplications do not alter the GBA gene itself, duplications may result in a higher chance of acquiring a second GBA variant. This is consistent with what was found in the RAPSODI and QSBB ONT data.

[0078] [Table 3]

[0079] In addition to benign CNVs, recombination variants were detected in all three cohorts (Table 4). GBA recombination variants are more common in LBD than in PD (50 / 2598, 1.92% vs. 19 / 2325, 0.82%, p-value=0.0009). GATK variant calling was available for PD and LBD samples from AMP-PD. Due to sequence homology in exons 9-11, GATK undercalled all recombination variants except D448H. For D448H, GATK called two false positives due to a GBAP1 haplotype in which the base was converted to GBA (see Figure 2C). For all PD and LBD case+control populations, GATK called 35 recombination variants and Gauchian called 77, a more than double variant call.

[0080] Gauchian also detected simple non-recombinant variants in three cohorts (Table 5). Again, GBA variants are more common in LBD than in PD.

[0081] [Table 4]

[0082] [Table 5]

[0083] Gauchian-WGS based GBA caller In Gauchian, the WGS-based GBA caller disclosed herein uses a novel approach to overcome this challenge based on a strategy to resolve closely related paralogs as described in the SMN1 / SMN2 caller (Chen et al. Spinal muscular atrophy diagnosis and carrier screening from genome sequencing data, Genet Med 22, 945-953 (2020), the contents of which are incorporated herein by reference in their entirety) and in the Cyrius CYP2D6 caller (Chen et al., Cyrius: accurate CYP2D6 genotyping using whole genome sequencing data, Pharmacogenomics J 21, 251-261 (2021), the contents of which are incorporated herein by reference in their entirety). In some embodiments, the Gauchian method can be applied to sequence reads from targeted sequencing, such as sequencing of 5, 10, 20, 30, 40, 50, 100, 200, or more genes.

[0084] First, Gauchian calculates the copy number of a 10 kb unique region between GBA and GBAP1 (chr1:155220429-155230539, hg38) following a similar targeted CNV calling method used by SMN1 / SMN2 callers and CyriusCYP2D6 callers. The number of reads aligned to this region was normalized and corrected for GC content, and the copy number was called from a Gaussian mixture model. Deviation of this copy number (CN) from the expected copy number (CN) of 2 indicates the presence of a CNV. For example, 1 copy indicates a deletion and 3 copies indicates a duplication. Thus, adding 2 to this number gives the total copies of both GBA and GBAP1 combined (compare Figure 2B1 and Figure 2B2). The total copies of both GBA and GBAP1 combined is abbreviated herein as CN(GBA+GBAP1).

[0085] Next, Gauchian identifies the CNV breakpoints following a similar approach as used by the CyriusCYP2D6 caller. To do this, it uses 82 reliable bases that differ between GBA and GBAP1. Gauchian estimates the GBACN at each of the 82 GBA / GBAP1 discriminatory base positions based on the CN (GBA+GBAP1) and the number of reads that support the GBA and GBAP1 specific bases. CNV breakpoints are identified when the CN of the GBA changes. For example, a switch from CN1 to CN2 indicates a breakpoint of a deletion, and a switch from CN3 to CN2 indicates a breakpoint of a duplication. The exact breakpoints are further refined by haplotype phasing, as described in the next paragraph.

[0086] To identify recombination variants, Gauchian analyzes a 1.1 kb region (Figure 2C) that contains the key GBA / GBAP1 recombination variants (p.L483P, p.D448H, c.1263del, RecNciI, RecTL, and c.1263del+RecTL). This region contains 10 GBA / GBAP1 base differences. Based on the reads and read pairs, Gauchian phases all haplotypes derived from either GBA or GBAP1 in this region and identifies hybrid haplotypes (i.e., a mixture of GBA and GBAP1 bases on the same haplotype). To assess the relative abundance of different haplotypes, Gauchian uses haplotype-supported read counts at the discriminating base to call the CN (GBA+GBAP1) as well as the CN of each haplotype. Gauchian compares two scenarios: one copy of the wild-type GBA haplotype versus two copies of the wild-type GBA haplotype. Gauchian determines which scenario is more likely given the number of supporting reads in the data. If Gauchian calls only one copy of the wild-type GBA haplotype, this indicates that the individual is a carrier of the disease-causing variant. If an individual is a carrier of two or more variant haplotypes and there is no haplotype with the GBA base at all variant sites of interest, Gauchian calls the sample as compound heterozygous. A homozygous variant is called if the CN of the GBA base is called 0. Finally, for simple small variants, Gauchian analyzes the read alignment and calls the CN of the variant used by the SMN1 / SMN2 caller and the CyriusCYP2D6 caller.

[0087] CYP21A2 variants The RCCX module of the human MHC class III region contains approximately 30 kb of tandem repeats with 99.6% similarity. The RCCX module encodes RP1, C4A / B, CYP21A2, and TNXB. CYP21A1P is a pseudogene of CYP21A2. Recombinant variants of CYP21A2 can cause 21-hydroxylase deficiency congenital adrenal hyperplasia (21-OHD CAH) with an incidence of 1:10,000 to 1:16,000 live births. C4A and C4B together form complement component 4 (C4). C4 deficiency is associated with autoimmune diseases such as lupus. Mutations in TNXB can cause Ehlers-Danlos syndrome.

[0088] The RCCX repeats in the subject's samples were as described above. The CN of the RCCX repeats in the samples was determined using a Gaussian mixture model. Figure 3A shows the race-dependent distribution of CN. The deletion breakpoints were identified by examining the switches in the CN of the discriminating SNP sites. C4A and C4B contain five SNPs that mark the functional difference between C4A / C4B. Reads and read pairs aligned to 14 differences between CYP21A2 and CYP21A1P, including the nine recombination variants shown in Figure 3B, were analyzed with read backphasing for haplotype phasing. Whether the CYP21A2 wild-type haplotype was present in only one copy in the subject's sample was tested based on the depth / number of supporting reads in the data. Figure 3C shows the distribution of CYP21A2 haplotypes other than the CYP21A2 wild-type haplotype that the subject had. The CYP21A2 wild-type haplotype is represented as 11111111111111, where all bases are CYP21A2 bases at position 14. The CYP21A2 haplotype shown in Figure 3C is not a CYP21A2 wild-type haplotype and would therefore be represented, for example, by 11111111111121, indicating that the 13th base is a CYP21A1P base and the remaining bases are CYP21A2 bases.

[0089] Determination of GBA mutants and mutant status FIG. 4 is a flow diagram illustrating an exemplary method 400 for determining or identifying one or more GBA variants or GBA variant states. The method 400 may be embodied in a set of executable program instructions stored on a computer-readable medium, such as one or more disk drives of a computing system. For example, a computing system 700, illustrated in FIG. 7 and described in more detail below, may execute a set of executable program instructions to perform the method 400. When the method 400 is initiated, the executable program instructions may be loaded into a memory, such as a RAM, and executed by one or more processors of the computing system 700. Although the method 400 is described with respect to the computing system 700 illustrated in FIG. 7, the description is merely exemplary and is not intended to be limiting. In some embodiments, the method 400, or portions thereof, may be executed serially or in parallel by multiple computing systems.

[0090] After the method 400 starts at block 404, the method 400 proceeds to block 408, where a computing system (e.g., computing system 700 described with reference to FIG. 7) aligns the first plurality of sequence reads to a reference sequence (e.g., a reference genome sequence such as hg19 or hg38) to obtain a second plurality of sequence reads aligned to the GBA gene or GBAP1 gene in the reference sequence (including an alignment of each of the second plurality of sequence reads to the GBA gene or GBAP1 gene in the reference sequence). The computing system can receive the first plurality of sequence reads generated from a sample obtained from the subject. The computing system can store the first plurality of sequence reads in a memory. The computing system can load the first plurality of sequence reads into the memory. The sequence reads can be generated by a technique such as sequencing by synthesis, sequencing by ligation, or sequencing by ligation. Sequence reads can be generated using instruments such as the MINISEQ, MISEQ, NEXTSEQ, HISEQ, and NOVASEQ sequencing instruments from Illumina, Inc. (San Diego, Calif.).

[0091] The sequence reads can be, for example, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1250, 1500, 1750, 2000 or more base pairs (bps) in length, respectively. For example, the sequence reads are about 100 base pairs to about 1000 base pairs in length, respectively. The sequence reads can include paired-end sequence reads. The sequence reads can include single-end sequence reads. The sequence reads can be generated by whole genome sequencing (WGS). The WGS can be clinical WGS (cWGS). The sequence reads can include single-end sequence reads. The sequence reads can be generated by targeted sequencing, such as sequencing of 5, 10, 20, 30, 40, 50, 100, 200 or more genes. The sample can include cells, cell-free DNA, cell-free fetal DNA, amniotic fluid, a blood sample, a biopsy sample, or a combination thereof.

[0092] The sequence reads can be aligned to the GBA or GBAP1 gene in the reference sequence with an alignment quality score of greater than or equal to 0. The sequence reads can be aligned to the GBA or GBAP1 gene in the reference sequence with an alignment quality score of about 0 (e.g., when sequences are aligned to regions where genes and gene paralogs are highly homologous). The computing systems used were Burrows-Wheeler Aligner (BWA), iSAAC, BarraCUDA, BFAST, BLASTN, BLAT, Bowtie, CASHX, Cloudburst, CUDA-EC, CUSHAW, CUSHAW2, CUSHAW2-GPU, drFAST, ELAND, ERNE, GNUMAP, GEM, GensearchNGS, GMAP and GSNAP, Geneious Assembler, LAST, MAQ, mrFAST and mrsFAST, MOM, MOSAIK, MPscan, Novoaligh and NovoalignCS, NextGENe, Omixon, PALMapper, Partek, PASS, PerM, PRIMEX, QPalma, RazerS, REAL, cREAL, RMAP, rNA, RT Sequence reads can be aligned to a reference sequence using an aligner or alignment method such as Investigator, Segemehl, SeqMap, Shrec, SHRiMP, SLIDER, SOAP, SOAP2, SOAP3 and SOAP3-dp, SOCS, SSAHA and SSAHA2, Stampy, SToRM, Subread and Subjunc, Taipan, UGENE, VelociMapper, XpressAlign and ZOOM.

[0093] The method 400 proceeds from block 408 to block 412, where the computing system determines a number (e.g., normalized and / or corrected number) of sequence reads of the second plurality of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference sequence. The unique region between the GBA gene and the GBAP1 gene in the reference sequence may comprise a unique region of about 10 kilobases in length. The unique region between the GBA gene and the GBAP1 gene in the reference sequence may comprise chr1:155220429-155230539 of hg38 or a corresponding region of the reference human genome sequence.

[0094] The computing system can determine a normalized number of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference sequence. The computing system can determine the normalized number of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference sequence using (1a) a depth of the sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene, (1b) a length of the unique region, (2a) a depth of the sequence reads of a first plurality of sequence reads aligned to each of a plurality of regions of the reference sequence other than the locus including the GBA gene and the GBAP1 gene, and / or (2b) a length of each of a plurality of regions of the reference other than the locus including the GBA gene and the GBAP1 gene.

[0095] The computing system can determine the normalized, corrected number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference sequence from the normalized number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference sequence. To determine the normalized, corrected number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference sequence, the computing system can determine the normalized, GC content-corrected number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference sequence from the normalized number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference sequence. The computing system can determine the normalized, GC content-corrected number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference sequence from the normalized number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference sequence using (1) the GC content of the unique region between the GBA gene and the GBAP1 gene. The computing system can determine the normalized GC content corrected number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference sequence from the normalized number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference sequence using (1) the GC content of the unique region between the GBA gene and the GBAP1 gene and / or (2) the GC content of one or more regions of the reference sequence other than the locus containing the GBA gene and the GBAP1 gene (or one or more regions of the reference sequence not containing the GBA gene and the GBAP1 gene). For example, the computer system can determine the normalized GC content corrected number of sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference sequence using (1) the GC content of the unique region between the GBA gene and the GBAP1 gene, and (2) the GC content of the region of the reference sequence other than the locus containing the GBA gene and the GBAP1 gene.As another example, the computer system can determine a normalized GC content corrected number of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in a reference sequence using (1) the GC content of the unique region between the GBA gene and the GBAP1 gene, and (2) the GC content of multiple regions (e.g., 2, 3, 4, 5, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500, 1000, 2000, 3000, 4000, 5000, 10000, or more regions) of the reference sequence other than the locus including the GBA gene and the GBAP1 gene.

[0096] Method 400 proceeds from block 412 to block 416, where the computing system determines the total copy number of the GBA and GBAP1 genes using a Gaussian mixture model including multiple Gaussians, each representing a different integer copy number, given the number of sequence reads aligned between the GBA and GBAP1 genes (e.g., normalized and / or corrected sequence reads). The computing system can determine the total copy number of the GBA and GBAP1 genes using a Gaussian mixture model, given the normalized number of sequence reads aligned to the region between the GBA and GBAP1 genes. The computing system can determine the total copy number of the GBA and GBAP1 genes using a Gaussian mixture model, given the normalized, corrected number of sequence reads aligned to the region between the GBA and GBAP1 genes.

[0097] The total copy number can be, for example, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more. The Gaussian mixture model can include a one-dimensional Gaussian mixture model. The Gaussians of the Gaussian mixture model can represent integer copy numbers, for example, 0 to 5, 0 to 6, 0 to 7, 0 to 8, 0 to 9, 0 to 10, 0 to 11, 0 to 12, 0 to 13, 0 to 14 or 0 to 15. For example, the Gaussians of the Gaussian mixture model can represent integer copy numbers 0 to 10. The average of each of the multiple Gaussians can be the integer copy number represented by the Gaussian. The average of each of the multiple Gaussians can be the integer copy number represented by the Gaussian (for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 or more copies). The standard deviation of the Gaussians can be, for example, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 or more, or can be about 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 or more. The multiple Gaussians of the Gaussian mixture model can include, for example, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, or more Gaussians. For example, the multiple Gaussians of the Gaussian mixture model can include 5 Gaussians.

[0098] To determine the total copy number of the GBA and GBAP1 genes, a computing system can use a Gaussian mixture model to determine the copy number of the region between the GBA and GBAP1 genes given the normalized number of sequence reads aligned to the region between the GBA and GBAP1 genes. The total copy number of the GBA and GBAP1 genes can be the copy number of the region between the GBA and GBAP1 genes plus 2.

[0099] The computing system can determine the total copy number of the GBA gene and the GBAP1 gene using a Gaussian mixture model and a predetermined posterior probability threshold given the normalized number of sequence reads aligned to the region between the GBA gene and the GBAP1 gene, for example, the predetermined posterior probability threshold can be greater than or equal to 0.70, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.80, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.90, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99. 0.70, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.80, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.90, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99 or greater. For example, the predetermined posterior probability threshold is 0.95.

[0100] Method 400 proceeds from block 416 to block 420, where the computing system uses sequence reads of the second plurality of sequence reads aligned to a region including a plurality of GBA / GBAP1 discriminating bases or a corresponding region to phase one or more haplotypes derived from the GBA gene or the GBAP1 gene in a region of the GBA gene or a corresponding region of the GBAP1 gene that includes a plurality of GBA / GBAP1 discriminating bases (or positions or sites of discriminating bases). For example, the sequence reads can be aligned to the reference sequence such that the sequence reads overlap with the GBA / GBAP1 discriminating bases (or sites of the GBA / GBAP1 discriminating bases) or such that the bases of the sequence reads are aligned to the GBA / GBAP1 paralog discriminating bases (or sites of the GBA / GBAP1 paralog discriminating bases). The sequence reads of the second plurality of sequence reads can be aligned to a region of the GBA gene including a plurality of GBA / GBAP1 discriminating bases or a corresponding region of the GBAP1 gene with an alignment quality score of 0 or greater.

[0101] The one or more haplotypes include a wild-type GBA haplotype, a wild-type GBAP1 haplotype, and / or a GBA / GBAP1 hybrid haplotype. The GBA / GBAP1 hybrid haplotype may include both GBA and GBAP1 bases. The GBA / GBAP1 hybrid haplotype may be a recombination variant. The GBA / GBAP1 hybrid haplotype may include a GBA mutant haplotype or a GBAP1 mutant haplotype. The haplotype may include a reciprocal recombination variant. The haplotype may include a non-reciprocal recombination variant or a gene conversion variant. The reference sequence may include a reference genome sequence.

[0102] To phase one or more haplotypes derived from the GBA gene or GBAP1 gene, the computing system can use sequence reads of the second plurality of sequence reads aligned to a region including or corresponding to the plurality of GBA / GBAP1 discriminating bases to analyze linkage information between the GBA / GBAP1 discriminating bases of the plurality of GBA / GBAP1 discriminating bases. The computing system can use sequence reads of the second plurality of sequence reads aligned to two or more of the plurality of GBA / GBAP1 discriminating bases, respectively, to phase one or more haplotypes derived from the GBA gene or GBAP1 gene. For example, referring to FIG. 1, assuming that gene A and gene B shown in the figure are GBA gene or GBAP1 gene, respectively, one read pair covering site 1 and site 4 of the discriminating base can indicate a haplotype (haplotype x) derived from the read pair having the GBA gene base at site 1 and the GBA gene base at site 4. One read pair covering site 3 and site 5 of the distinguishing bases can indicate a haplotype (haplotype y) resulting from the read pair having a GBA gene base at site 3 and a GBAP1 gene base at site 5. One read covering site 4 and site 5 of the distinguishing bases can indicate a haplotype (haplotype y) resulting from the read having a GBA gene base at site 4 and a GBAP1 gene base at site 5. The computing system can phase all haplotypes resulting from either the GBA gene or the GBAP1 gene in a region having a set of five credible base differences and identify haplotypes of the GBA gene (haplotype 1) and haplotypes of the GBAP1 gene (haplotype 2) as well as hybrid haplotypes (haplotypes 3 and 4). The number of haplotypes, bases of the haplotypes at sites, number of sites, and site (e.g., sites 1 and 4, or sites 3 and 5) sequence read coverage are shown in FIG. 1 for illustration only and are not intended to be limiting.

[0103] The region of the GBA gene or the corresponding region of the GBAP1 gene that includes the plurality of GBA / GBAP1 discriminator bases may be about 1.1 (or 0.8, 0.9, 1, 1.2, 1.3, or more) kilobases in length. The region of the GBA gene or the corresponding region of the GBAP1 gene that includes the plurality of GBA / GBAP1 discriminator bases may include exons 9-11 of the GBA gene or the GBAP1 gene, respectively. The region of the GBA gene or the corresponding region of the GBAP1 gene that includes the plurality of GBA / GBAP1 discriminator bases may include p.L483P, p.D448H, c.1263del, RecNciI, RecTL, and c.1263del+RecTL. The plurality of GBA / GBAP1 discriminator bases may include 10 GBA / GBAP1 discriminator bases.

[0104] From block 420, method 400 proceeds to block 424, where the computing system determines a copy number for each of the one or more haplotypes using the total copy number of the GBA and GBAP1 genes and the number of sequence reads of the second plurality of sequence reads that each include one or more of the plurality of GBA / GBAP1 discriminating bases that support the haplotype. The copy number of the haplotype can be, for example, 1, 2, 3, 4 or more.

[0105] The computing system can determine the subject's GBA status (e.g., carrier, compound heterozygote, or homozygote) using one or more haplotypes from the GBA gene or the GBAP1 gene and / or copy numbers of each of the one or more haplotypes in a region of the GBA gene or a corresponding region of the GBAP1 gene. The computing system can generate a user interface (UI), such as a graphical user interface, that includes UI elements that represent or include the GBA status. The UI can include the GBA status as part of the UI element. The UI element can be a window (e.g., a container window, a browser window, a text terminal, a child window, or a message window), a menu (e.g., a menu bar, a context menu, or a menu extra), an icon, or a tab. The UI element can be for an input control (e.g., a check box, a radio button, a drop-down list, a list box, a button, a toggle, a text field, or a date field). The UI element can be for navigation (e.g., a breadcrumb, a slider, a search field, pagination, a slider, a tag, an icon). A UI element can provide information (e.g., a tooltip, an icon, a progress bar, a notification, a message box, or a modal window). A UI element can be a container (e.g., an accordion).

[0106] Carrier. To determine the copy number of each of the one or more haplotypes, the computing system can determine that the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype given the number of sequence reads of the second plurality of sequence reads, each of which includes one or more of the plurality of GBA / GBAP1 discriminating bases supporting the wild-type GBA haplotype. The computing system can determine that for each of one, one or more (e.g., two, three, or four), or one or more haplotypes (e.g., using the sequence reads thereof), the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype. The likelihood difference can be, for example, 1%, 2%, 3%, 5%, 10%, 15%, 20%, or more. The computing system can determine that the copy number of the wild-type GBA haplotype is 1. The computer system can determine that the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype for each of one or more pairs (or all pairs) of a plurality of consecutive GBA / GBAP1 discriminatory bases of the GBA / GBAP1 discriminatory bases, where a first haplotype of the one or more haplotypes includes a GBA base in the consecutive GBA / GBAP1 discriminatory bases and a second haplotype of the one or more haplotypes includes a GBA base and a GBAP1 base (or a GBAP1 base and a GBA base) in the consecutive GBA / GBAP1 discriminatory bases. The first haplotype can include a GBA base in the consecutive GBA / GBAP1 discriminatory bases. The second haplotype can include a conversion of a GBA base to a GBAP1 base (or a conversion of a GBAP1 base to a GBA base) between the consecutive GBA / GBAP1 discriminatory bases. Contiguous GBA / GBAP1 discriminator bases are contiguous within a plurality of GBA / GBAP1 discriminator bases, regardless of whether the contiguous GBA / GBAP1 discriminator bases are adjacent bases in the reference sequence.For example, the GBA / GBAP1 discriminating bases at positions 5 and 6 in the examples below are consecutive GBA / GBAP1 discriminating bases, regardless of whether the consecutive GBA / GBAP1 discriminating bases are adjacent bases in the reference sequence. The computing system may combine (e.g., average or weighted average) the likelihoods determined for each of one or more pairs (or all pairs) of consecutive GBA / GBAP1 discriminating bases to determine that the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype.

[0107] The computing system can determine that (1) given a number of sequence reads of the second plurality of sequence reads, each of which includes a GBA base in consecutive GBA / GBAP1 discriminating bases, for each of one or more pairs (or all pairs) of consecutive GBA / GBAP1 discriminating bases, the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype. The computing system can determine that, given (2) the number of sequence reads of the second plurality of sequence reads that each comprise a GBA base and a GBAP1 base, (or a GBAP1 base and a GBA base), of consecutive GBA / GBAP1 discriminating bases, and / or (3) the number of sequence reads of the second plurality of sequence reads that each comprise a GBAP1 base and a GBA base, of consecutive GBA / GBAP1 discriminating bases, the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype for each of one or more pairs (or all pairs) of consecutive GBA / GBAP1 discriminating bases. In some embodiments, the computing system can determine that (4) given a number of sequence reads of the second plurality of sequence reads, each of which comprises a GBAP1 base in consecutive GBA / GBAP1 discriminating bases, for each of one or more pairs (or all pairs) of consecutive GBA / GBAP1 discriminating bases, the likelihood of one copy of the wildtype GBA haplotype is higher than the likelihood of two copies of the wildtype GBA haplotype.

[0108] For example, by analyzing the linkage information between the GBA / GBAP1 discriminator bases, the following haplotypes can be determined for a subject (for illustrative purposes only, at 6 bases or sites or positions of the discriminator base):

[0109] [Table 6] Conversion of the GBA gene base to the GBAP1 gene base occurs between sites 5 and 6 for haplotype 2. For example, the number of reads with a GBA base at site 5 and site 6 is 98, the number of reads with a GBA base at site 5 and a GBAP1 base at site 6 is 105, and the number of reads with a GBAP1 base at site 5 and site 6 is 190. Given that the number of reads with a GBA base at sites 5 and 6 is 98 and the number of reads with a GBA base at site 5 and a GBAP1 base at site 6 is 105, the computing system can determine that the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype. The computing system can determine that the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype without using reads from the wild-type GBAP1 haplotype.

[0110] As another example, by analyzing linkage information between the GBA / GBAP1 discriminator bases, the following haplotypes can be determined for a subject (at six sites for illustrative purposes only):

[0111] [Table 7] A conversion from a GBA base to a GBAP1 base occurs between sites 5 and 6 of haplotype 2. For example, the number of reads with a GBA base at sites 5 and 6 is 98, the number of reads with a GBA base at site 5 and a GBAP1 base at site 6 is 105, the number of reads with a GBAP1 base at site 5 and a GBA base at site 6 is 95, and the number of reads with a GBAP1 base at sites 5 and 6 is 104. Given that the number of reads with a GBA base at sites 5 and 6 is 98 and the number of reads with a GBA base at site 5 and a GBAP1 base at site 6 is 105, a computing system can determine that the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype. A computing system can determine that the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype without using reads from the wild-type GBAP1 haplotype. The computing system can determine, without using reads with a GBAP1 base at site 5 and a GBA base at site 6, that the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype because that haplotype (haplotype 3) has most of the GBAP1 bases at the differentiating base sites / positions and is therefore less likely / more likely to be a GBA variant haplotype.

[0112] The copy number of the wild-type GBA haplotype can be 1 (i.e., carrier of the GBA variant haplotype). The computing system can determine the GBA status of the subject as a carrier of the GBA variant haplotype. The one or more haplotypes can include four haplotypes. The total copy number of the GBA gene and the GBAP1 gene can be four. The copy number of each of the four haplotypes can be one (e.g., one copy of the wild-type GBA haplotype, one copy of the GBA variant haplotype, one copy of the GBAP1 wild-type haplotype, and one copy of a haplotype having a high percentage (such as 80%, 85%, 90%, 95%, or more) of the discriminatory bases) and thus unlikely to be a GBA variant haplotype / probable to be a GBAP1 variant haplotype. The computing system can determine the GBA status of the subject as a carrier of the GBA variant haplotype. For example, by analyzing the linkage information between the GBA / GBAP1 discriminatory bases, the following haplotypes can be determined for the subject (at six sites for illustrative purposes only):

[0113] [Table 8] The subject is a carrier of a GBA mutant haplotype because the subject has one copy of a wild-type GBA haplotype and one copy of a GBA mutant haplotype (and one copy of a wild-type GBAP1 haplotype and one copy of a haplotype that has predominantly GBAP1 bases at the discriminatory base site / position and is therefore unlikely to be a GBA mutant haplotype / likely to be a GBAP1 mutant haplotype).

[0114] The one or more haplotypes may include three haplotypes. The total copy number of the GBA gene and the GBAP1 gene may be four. The copy numbers of the wild-type GBA haplotype, the GBA mutant haplotype, and the wild-type GBAP1 haplotype (or haplotypes having GBAP1 bases at a high percentage of discriminant bases, such as 80%, 85%, 90%, 95%, or more) may be one, one, and two, respectively. The computing system may determine the subject's GBA status as a carrier of a GBA mutant haplotype. For example, by analyzing the linkage information between the GBA / GBAP1 discriminant bases, the following haplotypes may be determined for the subject (at six sites for illustrative purposes only):

[0115] [Table 9] Because the subject has one copy of the wild-type GBA haplotype and one copy of the GBA mutant haplotype, the subject is a carrier of the GBA mutant haplotype.

[0116] Compound heterozygosity. The one or more haplotypes may include two or more GBA variant haplotypes. None of the two or more GBA variant haplotypes may include a GBA base at each of the plurality of GBA / GBAP1 discriminating bases. None of the two or more GBA variant haplotypes may include a GBA base at all of the plurality of GBA / GBAP1 discriminating bases. Each of the two or more GBA variant haplotypes may include one or more GBAP1 bases at one or more of the plurality of GBA / GBAP1 discriminating bases. The computing system may determine the subject's GBA status as a compound heterozygosity of GBA variant haplotypes. For example, the following haplotypes may be determined for the subject by analyzing linkage information between the GBA / GBAP1 discriminating bases (at six sites for illustrative purposes only):

[0117] [Table 10] Because the subject does not have any copies of the wild-type GBA haplotype, and has one copy of each of the two GBA mutant haplotypes, the subject is a compound heterozygote for the GBA mutant haplotypes.

[0118] Homozygous. One or more haplotypes can include an identical base (e.g., a GBA base or a GBAP1 base) at the GBA / GBAP1 discriminating base, or at each of two or more of the plurality of GBA / GBAP1 discriminating bases. The computing system can determine that the subject is homozygous at one or more of the plurality of GBA / GBAP1 discriminating bases (e.g., homozygous for a wild-type GBAP1 gene haplotype or homozygous for a GBA mutant haplotype).

[0119] The computing system can use sequence reads of the second plurality of sequence reads, each including a base in the GBA / GBAP1 discriminating base that is not a GBA base, to determine that the copy number of the GBA base in each of one or more of the plurality of GBA / GBAP1 discriminating bases is 0. The base in the GBA / GBAP1 discriminating base that is not a GBA base can be a GBAP1 base. The computing system can determine that the subject's GBA status is homozygous for the GBA variant haplotype in one, one or more, or each of the plurality of GBA / GBAP1 discriminating bases. For example, based on the plurality of GBA / GBAP1 discriminating bases, the computing system can determine the copy number (CN) of the GBA base. The number of reads supporting the GBA base or the GBAP1 base, and the total CN of the GBA gene and the GBAP1 gene can be used to determine the most likely combination of the CN of the GBA base and the CN of the GBAP1 base. If the CN of the GBA base is determined to be 0, this indicates that the subject has no copies of the wild-type GBA gene haplotype (a haplotype that carries the GBA base at the mutation site of interest) and is homozygous for the GBA gene mutant haplotype.

[0120] The method 400 ends at block 428.

[0121] Determination of CYP21A2 variants and variant status FIG. 5 is a flow diagram illustrating an exemplary method 500 for determining or identifying one or more CYP21A2 variants or variant status. Method 500 may be embodied in a set of executable program instructions stored on a computer-readable medium, such as one or more disk drives of a computing system. For example, a computing system 700, shown in FIG. 7 and described in more detail below, may execute a set of executable program instructions to perform method 500. When method 500 is initiated, the executable program instructions may be loaded into a memory, such as a RAM, and executed by one or more processors of computing system 700. Although method 500 is described with respect to computing system 700 shown in FIG. 7, the description is merely exemplary and not intended to be limiting. In some embodiments, method 500 or portions thereof may be performed serially or in parallel by multiple computing systems.

[0122] After the method 500 starts at block 504, the method 500 proceeds to block 508, where a computing system (e.g., computing system 700 described with reference to FIG. 7) aligns the first plurality of sequence reads to a reference sequence (e.g., a reference genome sequence such as hg19 or hg38) to obtain a second plurality of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference sequence (including an alignment of each of the second plurality of sequence reads to the CYP21A2 gene or the CYP21A1P pseudogene in the reference sequence). The computing system can receive the first plurality of sequence reads generated from a sample obtained from the subject. The computing system can store the first plurality of sequence reads in a memory. The computing system can load the first plurality of sequence reads into the memory. The sequence reads can be generated by a technique such as sequencing by synthesis, sequencing by ligation, or sequencing by ligation. Sequence reads can be generated using instruments such as the MINISEQ, MISEQ, NEXTSEQ, HISEQ, and NOVASEQ sequencing instruments from Illumina, Inc. (San Diego, Calif.).

[0123] The sequence reads can be, for example, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1250, 1500, 1750, 2000 or more base pairs (bps) in length, respectively. For example, the sequence reads are about 100 base pairs to about 1000 base pairs in length, respectively. The sequence reads can include paired-end sequence reads. The sequence reads can include single-end sequence reads. The sequence reads can be generated by whole genome sequencing (WGS). The WGS can be clinical WGS (cWGS). The sample can include cells, cell-free DNA, cell-free fetal DNA, amniotic fluid, a blood sample, a biopsy sample, or a combination thereof.

[0124] The sequence reads can be aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference sequence with an alignment quality score of at least 0. The sequence reads can be aligned to the CYP21A2 gene or the CYP21A1P gene in the reference sequence with an alignment quality score of about 0 (e.g., when the sequences are aligned to regions where genes and gene paralogs are highly homologous). The computing systems used were Burrows-Wheeler Aligner (BWA), iSAAC, BarraCUDA, BFAST, BLASTN, BLAT, Bowtie, CASHX, Cloudburst, CUDA-EC, CUSHAW, CUSHAW2, CUSHAW2-GPU, drFAST, ELAND, ERNE, GNUMAP, GEM, GensearchNGS, GMAP and GSNAP, Geneious Assembler, LAST, MAQ, mrFAST and mrsFAST, MOM, MOSAIK, MPscan, Novoaligh and NovoalignCS, NextGENe, Omixon, PALMapper, Partek, PASS, PerM, PRIMEX, QPalma, RazerS, REAL, cREAL, RMAP, rNA, RT Sequence reads can be aligned to a reference sequence using an aligner or alignment method such as Investigator, Segemehl, SeqMap, Shrec, SHRiMP, SLIDER, SOAP, SOAP2, SOAP3 and SOAP3-dp, SOCS, SSAHA and SSAHA2, Stampy, SToRM, Subread and Subjunc, Taipan, UGENE, VelociMapper, XpressAlign and ZOOM.

[0125] Method 500 proceeds from block 508 to block 512, where the computing system determines the number (e.g., normalized and / or corrected number) of sequence reads of the second plurality of sequence reads that are aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference sequence.

[0126] The computing system can determine a normalized number of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference sequence. The computing system can determine the normalized number of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference sequence using (1a) the depth of the sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene, (1b) the length of the unique region, (2a) the depth of the sequence reads of a first plurality of sequence reads aligned to each of a plurality of regions of the reference sequence other than the locus including the CYP21A2 gene and the CYP21A1P pseudogene, and / or (2b) the length of each of a plurality of regions of the reference other than the locus including the CYP21A2 gene and the CYP21A1P pseudogene.

[0127] The computing system can determine a normalized, corrected number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference sequence from the normalized number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference sequence. The computing system can determine a normalized, GC content corrected number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference sequence from the normalized number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference sequence to determine the normalized, corrected number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference sequence. The computing system can determine a normalized GC content corrected number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference sequence from the normalized number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference sequence using (1) the GC content of the CYP21A2 gene or CYP21A1P pseudogene. The computing system can use (1) the GC content of the CYP21A2 gene or CYP21A1P pseudogene, and / or (2) the respective GC contents of one or more regions of the reference sequence other than the locus including the CYP21A2 gene and the CYP21A1P pseudogene (or one or more regions of the reference sequence not including the CYP21A2 gene and the CYP21A1P pseudogene) to determine a normalized GC content-corrected number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference sequence from the normalized number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference sequence.For example, the computing system can determine a normalized GC content corrected number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference sequence from the normalized number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene in the reference sequence using (1) the GC content of the CYP21A2 gene or CYP21A1P pseudogene, and (2) the GC content of a region of the reference sequence other than the locus including the CYP21A2 gene and the CYP21A1P pseudogene. As another example, the computing system can use the GC content of multiple regions (e.g., 2, 3, 4, 5, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500, 1000, 2000, 3000, 4000, 5000, 10000, or more regions) of the reference sequence other than the locus including (1) the CYP21A2 gene and the CYP21A1P pseudogene, and (2) the CYP21A2 gene and the CYP21A1P pseudogene to determine the number of sequence reads corrected for the normalized GC content of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference sequence from the normalized number of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference sequence.

[0128] Method 500 proceeds from block 512 to block 516, where the computing system determines the total copy number of the CYP21A2 gene or CYP21A1P pseudogene using a Gaussian mixture model including multiple Gaussians, each representing a different integer copy number, given the normalized number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene. The computing system can determine the total copy number of the CYP21A2 gene and the CYP21A1P pseudogene using a Gaussian mixture model given the normalized number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene. The computing system can determine the total copy number of the CYP21A2 gene and the CYP21A1P pseudogene using a Gaussian mixture model given the normalized, corrected number of sequence reads aligned to the CYP21A2 gene or CYP21A1P pseudogene.

[0129] The total copy number can be, for example, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more. The Gaussian mixture model can include a one-dimensional Gaussian mixture model. The Gaussians of the Gaussian mixture model can represent integer copy numbers, for example, 0 to 5, 0 to 6, 0 to 7, 0 to 8, 0 to 9, 0 to 10, 0 to 11, 0 to 12, 0 to 13, 0 to 14 or 0 to 15. For example, the Gaussians of the Gaussian mixture model can represent integer copy numbers 0 to 10. The average of each of the multiple Gaussians can be the integer copy number represented by the Gaussian. The average of each of the multiple Gaussians can be the integer copy number represented by the Gaussian (for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 or more copies). The standard deviation of the Gaussians can be, for example, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 or more, or can be about 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 or more. The multiple Gaussians of the Gaussian mixture model can include, for example, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, or more Gaussians. For example, the multiple Gaussians of the Gaussian mixture model can include 5 Gaussians.

[0130] Given the normalized number of sequence reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene, a computing system can determine the total copy number of the CYP21A2 gene and the CYP21A1P pseudogene using a Gaussian mixture model and a predefined posterior probability threshold. The predetermined posterior probability threshold may be, for example, 0.70, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.80, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.90, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, or more, 0.7, 0.75, 0.8, 0.85, 0.95, or more. or can be about 0.70, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.80, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.90, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, or more, 0.7, 0.75, 0.8, 0.85, 0.95, or more. For example, the predetermined posterior probability threshold is 0.95.

[0131] Method 500 proceeds from block 516 to block 520, where the computing system uses sequence reads of a second plurality of sequence reads aligned to a region including or corresponding to a plurality of CYP21A2 / CYP21A1P discriminant bases to phase one or more haplotypes derived from the CYP21A2 gene and the CYP21A1P pseudogene in a region of the CYP21A2 gene or a corresponding region of the CYP21A1P pseudogene that includes a plurality of CYP21A2 / CYP21A1P discriminant bases (or positions or sites of the discriminant bases). For example, the sequence reads can be aligned to the reference sequence such that the sequence reads overlap with the CYP21A2 / CYP21A1P discriminatory bases (or the sites of the CYP21A2 / CYP21A1P discriminatory bases) or such that the bases of the sequence reads are aligned to the CYP21A2 / CYP21A1P paralog discriminatory bases (or the sites of the CYP21A2 / CYP21A1P paralog discriminatory bases). The sequence reads of the second plurality of sequence reads are aligned to a region of the CYP21A2 gene or a corresponding region of the CYP21A1P pseudogene that includes the plurality of CYP21A2 / CYP21A1P discriminatory bases with an alignment quality score of 0 or greater. The haplotypes can include reciprocal recombination variants. The haplotypes can include non-reciprocal recombination variants or gene conversion variants.

[0132] The one or more haplotypes may include a wild-type CYP21A2 haplotype, a wild-type CYP21A1P, and / or a CYP21A2 / CYP21A1P hybrid haplotype. The CYP21A2 / CYP21A1P hybrid haplotype may include both a CYP21A2 base and a CYP21A1P base. The CYP21A2 / CYP21A1P hybrid haplotype may be a recombinant variant. The CYP21A2 / CYP21A1P hybrid haplotype may include a CYP21A2 variant haplotype or a CYP21A1P variant haplotype.

[0133] To determine the phase of one or more haplotypes derived from the CYP21A2 gene or the CYP21A1P pseudogene, the computing system can analyze linkage information between the CYP21A2 / CYP21A1P discriminating bases of the plurality of CYP21A2 / CYP21A1P discriminating bases using sequence reads of the second plurality of sequence reads aligned to a region including or corresponding to the plurality of CYP21A2 / CYP21A1P discriminating bases. The computing system can phase one or more haplotypes derived from the CYP21A2 gene or the CYP21A1P pseudogene using sequence reads of the second plurality of sequence reads aligned to two or more of the plurality of CYP21A2 / CYP21A1P discriminating bases, respectively. For example, with reference to Figure 1, assuming that gene A and gene B shown in the figure are the CYP21A2 gene or the CYP21A1P pseudogene, respectively, one read pair covering site 1 and site 4 may indicate the haplotype (haplotype x) from which the read pair is derived having a CYP21A2 gene base at site 1 and a CYP21A2 gene base at site 4. One read pair covering site 3 and site 5 may indicate the haplotype (haplotype y) from which the read pair is derived has a CYP21A2 gene base at site 3 and a CYP21A1P pseudogene base at site 5. One read covering site 4 and site 5 may indicate the haplotype (haplotype y) from which the read pair is derived has a CYP21A2 gene base at site 4 and a CYP21A1P pseudogene base at site 5. The computer system can phase all haplotypes derived from either the CYP21A2 gene or the CYP21A1P pseudogene in regions with a set of five reliable base differences and identify haplotypes of the CYP21A2 gene (haplotype 1) and haplotypes of the CYP21A1P pseudogene (haplotype 2) as well as hybrid haplotypes (haplotypes 3 and 4). The number of haplotypes, haplotype bases at sites, number of sites, and site (e.g., sites 1 and 4, or sites 3 and 5) sequence read coverage are shown in FIG. 1 for illustrative purposes only and are not intended to be limiting.

[0134] The plurality of CYP21A2 / CYP21A1P discriminant bases can include 14 (or 11, 12, 13, 15, 16, 17, or more) CYP21A2 / CYP21A1P discriminant bases. The 14 CYP21A2 / CYP21A1P discriminant bases can include 9 (or 6, 7, 8, 10, 11, 12, or more) CYP21A2 / CYP21A1P recombinant variants. The CYP21A2 / CYP21A1P identification bases are chr6 of hg38: 32039081 / 32006353, 32039128 / 32006400, 32039132 / 32006404, 32039143 / 32006407, 32039426 / 32006690, 32039548 / 32006812, 32039802 / 320070 66, 32039807 / 32007071, 32039810 / 32007074, 32039816 / 32007080, 32040182 / 32007446, 32040216 / 32007481, 32040421 / 32007686, and 32040535 / 32007800, or these bases of a reference human genome sequence.

[0135] Method 500 proceeds from block 520 to block 524, where the computing system determines a copy number for each of the one or more haplotypes using the total copy number of the CYP21A2 gene and the CYP21A1P pseudogene and the number of sequence reads of the second plurality of sequence reads that each include one or more of the plurality of CYP21A2 / CYP21A1P discriminatory bases that support the haplotype. The copy number of the haplotype can be, for example, 1, 2, 3, 4 or more.

[0136] The computing system can determine the subject's CYP21A2 status (e.g., carrier, compound heterozygous, or homozygous) using one or more haplotypes from the CYP21A2 gene or the CYP21A1P pseudogene and / or the copy number of each of the one or more haplotypes in a region of the CYP21A2 gene or a corresponding region of the CYP21A1P pseudogene. The computing system can generate a user interface (UI), such as a graphical user interface, that includes a UI element that represents or includes the CYP21A2 status. The UI can include the CYP21A2 status as part of the UI element. The UI element can be a window (e.g., a container window, a browser window, a text terminal, a child window, or a message window), a menu (e.g., a menu bar, a context menu, or a menu extra), an icon, or a tab. The UI element can be for an input control (e.g., a check box, a radio button, a drop-down list, a list box, a button, a toggle, a text field, or a date field). A UI element can be navigational (e.g., breadcrumbs, sliders, search fields, pagination, sliders, tags, icons). A UI element can provide information (e.g., tooltips, icons, progress bars, notifications, message boxes, or modal windows). A UI element can be a container (e.g., an accordion).

[0137] Carrier. To determine the copy number of each of the one or more haplotypes, the computing system can determine that the likelihood of one copy of the wild-type CYP21A2 haplotype is higher than the likelihood of two copies of the wild-type CYP21A2 haplotype given the number of sequence reads of the second plurality of sequence reads, each of which includes one or more of the plurality of CYP21A2 / CYP21A1P discriminatory bases supporting the wild-type CYP21A2 haplotype. The computing system can determine (e.g., using the sequence reads) for one, one or more (e.g., two, three, or four), or each of one or more haplotypes that the likelihood of one copy of the wild-type CYP21A2 haplotype is higher than the likelihood of two copies of the wild-type CYP21A2 haplotype. The likelihood difference can be, for example, 1%, 2%, 3%, 5%, 10%, 15%, 20%, or more. The computing system can determine that the copy number of the wild-type CYP21A2 haplotype is 1. The computing system can determine that the likelihood of one copy of the wild-type CYP21A2 haplotype is higher than the likelihood of two copies of the wild-type CYP21A2 haplotype for each of one or more pairs (or all pairs) of a plurality of consecutive CYP21A2 / CYP21A1P discriminatory bases, where a first haplotype of the one or more haplotypes includes a CYP21A2 base in the consecutive CYP21A2 / CYP21A1P discriminatory bases and a second haplotype of the one or more haplotypes includes a CYP21A2 base and a CYP21A1P base (or a CYP21A1P base and a CYP21A2 base) in the consecutive CYP21A2 / CYP21A1P discriminatory bases. The first haplotype can include a CYP21A2 base in consecutive CYP21A2 / CYP21A1P discriminator bases. The second haplotype can include a transition from a CYP21A2 base to a CYP21A1P base (or a CYP21A1P base to a CYP21A2 base) between the consecutive CYP21A2 / CYP21A1P discriminator bases.Consecutive CYP21A2 / CYP21A1P discriminatory bases are contiguous within the plurality of CYP21A2 / CYP21A1P discriminatory bases, regardless of whether the contiguous CYP21A2 / CYP21A1P discriminatory bases are adjacent bases in the reference sequence. The computing system can combine (e.g., average or weighted average) the likelihoods determined for each of one or more pairs (or all pairs) of contiguous CYP21A2 / CYP21A1P discriminatory bases to determine that the likelihood of one copy of a wild-type CYP21A2 haplotype is higher than the likelihood of two copies of a wild-type CYP21A2 haplotype.

[0138] The computing system can (1) determine that, given a number of sequence reads of a second plurality of sequence reads, each of which includes a CYP21A2 base in consecutive CYP21A2 / CYP21A1P discriminant bases, for each of one or more pairs (or all pairs) of consecutive CYP21A2 / CYP21A1P discriminant bases, the likelihood of one copy of the wild-type CYP21A2 haplotype is higher than the likelihood of two copies of the wild-type CYP21A2 haplotype. The computing system can determine that, given (2) the number of sequence reads in the second plurality of sequence reads that include a CYP21A2 base and a CYP21A1P base, respectively, in consecutive CYP21A2 / CYP21A1P discriminant bases, and / or (3) the number of sequence reads in the second plurality of sequence reads that include a CYP21A1P base and a CYP21A2 base, respectively, in consecutive CYP21A2 / CYP21A1P discriminant bases, the likelihood of one copy of the wild-type CYP21A2 haplotype is higher than the likelihood of two copies of the wild-type CYP21A2 haplotype in each of one or more pairs (or all pairs) of consecutive CYP21A2 / CYP21A1P discriminant bases. The computing system can (4) determine that, given a number of sequence reads of a second plurality of sequence reads, each of which includes a CYP21A1P base in consecutive CYP21A2 / CYP21A1P distinguishing bases, for each of one or more pairs (or all pairs) of consecutive CYP21A2 / CYP21A1P distinguishing bases, the likelihood of one copy of the wild-type CYP21A2 haplotype is higher than the likelihood of two copies of the wild-type CYP21A2 haplotype.

[0139] The copy number of the wild-type CYP21A2 haplotype can be 1. The computing system can determine that the subject is a carrier of a CYP21A2 variant haplotype. The one or more haplotypes can include four haplotypes. The total copy number of the CYP21A2 and CYP21A1P genes can be four. The copy number of each of the four haplotypes can be one (e.g., one copy of the wild-type CYP21A2 haplotype, one copy of the CYP21A2 variant haplotype, one copy of the CYP21A1P wild-type haplotype, and one copy of a haplotype having a high percentage (e.g., 80%, 85%, 90%, 95%, or more) of CYP21A1P bases at the discriminating bases), and thus unlikely to be a CYP21A2 variant haplotype / likely to be a CYP21A1P variant haplotype. The computing system can determine the subject's CYP21A2 status as a carrier of a CYP21A2 variant haplotype. The one or more haplotypes can include three haplotypes. The total copy number of the CYP21A2 and CYP21A1P genes can be four. The copy numbers of the wild-type CYP21A2 haplotype, the CYP21A2 variant haplotype, and the CYP21A1P wild-type haplotype (or a haplotype having CYP21A1P bases at a high percentage, such as 80%, 85%, 90%, 95%, or more, of the discriminating bases) can be one, one, and two, respectively. The computing system can determine that the subject is a carrier of a CYP21A2 variant haplotype.

[0140] Compound heterozygotes. The one or more haplotypes may include two or more haplotypes. None of the two or more haplotypes may include a CYP21A2 base at each of the plurality of CYP21A2 / CYP21A1P discriminatory bases. None of the two or more haplotypes may include a CYP21A2 base at all of the plurality of CYP21A2 / CYP21A1P discriminatory bases. Each of the two or more haplotypes may include a CYP21A1P base at one or more of the plurality of CYP21A2 / CYP21A1P discriminatory bases. The computing system may determine that the subject is a compound heterozygote for a CYP21A2 variant haplotype.

[0141] Homozygous. One or more haplotypes can include identical bases (e.g., CYP21A bases or CYP21A1P bases) at the CYP21A2 / CYP21A1P discriminatory base, or at each of two or more of the plurality of CYP21A2 / CYP21A1P discriminatory bases. The computing system can determine that the subject is homozygous at one or more of the plurality of CYP21A2 / CYP21A1P discriminatory bases (e.g., homozygous for a wild-type CYP21A1P gene haplotype or homozygous for a CYP21A2 gene variant).

[0142] The one or more haplotypes may include only one haplotype. The only haplotype may not include a CYP21A2 base at one, one or more, or each of the multiple CYP21A2 / CYP21A1P discriminatory bases. The computing system can determine that the subject is homozygous for a CYP21A2 variant haplotype. For example, based on the multiple CYP21A2 / CYP21A1P discriminatory bases, the computing system can determine the copy number (CN) of the CYP21A2 base. The number of reads supporting the CYP21A2 base or the CYP21A1P gene base, and the total CN of the CYP21A2 gene and the CYP21A1P gene can be used to determine the most likely combination of the CN of the CYP21A2 base and the CN of the CYP21A1P gene base. If the CN of a CYP21A2 base is determined to be 0, this indicates that the subject has no copies of the wild-type CYP21A2 gene haplotype (a haplotype that has the CYP21A2 gene base at the discriminatory base of interest) and is homozygous for the CYP21A2 gene variant haplotype.

[0143] The method 500 ends at block 528.

[0144] Determination of transgenic variants and genetic variant status FIG. 6 is a flow diagram illustrating an exemplary method 600 for determining or identifying one or more genetic variants or genetic variant status. The method 600 may be embodied in a set of executable program instructions stored on a computer-readable medium, such as one or more disk drives of a computing system. For example, a computing system 700, shown in FIG. 7 and described in more detail below, may execute a set of executable program instructions to perform the method 600. When the method 600 is initiated, the executable program instructions may be loaded into a memory, such as a RAM, and executed by one or more processors of the computing system 700. Although the method 600 is described with respect to the computing system 700 shown in FIG. 7, the description is merely exemplary and is not intended to be limiting. In some embodiments, the method 600 or portions thereof may be executed by multiple computing systems, either serially or in parallel.

[0145] After the method 600 starts at block 604, the method 600 proceeds to block 608, where a computing system (such as the computing system 700 described with reference to FIG. 7) aligns the first plurality of sequence reads to a reference sequence (e.g., a reference genome sequence such as hg19 or hg38) to obtain a second plurality of sequence reads aligned to genes or gene paralogs (or regions therebetween) in the reference sequence (including alignment of each of the second plurality of sequence reads to genes or gene paralogs in the reference sequence). The computing system can receive the first plurality of sequence reads generated from a sample obtained from the subject. The computing system can store the first plurality of sequence reads in a memory. The computing system can load the first plurality of sequence reads into the memory. The sequence reads can be generated by a technique such as sequencing-by-synthesis, sequencing-by-ligation, or sequencing-by-ligation. Sequence reads can be generated using instruments such as the MINISEQ, MISEQ, NEXTSEQ, HISEQ, and NOVASEQ sequencing instruments from Illumina, Inc. (San Diego, Calif.).

[0146] The sequence reads can be, for example, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1250, 1500, 1750, 2000 or more base pairs (bps) in length, respectively. For example, the sequence reads are about 100 base pairs to about 1000 base pairs in length, respectively. The sequence reads can include paired-end sequence reads. The sequence reads can include single-end sequence reads. The sequence reads can be generated by whole genome sequencing (WGS). The WGS can be clinical WGS (cWGS). The sample can include cells, cell-free DNA, cell-free fetal DNA, amniotic fluid, a blood sample, a biopsy sample, or a combination thereof.

[0147] The sequence reads may be aligned to genes or pseudogenes in the reference sequence with an alignment quality score of greater than or equal to 0. The sequence reads may be aligned to genes or pseudogenes in the reference sequence with an alignment quality score of about 0 (e.g., when sequences are aligned to regions where genes and gene paralogs are highly homologous). The computing systems used were Burrows-Wheeler Aligner (BWA), iSAAC, BarraCUDA, BFAST, BLASTN, BLAT, Bowtie, CASHX, Cloudburst, CUDA-EC, CUSHAW, CUSHAW2, CUSHAW2-GPU, drFAST, ELAND, ERNE, GNUMAP, GEM, GensearchNGS, GMAP and GSNAP, Geneious Assembler, LAST, MAQ, mrFAST and mrsFAST, MOM, MOSAIK, MPscan, Novoaligh and NovoalignCS, NextGENe, Omixon, PALMapper, Partek, PASS, PerM, PRIMEX, QPalma, RazerS, REAL, cREAL, RMAP, rNA, RT Sequence reads can be aligned to a reference sequence using an aligner or alignment method such as Investigator, Segemehl, SeqMap, Shrec, SHRiMP, SLIDER, SOAP, SOAP2, SOAP3 and SOAP3-dp, SOCS, SSAHA and SSAHA2, Stampy, SToRM, Subread and Subjunc, Taipan, UGENE, VelociMapper, XpressAlign and ZOOM.

[0148] The gene paralog can be a gene. The gene paralog can be a pseudogene. The gene and the gene paralog have at least 70%, 75%, 80%, 85%, 90%, 95%, 96%, 97%, 98%, 99% or more sequence identity. In some embodiments, the gene is the GBA gene and the gene paralog is the GBAP1 gene. If the gene is the GBA gene and the gene paralog is the GBAP1 gene, the computing system can perform the method 400 (or one or more things of the method 400) described with reference to FIG. 4. In some embodiments, the gene is the CYP21A2 gene and the gene paralog is the CYP21A1P pseudogene. If the gene is the CYP21A2 gene and the gene paralog is the CYP21A1P pseudogene, the computing system can perform the method 500 (or one or more things of the method 500) described with reference to FIG. 5. In some embodiments, the genes are ABCC6, ABCD1, ACTB, ACTG1, ACTN4, ADAMTSL2, ADIPOR1, AFG3L2, AGK, ALG1, ALMS1, ANKRD11, ANOS1, AP4S1, ARMC4, ARSE, ASNS, ATAD3A, B3GAT3, BCAP31, BDP1, BMPR1A, BRAF, BRCA1, C2, CACNA1C, CALM1, CD46, CEP290, CFH, CFH, CFH, CHEK2, CISD2, CLCNKA, CLCNKB, CORO1A, COX10, CP, CRYBB2, CSF2RA, CUBN, CUBN, CYCS, CYP11B1, CYP2 1A2, DCLRE1C, DHFR, DICER1, DIS3L2, DNAH11, DNAH11, DNM1, DSE, DUOX2, EGLN1, ELK1, ELMO2, ERCC6, ESPN, EYS, F8, FANCD2, FANCD2, FAR1, FHL1, FLG, FLNC, FOXD4, FXN, GBA, GH 1, GJA1, GK, GLUD1, GLUD1, GOSR2, GUSB, HBA1, HBA2, HNRNPA1, HPS1, HSPD1, HYDIN, IDS, IFT122, IGLL1, KANSL1, KCTD1, KIF1C, KRAS, KRT14, KRT16, KRT17, KRT6A, KRT6B, KRT6C,LEFTY2, LRP5, LRP5, MAT2A, MID1, MOCS1, MSN, MSX2, MYO5B, NCF1, NEB, NECAP1, NEFH, NF1, NF1, NF1, NOTCH2, NXF5, OCLN, OTOA, PARN, PBX1, PIGA, PIGN, PIK3CA, PIK3CD, PKD1, PKP2, PMS2, PMS2, PMS2, PNPT1, POLH, PRODH, PRODH, PROS1, PRPS1, PRSS1, PTEN, RAD21, RBM8A, RBPJ, RDX, RMND1, RNF216, RNF216, RPL15, SALL1, SBDS, SDHA, SHOX, SLC25A15, SLC25A15, SLC33A1, SLC6A8, SMN1, SMN2, SOX2, SPTLC1, SRD5A3, SRP72, STAT5B, STRC, SYT14, TARDBP, TBL1XR1, TBX20, TIMM8A, TP M3, TPMT, TRAPPC2, TRIP11, TTN, TUBA1A, TUBB2A, TUBB2B, TUBB3, TUBB4A, TUBG1, TYR, UBA5, UBE3A, UNC93B1, USP18, VPS35, VWF, WRN, XIAP, ZEB2, or ZNF341. ,

[0149] In some embodiments, the computing system can determine the number of sequence reads aligned to a gene or gene paralog (or a region therebetween). The number of sequence reads aligned to a gene or gene paralog, or a region therebetween, includes a normalized and / or GC-corrected number of sequence reads aligned to a gene or gene paralog, or a region therebetween.

[0150] The computing system can determine a normalized number of sequence reads aligned to genes or gene paralogs in the reference sequence using (1a) the depth of sequence reads aligned to genes or gene paralogs, (1b) the length of the unique region, (2a) the depth of sequence reads of the first plurality of sequence reads aligned to each of the plurality of regions of the reference sequence other than the locus containing the gene and gene paralog, and / or (2b) the length of each of the plurality of regions of the reference other than the locus containing the gene and gene paralog. The computing system can determine a normalized corrected number of sequence reads aligned to genes or gene paralogs in the reference sequence from the normalized number of sequence reads aligned to genes or gene paralogs in the reference sequence. To determine the normalized corrected number of sequence reads aligned to genes or gene paralogs in the reference sequence, the computing system can include determining a normalized GC content corrected number of sequence reads aligned to genes or gene paralogs in the reference sequence from the normalized number of sequence reads aligned to genes or gene paralogs in the reference sequence. The computing system can determine a normalized GC content-corrected number of sequence reads aligned to genes or gene paralogs in the reference sequence from the normalized number of sequence reads aligned to genes or gene paralogs in the reference sequence using (1) the GC content of the genes or gene paralogs, and / or (2) the GC content of each of one or more regions of the reference sequence other than the locus that includes the genes and gene paralogs (or one or more regions of the reference sequence that do not include the genes and gene paralogs).For example, the computing system can use (1) the GC content of the gene or gene paralog, and (2) the GC content of a region of the reference sequence other than the locus that contains the gene and gene paralog to determine a normalized GC content-corrected number of sequence reads aligned to the gene or gene paralog in the reference sequence from the normalized number of sequence reads aligned to the gene or gene paralog in the reference sequence. As another example, the computing system can determine a normalized GC content corrected number of sequence reads aligned to genes or gene paralogs in the reference sequence from a normalized number of sequence reads aligned to genes or gene paralogs in the reference sequence using (1) the GC content of the gene or gene paralog, and (2) the GC content of multiple regions (e.g., 2, 3, 4, 5, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500, 1000, 2000, 3000, 4000, 5000, 10000, or more regions) of the reference sequence other than the locus that includes the gene and gene paralog.

[0151] From block 608, the method 600 proceeds to block 612, where the computing system determines the total copy number of genes and gene paralogs using a Gaussian mixture model that includes multiple Gaussians, each representing a different integer copy number, given the number of sequence reads aligned to a gene or gene paralog. The computing system can determine the total copy number of genes and gene paralogs using a Gaussian mixture model, given the number of sequence reads aligned to a gene or gene paralog (e.g., the number of normalized and / or corrected sequence reads).

[0152] The total copy number can be, for example, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more. The Gaussian mixture model can include a one-dimensional Gaussian mixture model. The Gaussians of the Gaussian mixture model can represent integer copy numbers, for example, 0 to 5, 0 to 6, 0 to 7, 0 to 8, 0 to 9, 0 to 10, 0 to 11, 0 to 12, 0 to 13, 0 to 14 or 0 to 15. For example, the Gaussians of the Gaussian mixture model can represent integer copy numbers 0 to 10. The average of each of the multiple Gaussians can be the integer copy number represented by the Gaussian. The average of each of the multiple Gaussians can be the integer copy number represented by the Gaussian (for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 or more copies). The standard deviation of the Gaussians can be, for example, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 or more, or can be about 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 or more. The multiple Gaussians of the Gaussian mixture model can include, for example, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, or more Gaussians. For example, the multiple Gaussians of the Gaussian mixture model can include 5 Gaussians.

[0153] To determine the total copy number of genes and gene paralogues, a computing system can use a Gaussian mixture model to determine the copy number of the region between genes or gene paralogues given the normalized number of sequence reads aligned to the genes or gene paralogues. The total copy number of genes and gene paralogues can be the copy number of the region between genes or gene paralogues plus 2.

[0154] The computing system can determine the total copy number of genes and gene paralogs using a Gaussian mixture model and a predefined posterior probability threshold given the normalized number of sequence reads aligned to a gene or gene paralog, the predefined posterior probability threshold being greater than or equal to 0.70, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.80, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.90, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99. The predetermined posterior probability threshold can be, or can be about 0.70, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.80, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.90, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99 or greater. For example, the predetermined posterior probability threshold is 0.95.

[0155] Method 600 proceeds from block 612 to block 616, where the computing system phases one or more haplotypes of or derived from a gene (including recombinant variants of a gene) or gene paralog comprising a plurality of gene / gene paralog discriminating bases (or positions or sites of discriminating bases), or a region of a gene or a corresponding region of a gene paralog, using sequence reads of the second plurality of sequence reads aligned to the region comprising or corresponding to the plurality of gene / gene paralog discriminating bases. For example, the sequence reads can be aligned to the reference sequence such that the sequence reads overlap with the gene / gene paralog discriminating bases (or sites of the gene / gene paralog discriminating bases) or such that the bases of the sequence reads are aligned to the gene / gene paralog discriminating bases (or sites of the gene / gene paralog discriminating bases). The sequence reads can be aligned to the reference sequence of the second plurality of sequence reads and aligned to the region of a gene or a corresponding region of a gene paralog comprising a plurality of gene / gene paralog discriminating bases having an alignment quality score of 0 or greater.

[0156] The one or more haplotypes may include wild type gene haplotypes, wild type gene paralogs, and / or gene / gene paralog hybrid haplotypes. Gene / gene paralog hybrid haplotypes may include both gene bases and gene paralog bases. Gene / gene paralog hybrid haplotypes may be recombination variants. Gene / gene paralog hybrid haplotypes may include gene variant haplotypes or gene paralog variant haplotypes. Recombination variants may include reciprocal recombination variants. Recombination variants may include non-reciprocal recombination variants or gene conversion variants.

[0157] To phase one or more haplotypes derived from a gene or gene paralog, the computing system can use sequence reads of the second plurality of sequence reads aligned to a region including or corresponding to the plurality of gene / gene paralog discriminating bases to analyze linkage information between gene / gene paralog discriminating bases of the plurality of gene / gene paralog discriminating bases. The computing system can use sequence reads of the second plurality of sequence reads aligned to two or more of the plurality of gene / gene paralog discriminating bases, respectively, to phase one or more haplotypes derived from a gene or gene paralog. For example, referring to FIG. 1, assuming that gene A and gene B shown in the figure are genes or gene paralogs, respectively, one read pair covering site 1 and site 4 can indicate that the haplotype (haplotype x) from which the read pair is derived has a gene base at site 1 and a gene base at site 4. One read pair covering site 3 and site 5 can indicate that the haplotype (haplotype y) from which the read pair is derived has a gene base at site 3 and a gene paralog base at site 5. One read covering sites 4 and 5 can indicate that the haplotype from which the read originates (haplotype y) has a gene base at site 4 and a gene paralog base at site 5. The caller can phase all haplotypes that originate from either the gene or gene paralog in the region with 5 reliable base differences, identifying gene (haplotype 1) and gene paralog (haplotype 2) haplotypes as well as hybrid haplotypes (haplotypes 3 and 4). The number of haplotypes, haplotype bases at sites, number of sites, and site (e.g., sites 1 and 4, or sites 3 and 5) sequence read coverage are shown in FIG. 1 for illustrative purposes only and are not intended to be limiting.

[0158] Method 600 proceeds from block 616 to block 620, where the computing system determines a copy number of each of the one or more haplotypes using the total copy number of the gene and gene paralogs and the number of sequence reads of the second plurality of sequence reads each including one or more of the plurality of gene / gene paralog discriminating bases supporting the haplotype. The copy number of the haplotype may be, for example, one, two, three, four or more. The computing system may use the copy number of each of the one or more haplotypes derived from the gene or gene paralog, or a region of the gene or a corresponding region of the gene paralog, and / or the one or more haplotypes to determine the genetic variant status of the subject (e.g., carrier, compound heterozygous, or homozygous). The computing system may generate a user interface (UI), such as a graphical user interface, that includes UI elements that represent or include the genetic variant status. The UI may include the genetic variant status as part of the UI element. The UI element may be a window (e.g., a container window, a browser window, a text terminal, a child window, or a message window), a menu (e.g., a menu bar, a context menu, or a menu extra), an icon, or a tab. The UI element can be for input control (e.g., a checkbox, radio button, drop-down list, list box, button, toggle, text field, or date field). The UI element can be navigation (e.g., a breadcrumb, slider, search field, pagination, slider, tag, icon). The UI element can provide information (e.g., a tooltip, icon, progress bar, notification, message box, or modal window). The UI element can be a container (e.g., an accordion).

[0159] Carrier. To determine the copy number of each of the one or more haplotypes, the computing system can determine that the likelihood of one copy of the wild-type gene haplotype is higher than the likelihood of two copies of the wild-type gene haplotype given the number of sequence reads of the second plurality of sequence reads, each of which includes one or more of the plurality of gene / gene paralog discriminating bases supporting the wild-type gene haplotype. The computing system can determine that for each of one, one or more (e.g., two, three, or four), or one or more haplotypes (e.g., using their sequence reads), the likelihood of one copy of the wild-type gene haplotype is higher than the likelihood of two copies of the wild-type gene haplotype. The likelihood difference can be, for example, 1%, 2%, 3%, 5%, 10%, 15%, 20%, or more. The computing system can determine that the copy number of the wild-type gene haplotype is 1. If the likelihood of one copy of the wild-type gene haplotype is greater than the likelihood of two copies of the wild-type gene haplotype, then the copy number of the wild-type gene haplotype can be one.

[0160] In some embodiments, the computing system can determine that for each of one or more pairs (or all pairs) of consecutive gene / gene paralog distinguishing bases of a plurality of gene / gene paralog distinguishing bases, where a first haplotype of the one or more haplotypes includes a gene base in the consecutive gene / gene paralog distinguishing bases and a second haplotype of the one or more haplotypes includes a gene base and a gene paralog base (or a gene paralog base and a gene base) in the consecutive gene / gene paralog distinguishing bases, the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype.

[0161] The first haplotype may include gene bases in consecutive gene / gene paralogue discriminating bases. The second haplotype may include a gene base to gene paralogue base (or gene paralogue base to gene base) transition between consecutive gene / gene paralogue discriminating bases. The consecutive gene / gene paralogue discriminating bases are consecutive within the plurality of gene / gene paralogue discriminating bases, regardless of whether the consecutive gene / gene paralogue discriminating bases are adjacent bases in the reference sequence. For example, the gene / gene paralogue discriminating bases at sites 2 and 3 (or discriminating base positions) in the following examples are consecutive gene / gene paralogue discriminating bases, regardless of whether the consecutive gene / gene paralogue discriminating bases are adjacent bases in the reference sequence. The computing system may combine (e.g., average or weight average) the likelihoods determined for each of one or more pairs (or all pairs) of consecutive gene / gene paralogue discriminating bases to determine that the likelihood of one copy of a wild-type gene haplotype is higher than the likelihood of two copies of a wild-type gene haplotype. The likelihood of one copy of a wild-type gene haplotype may include the sum of the likelihoods of one copy of a wild-type gene haplotype determined for each of one or more pairs of consecutive gene / gene paralog discriminating bases. The likelihood of two copies of a wild-type gene haplotype may include the sum of the likelihoods of two copies of a wild-type gene haplotype determined for each of one or more pairs of consecutive gene / gene paralog discriminating bases.

[0162] The computing system can determine that (1) given the number of sequence reads of the second plurality of sequence reads, each of which includes a gene base in consecutive gene / gene paralog discriminating bases, the likelihood of one copy of the wild type gene haplotype is higher than the likelihood of two copies of the wild type gene haplotype for one or more (or all) pairs of consecutive gene / gene paralog discriminating bases. The computing system can determine that (2) given the number of sequence reads of the second plurality of sequence reads, each of which includes a gene base and a gene paralog base in consecutive gene / gene paralog discriminating bases, and / or the number of sequence reads of the second plurality of sequence reads, each of which includes a gene base and a gene base in consecutive gene / gene paralog discriminating bases, the likelihood of one copy of the wild type gene haplotype is higher than the likelihood of two copies of the wild type gene haplotype for one or more (or all) pairs of consecutive gene / gene paralog discriminating bases. In some embodiments, the computing system can determine that (4) given a number of sequence reads of the second plurality of sequence reads, each comprising a gene paralogue base in consecutive gene / gene paralogue discriminating bases, for each of one or more pairs (or all pairs) of consecutive gene / gene paralogue discriminating bases, the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype.

[0163] For example, by analyzing linkage information between gene / gene paralog discriminating bases, the following haplotypes (at six sites for illustrative purposes only) can be determined for a subject:

[0164] [Table 11] Conversion of gene bases to gene paralog bases occurs between sites 2 and 3 for haplotype 2. For example, the number of reads having a gene base at site 2 and site 3 is 103, the number of reads having a gene base at site 2 and a gene paralog base at site 3 is 99, and the number of reads having a gene paralog base at site 2 and site 3 is 210. The computing system can determine that the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype if the number of reads having a gene base at site 2 and site 3 is 103 and the number of reads having a gene base at site 2 and a gene paralog base at site 3 is 99. The computing system can determine that the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype without using reads from the wildtype gene paralog haplotype.

[0165] Continuing with the example, conversion of gene bases to gene paralogue bases occurs between site 3 and site 4 for haplotype 2. For example, the number of reads having a gene base at site 3 and site 4 is 100, the number of reads having a gene paralogue base at site 3 and a paralogue base at site 4 is 99, and the number of reads having a gene paralogue base at site 3 and site 4 is 190. The computing system can determine that the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype if the number of reads having a gene base at site 3 and site 4 is 100 and the number of reads having a gene paralogue base at site 3 and a gene base at site 4 is 99. The computing system can determine that the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype without using reads from the wildtype gene paralogue haplotype.

[0166] In some embodiments, the computing system can determine that the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype by combining (e.g., averaging or weighted averaging) (1) that the likelihood of one copy of the wildtype gene haplotype determined taking into account the number of reads having a gene base at site 2 and site 3, and the number of reads having a gene base at site 2 and a gene paralogous base at site 3, is higher than the likelihood of two copies of the wildtype gene haplotype, and (2) that the likelihood of one copy of the wildtype gene haplotype determined taking into account the number of reads having a gene base at site 3 and site 4, and the number of reads having a gene paralogous base at site 3 and a gene base at site 4, is higher than the likelihood of two copies of the wildtype gene haplotype.

[0167] In some embodiments, the computing system can determine that (1) the likelihood of one copy of the wild-type gene haplotype is higher than the likelihood of two copies of the wild-type gene haplotype given the total number of sequence reads for each of one or more pairs (or all pairs) of consecutive gene / gene paralogue discriminating bases of the plurality of gene / gene paralogue discriminating bases, where a first haplotype of the one or more haplotypes includes a gene base in consecutive gene / gene paralogue discriminating bases. The computing system can determine that (2) the likelihood of one copy of the wild-type gene haplotype is higher than the likelihood of two copies of the wild-type gene haplotype based on the total number of sequence reads for each of one or more pairs (or all pairs) of consecutive gene / gene paralogue discriminating bases of the plurality of gene / gene paralogue discriminating bases, where a second haplotype of the one or more haplotypes includes a gene base and a gene paralogue base (or a gene paralogue base and a gene base) in consecutive gene / gene paralogue discriminating bases. In the above example, the computing system can determine that given (1) the total number of reads that have a gene base at site 2 and site 3 and a gene base at site 3 and a gene base at site 4, and (2) the number of reads that have a gene base at site 2 and a gene paralogous base at site 3 and a gene paralogous base at site 3 and a gene base at site 4, the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype.

[0168] For example, by analyzing linkage information between gene / gene paralog discriminator bases, the following haplotypes (for illustrative purposes only, at six bases or sites or positions of the discriminator base) can be determined for a subject:

[0169] [Table 12] Conversion of gene bases to gene paralog bases occurs between sites 2 and 3 for haplotype 2. For example, the number of reads having a gene base at sites 2 and 3 is 103, the number of reads having a gene base at site 2 and a gene paralog base at site 3 is 99, the number of reads having a gene paralog base at site 2 and a gene base at site 3 is 90, and the number of reads having a gene paralog base at site 2 and a gene paralog base at site 3 is 104. The computing system can determine that the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype if the number of reads having a gene base at site 2 and site 3 is 103 and the number of reads having a gene base at site 2 and a gene paralog base at site 3 is 99. The computing system can determine that the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype without using reads from the wildtype gene paralog haplotype. The computing system can determine, without using reads having a gene paralog base at site 2 and a gene base at site 3, that the likelihood of one copy of the wild-type gene haplotype is higher than the likelihood of two copies of the wild-type gene haplotype because that haplotype (haplotype 3) has most gene paralog bases at the differentiating base sites / positions and is therefore less likely / more likely to be a gene variant haplotype.

[0170] The copy number of the wild type gene haplotype can be 1. The computing system can determine that the subject is a carrier of a gene variant haplotype. The one or more haplotypes can include four haplotypes (e.g., one copy of the wild type gene haplotype, one copy of the gene variant haplotype, one copy of the gene paralog wild type haplotype, and one copy of a haplotype having a gene paralog base at a high percentage (such as 80%, 85%, 90%, 95% or more) of the discriminating bases), and thus unlikely to be a gene variant haplotype / probable to be a gene paralog variant haplotype. The total copy number of the gene and gene paralog can be 4. The copy number of each of the four haplotypes can be 1. The computing system can determine the gene variant status of the subject as a carrier of a gene variant haplotype. For example, by analyzing the linkage information between the gene / gene paralog discriminating bases, the following haplotypes (at six sites for illustrative purposes only) can be determined for the subject:

[0171] [Table 13] The computing system can determine, given a number of reads having a gene base at site 2 and site 3 (haplotype 1), and a number of reads having a gene base at site 2 and a gene paralogous base at site 3 (haplotype 2), that the likelihood of one copy of the wild type gene base at site 2 and site 3 is higher than the likelihood of two copies of the wild type gene base at site 2 and site 3. The computing system can determine, given a number of reads having a gene base at site 4 and site 5 (haplotype 1), and a number of reads having a gene paralogous base at site 4 and a gene base at site 5 (haplotype 2), that the likelihood of one copy of the wild type gene haplotype base at site 4 and site 5 is higher than the likelihood of two copies of the wild type gene base at site 4 and site 5. The computing system can determine that given the number of reads with a gene base at site 2 and site 3 (haplotype 1) and the number of reads with a gene base at site 2 and a gene paralog base at site 3 (haplotype 2), and / or the number of reads with a gene base at site 4 and site 5 (haplotype 1) and the number of reads with a gene paralog base at site 4 and a gene base at site 5 (haplotype 2), the likelihood of one copy of a wild-type gene haplotype is higher than the likelihood of two copies of a wild-type gene haplotype. The subject is a carrier of a gene variant haplotype because the subject has one copy of a wild-type gene haplotype and one copy of a gene variant haplotype (and one copy of a wild-type gene paralog haplotype and one copy of a haplotype with most gene paralog bases at the discriminatory base site / position, and therefore less likely to be a gene variant haplotype / more likely to be a gene paralog variant haplotype).

[0172] The one or more haplotypes may include three haplotypes. The total copy number of the gene and the gene paralog gene may be four. The copy numbers of the wild type gene haplotype, the gene variant haplotype, and the wild type gene paralog haplotype (or haplotypes having gene paralog bases at a high percentage of discriminating bases, such as 80%, 85%, 90%, 95%, or more) may be one, one, and two, respectively. The computing system may determine the genetic status of the subject as a carrier of a gene variant haplotype. For example, by analyzing the linkage information between the gene / gene paralog discriminating bases, the following haplotypes (at six sites for illustrative purposes only) may be determined for the subject:

[0173] [Table 14] The computing system can determine, given a number of reads having a gene base at site 2 and site 3 (haplotype 1), and a number of reads having a gene base at site 2 and a gene paralogous base at site 3 (haplotype 2), that the likelihood of one copy of the wild type gene base at site 2 and site 3 is higher than the likelihood of two copies of the wild type gene base at site 2 and site 3. The computing system can determine, given a number of reads having a gene base at site 4 and site 5 (haplotype 1), and a number of reads having a gene paralogous base at site 4 and a gene base at site 5 (haplotype 2), that the likelihood of one copy of the wild type gene haplotype base at site 4 and site 5 is higher than the likelihood of two copies of the wild type gene base at site 4 and site 5. The computing system can determine that given the number of reads having a gene base at site 2 and site 3 (haplotype 1), and the number of reads having a gene base at site 2 and a gene paralogous base at site 3 (haplotype 2), and / or the number of reads having a gene base at site 4 and site 5 (haplotype 1), and the number of reads having a gene paralogous base at site 4 and a gene base at site 5 (haplotype 2), the likelihood of one copy of the wildtype gene haplotype is higher than the likelihood of two copies of the wildtype gene haplotype. Because the subject has one copy of the wildtype gene haplotype and one copy of the gene variant haplotype, the subject is a carrier of the gene variant haplotype.

[0174] Compound heterozygosity. The one or more haplotypes may include two or more haplotypes. None of the two or more haplotypes may include a gene base at each of the multiple gene / gene paralogue distinguishing bases. Each of the two or more haplotypes may not include a gene base at all of the multiple gene / gene paralogue distinguishing bases. Each of the two or more haplotypes may include a gene paralogue base at one or more of the multiple gene / gene paralogue distinguishing bases. The computing system can determine that the subject is a compound heterozygosity of a gene variant haplotype. For example, by analyzing the linkage information between the gene / gene paralogue distinguishing bases, the following haplotypes (at six sites for illustrative purposes only) may be determined for the subject:

[0175] [Table 15] Because the subject does not have any copies of the wild-type gene haplotype, and has one copy of each of the two gene variant haplotypes, the subject is a compound heterozygote for the gene variant haplotypes.

[0176] Homozygous. One or more haplotypes may include an identical base (e.g., a gene base or a gene paralog base) at a gene / gene paralog distinguishing base, or at each of two or more of the plurality of gene / gene paralog distinguishing bases. The computing system can determine that the subject is homozygous (e.g., homozygous for a wild type gene paralog haplotype or homozygous for a gene variant haplotype) at one or more of the plurality of gene / gene paralog distinguishing bases.

[0177] The one or more haplotypes may include only one haplotype. The only haplotype may not include a gene base at one, one or more, or each of the multiple gene / gene paralog distinguishing bases. The only haplotype may include a gene paralog base at one, one or more, or each of the multiple gene / gene paralog distinguishing bases. The computing system can determine that the subject is homozygous for a gene variant haplotype at one, one or more, or each of the multiple gene / gene paralog distinguishing bases. For example, based on the multiple gene / gene paralog distinguishing bases, the computing system can determine the copy number (CN) of the gene base. The number of reads supporting the gene base or gene paralog base, and the total CN of the gene and gene paralog can be used to determine the most likely combination of the CN of the gene base and the CN of the gene paralog base. If the CN of the gene base is determined to be 0, this indicates that the subject does not have a copy of the wild type gene haplotype (a haplotype carrying a gene A base at the modification site of interest) and is homozygous for the gene variant haplotype.

[0178] The method 600 ends at block 624.

[0179] Execution environment FIG. 7 illustrates a general architecture of an exemplary computing device 700 configured to determine or identify one or more genetic variants (e.g., GBA variants, CYP21A2 variants) or genetic variant status (e.g., carrier, compound heterozygous, or homozygous). The general architecture of the computing device 700 illustrated in FIG. 7 includes an arrangement of computer hardware and software components. The computing device 700 may include more (or fewer) elements than those illustrated in FIG. 7. However, not all of these general conventional elements need be illustrated to provide a useful disclosure. As illustrated, the computing device 700 includes a processing unit 710, a network interface 720, a computer-readable medium drive 730, an input / output device interface 740, a display 750, and an input device 760, all of which may communicate with each other via a communication bus. The network interface 720 may provide connectivity to one or more networks or computing systems. The processing unit 710 may thus receive information and instructions from other computing systems or services via a network. The processing unit 710 may also communicate with memory 770 and further provide output information for an optional display 750 via an input / output device interface 740. The input / output device interface 740 may also receive input from any input device 760, such as a keyboard, mouse, digital pen, microphone, touch screen, gesture recognition system, voice recognition system, game pad, accelerometer, gyroscope, or other input device.

[0180] Memory 770 may include computer program instructions (grouped into modules or components in some embodiments) that processing unit 710 executes to implement one or more embodiments. Memory 770 generally comprises RAM, ROM, and / or other persistent, secondary, or non-transitory computer-readable media. Memory 770 may store an operating system 772 that provides computer program instructions for use by processing unit 710 in the overall management and operation of computing device 700. Memory 770 may further comprise computer program instructions and other information for implementing aspects of the present disclosure.

[0181] For example, in one embodiment, the memory 770 includes a genetic variant or genetic variant status determination module 774 for determining or identifying one or more genetic variants or genetic variant status (e.g., carrier, compound heterozygous, or homozygous), such as method 400 described with reference to Figure 4, method 500 described with reference to Figure 5, or method 600 described with reference to Figure 6. Additionally, the memory 770 may include or be in communication with a data store 790 and / or one or more other data stores that store processed sequence reads, determined read counts, Gaussian mixture models, determined recombinant variants, determined recombinant variant copy numbers, or determined genetic variant status.

[0182] Additional Considerations In at least some of the foregoing embodiments, one or more elements used in one embodiment may be used interchangeably in another embodiment, except where such an exchange is not technically feasible. Those skilled in the art will appreciate that various other omissions, additions, and modifications may be made to the methods and structures described above without departing from the scope of the claimed subject matter. All such modifications and variations are intended to be included within the scope of the subject matter, as defined by the appended claims.

[0183] Those skilled in the art will appreciate that, for this and other processes and methods disclosed herein, the functions performed in the processes and methods may be performed in differing orders. Moreover, the outlined steps and operations are provided only as examples, and some of the steps and operations may be optional, may be combined into fewer steps and operations, or may be expanded to additional steps and operations without departing from the essence of the disclosed embodiments.

[0184] For the use of substantially any plural and / or singular term herein, one of ordinary skill in the art may substitute the plural for the singular and / or the singular for the plural as appropriate to the context and / or application. For clarity, various singular / plural permutations may be expressly set forth herein. As used herein and in the appended claims, the singular forms "a," "an," and "the" include plural referents unless the context clearly dictates otherwise. Thus, phrases such as "a device configured to" are intended to include one or more enumerated devices. Such one or more enumerated devices may also be collectively configured to execute the recited detailed description. For example, "a processor configured to execute detailed description A, B, and C" may include a first processor configured to execute detailed description A and to perform operations in conjunction with a second processor configured to execute detailed description B and C. Any reference to "or" herein is intended to encompass "and / or" unless otherwise indicated.

[0185] In general, those of skill in the art will understand that the terms used herein, and particularly in the appended claims (e.g., the body of the appended claims), are generally intended as "open" terms (e.g., the term "including" should be interpreted as "including but not limited to," the term "having" should be interpreted as "having at least," the term "includes" should be interpreted as "includes but is not limited to," etc.). Those of skill in the art will further understand that where a specific number of introduced claim recitations are intended, such intent will be expressly recited in the claim, and in the absence of such recitation, no such intent exists. For example, to aid in understanding, the following appended claims may include the use of the introductory phrases "at least one" and "one or more" to introduce the claim recitations. However, the use of such phrases should not be construed as limiting any particular claim that includes such an introduced claim statement to embodiments that include only one of such statements, even when the same claim includes "one or more" or "at least one" and an indefinite article such as "a" or "an" (e.g., "a" and / or "an" should be interpreted to mean "at least one" or "one or more"), and the same applies to the use of indefinite articles used to introduce claim statements. Moreover, even when a specific number of introduced claim statements is explicitly recited, one of skill in the art will recognize that such a recitation should be interpreted in the sense of at least the recited number (e.g., the unmodified recitation of "two statements" means, in the absence of other modifications, at least two statements or more than two statements).Furthermore, when phrases similar to "at least one of A, B, and C, etc." are used, generally such structures are intended in the sense that one of ordinary skill in the art would understand the phrase (e.g., "a system having at least one of A, B, and C" includes, but is not limited to, systems having only A, only B, only C, both A and B, both A and C, both B and C, and / or both A, B, and C, etc.). When phrases similar to "at least one of A, B, or C, etc." are used, generally such structures are intended in the sense that one of ordinary skill in the art would understand the phrase (e.g., "a system having at least one of A, B, or C" includes, but is not limited to, systems having only A, only B, only C, both A and B, both A and C, both B and C, and / or both A, B, and C, etc.). Those skilled in the art will further appreciate that virtually any disjunction and / or phrase presenting two or more alternative terms, whether in the specification, claims, or drawings, should be understood to contemplate the possibility of including one of the terms, either of the terms, or both terms. For example, the phrase "A or B" will be understood to include the possibilities of "A" or "B" or "A and B."

[0186] Additionally, where features or aspects of the disclosure are described in terms of a Markush group, one of skill in the art will thereby recognize that the disclosure is also described in terms of any individual members or subgroups of the Markush group members.

[0187] As will be understood by those skilled in the art, for any and all purposes, such as in terms of providing a written description, all ranges disclosed herein also encompass any and all possible subranges and combinations of those subranges. Any recited range can be readily recognized as fully descriptive and allowing the same range to be broken down into at least equal halves, thirds, quarters, fifths, tenths, etc. As a non-limiting example, each range described herein can be readily broken down into a lower third, a middle third, and an upper third, etc. Also, as will be understood by those skilled in the art, all language such as "up to," "at least," "greater than," "less than" and the like refers to a range that includes the recited numbers and can then be broken down into subranges as described above. Finally, as will be understood by those skilled in the art, a range includes each individual component. Thus, for example, a group having 1 to 3 items means a group having 1, 2, or 3 items. Similarly, a group having 1-5 items means a group having 1, 2, 3, 4, or 5 items, etc.

[0188] Various embodiments of the present disclosure have been described herein for illustrative purposes, and it will be understood that various modifications may be made without departing from the scope and spirit of the present disclosure. Accordingly, the various embodiments disclosed herein are not intended to be limiting, the true scope and spirit of which is set forth by the following claims.

[0189] It is to be understood that not necessarily all objectives or advantages may be achieved in accordance with any particular embodiment described herein. Thus, for example, those skilled in the art will recognize that a particular embodiment may be configured to operate in a manner that achieves or optimizes an advantage or set of advantages as taught herein, without necessarily achieving other objectives or advantages that may be taught or suggested herein.

[0190] All of the processes described herein may be embodied in, and may be fully automated via, software code modules executed by a computing system including one or more computers or processors. The code modules may be stored on any type of non-transitory computer-readable medium or other computer storage device. Some or all of the methods may be embodied in dedicated computer hardware.

[0191] Many other variations beyond those described herein will be apparent from this disclosure. For example, depending on the embodiment, certain acts, events, or functions of any of the algorithms described herein may be performed in a different order, added, combined, or omitted entirely (e.g., not all described acts or events are necessary to implement the algorithm). Furthermore, in certain embodiments, acts or events may be performed simultaneously rather than sequentially, e.g., via multi-threading, interrupt processing, or multiple processors or processor cores, or on other parallel systems. In addition, different tasks or processes may be performed by different machines and / or computing systems that can function together.

[0192] The various example logic blocks and modules described in connection with the embodiments disclosed herein may be implemented or performed by mechanical devices such as processing units or processors, digital signal processors (DSPs), application specific integrated circuits (ASICs), field programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic, discrete hardware components, or any combination thereof designed to perform the functions described herein. The processor may be a microprocessor, but alternatively the processor may be a controller, microcontroller, or state machine, combinations thereof, and the like. The processor may include electrical circuitry configured to process computer-executable instructions. In another embodiment, the processor includes an FPGA or other programmable device that performs logical operations without processing computer-executable instructions. The processor may also be implemented as a combination of computing devices, such as a combination of a DSP and a microprocessor, multiple microprocessors, one or more microprocessors in association with a DSP core, or any other such configuration. Although primarily digital technology is described herein, the processor may also include primarily analog components. For example, some or all of the signal processing algorithms described herein may be implemented in analog circuitry or mixed analog and digital circuitry. The computing environment may include any type of computer system, including, but not limited to, a computer system based on a computing engine within a microprocessor, mainframe computer, digital signal processor, portable computing device, device controller, or appliance, to name a few.

[0193] Any process illustrations, elements or blocks in the flow diagrams described herein and / or shown in the accompanying drawings should be understood as potentially representing modules, segments or portions of code that include one or more executable instructions for implementing a particular logical function or element in the process. Alternate embodiments are included within the scope of the embodiments described herein, and depending on the functionality involved, elements or functions may be omitted or removed from the order from that shown or discussed, including substantially simultaneously or in reverse order, as will be understood by those skilled in the art.

[0194] It should be emphasized that many variations and modifications may be made to the above-described embodiments, and that the elements thereof should be understood to be among the other acceptable examples. All such modifications and variations are intended to be included herein within the scope of this disclosure and protected by the following claims.

Claims

A system for determining a GBA state, comprising: A non-transitory memory configured to store executable instructions; and A hardware processor in communication with the non-transitory memory, the hardware processor being caused by the executable instructions to Receive a first plurality of sequence reads generated from a sample obtained from a subject; Align the first plurality of sequence reads to a reference genome sequence to obtain a second plurality of sequence reads aligned to a GBA gene or a GBAP1 gene in the reference genome sequence; Determine the number of sequence reads of the second plurality of sequence reads aligned to a unique region between the GBA gene and the GBAP1 gene in the reference genome sequence; Determine the normalized number of sequence reads of the sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference genome sequence; Given the normalized number of sequence reads aligned to a region between the GBA gene and the GBAP1 gene, determine the total copy number of the GBA gene and the GBAP1 gene using a mixture Gaussian model that includes a plurality of Gaussians each representing a different integer copy number; Phase one or more haplotypes derived from the GBA gene or the GBAP1 gene in a region of the GBA gene or the corresponding region of the GBAP1 gene that includes a plurality of GBA / GBAP1 discriminant bases, using the sequence reads of the second plurality of sequence reads aligned to the region or the corresponding region that includes the plurality of GBA / GBAP1 discriminant bases; Determine the copy number of each of the one or more haplotypes using the total copy number of the GBA gene and the GBAP1 gene and the number of sequence reads of the second plurality of sequence reads each including one or more of the plurality of GBA / GBAP1 discriminant bases that support the haplotype; and Determine the GBA state of the subject using the one or more haplotypes derived from the GBA gene or the GBAP1 gene in a region of the GBA gene or the corresponding region of the GBAP1 gene, and / or the copy number of each of the one or more haplotypes, the system being programmed to perform the foregoing. Claim 2 The unique region between the GBA gene and the GBAP1 gene in the reference genomic sequence includes a unique region having a length of about 10 kilobases, and / or the unique region between the GBA gene and the GBAP1 gene in the reference genomic sequence includes chr1:155220429-155230539 of hg38 or a corresponding region of the reference human genomic sequence. The system according to claim 1.

3. Determining the normalized number of the sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference genomic sequence includes (1a) the depth of the sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene, (1b) the length of the unique region, (2a) the depth of the sequence reads of the first plurality of sequence reads aligned to each of a plurality of regions in the reference genomic sequence other than the locus including the GBA gene and the GBAP1 gene, and (2b) the length of each of the plurality of regions of the reference genome other than the locus including the GBA gene and the GBAP1 gene, and using the above to determine the normalized number of the sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference genomic sequence. The system according to claim 1. **Claim 4**: The hardware processor, by the executable instructions, determines a normalized and corrected number of the sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference genome sequence from the normalized number of the sequence reads aligned to the unique region between the GBA gene and the GBAP1 gene in the reference genome sequence, using (1) the GC content of the unique region between the GBA gene and the GBAP1 gene and, optionally, (2) the GC content of each of one or more regions of the reference genome sequence other than the locus containing the GBA gene and the GBAP1 gene, and determining the total copy number of the GBA gene and the GBAP1 gene includes determining the total copy number of the GBA gene and the GBAP1 gene using the mixture Gaussian model with the normalized and corrected number of the sequence reads aligned to the region between the GBA gene and the GBAP1 gene given, and the system according to claim 1, further programmed to perform this. **Claim 5** Determining the total copy number of the GBA gene and the GBAP1 gene includes determining the copy number of the region between the GBA gene and the GBAP1 gene using the mixture Gaussian model with the normalized number of the sequence reads aligned to the region between the GBA gene and the GBAP1 gene given, and the total copy number of the GBA gene and the GBAP1 gene is the copy number of the region between the GBA gene and the GBAP1 gene + 2, and the system according to claim 1. **Claim 6** Determining the total copy number of the GBA gene and the GBAP1 gene includes determining the total copy number of the GBA gene and the GBAP1 gene using the mixture Gaussian model and a predetermined posterior probability threshold with the normalized number of the sequence reads aligned to the region between the GBA gene and the GBAP1 gene given, and optionally, the predetermined posterior probability threshold is 0.95, and the system according to claim 1. **Claim 7** Phasing the one or more haplotypes derived from the GBA gene or the GBAP1 gene involves analyzing the linkage information between the GBA / GABP1 discriminatory bases of the plurality of GBA / GABP1 discriminatory bases using the sequence reads of the second plurality of sequence reads aligned to the region containing the plurality of GBA / GABP1 discriminatory bases or the corresponding region. The system according to claim 1.

8. Phasing the one or more haplotypes derived from the GBA gene or the GBAP1 gene involves phasing the one or more haplotypes derived from the GBA gene or the GBAP1 gene using the sequence reads of the second plurality of sequence reads aligned to two or more of the plurality of GBA / GABP1 discriminatory bases respectively. The system according to claim 1.

9. The sequence reads of the second plurality of sequence reads are aligned to the region of the GBA gene containing the plurality of GBA / GABP1 discriminatory bases or the corresponding region of the GBAP1 gene with an alignment quality score of 0 or more. The system according to claim 1.

10. The region of the GBA gene containing the plurality of GBA / GABP1 discriminatory bases or the corresponding region of the GBAP1 gene is about 1.1 kilobases in length. The region of the GBA gene containing the plurality of GBA / GABP1 discriminatory bases or the corresponding region of the GBAP1 gene contains exons 9-11 of the GBA gene or the GBAP1 gene respectively, and / or The region of the GBA gene containing the plurality of GBA / GABP1 discriminatory bases or the corresponding region of the GBAP1 gene contains p.L483P, p.D448H, c.1263del, RecNciI, RecTL, and c.1263del+RecTL. The system according to claim 1.

11. The plurality of GBA / GABP1 discriminatory bases contains 10 GBA / GABP1 discriminatory bases. The system according to claim 1.

12. The one or more haplotypes include a wild-type GBA haplotype, a wild-type GBAP1 haplotype, and / or a GBA / GBAP1 hybrid haplotype. Optionally, the GBA / GBAP1 hybrid haplotype includes a GBA mutant haplotype or a GBAP1 mutant haplotype. The system according to claim 1.

13. determining the copy number of each of the one or more haplotypes, determining that the likelihood of one copy of the wild-type GBA haplotype is higher than the likelihood of two copies of the wild-type GBA haplotype, given the number of sequence reads of the second plurality of sequence reads, each of which includes one or more of the plurality of GBA / GBA P1 discriminatory bases that support the wild-type GBA haplotype, determining that the copy number of the wild-type GBA haplotype is 1, the system of claim 1, comprising:

14. The system of claim 1, wherein the copy number of the wild-type GBA haplotype is 1 and the GBA status of the subject comprises a carrier of a GBA variant haplotype.

15. The system of claim 1, wherein the one or more haplotypes include four haplotypes, the total copy number of the GBA gene and the GBA P1 gene is 4, the copy number of each of the four haplotypes is 1, and the GBA status of the subject comprises a carrier of a GBA variant haplotype.

16. The system of claim 1, wherein the one or more haplotypes include two or more GBA variant haplotypes, none of the two or more GBA variant haplotypes includes a GBA base at each of the plurality of GBA / GBA P1 discriminatory bases, and the GBA status of the subject comprises a compound heterozygosity of a GBA variant haplotype.

17. The system of claim 1, wherein the hardware processor is further programmed by the executable instructions to determine that the copy number of the GBA base at each of one or more of the plurality of GBA / GBA P1 discriminatory bases that are not the GBA base is 0, using the sequence reads of the second plurality of sequence reads, each of which includes a base at the GBA / GBA P1 discriminatory base that is not the GBA base, optionally, the base at the GBA / GBA P1 discriminatory base that is not the GBA base is a GBA P1 base, and optionally, determining the GBA status comprises determining that the subject is homozygous at each of the one or more of the plurality of GBA / GBA P1 discriminatory bases. The system according to claim 1, wherein the hardware processor is further programmed to execute, by the executable instructions, to generate a user interface (UI) including a UI element representing or including a GBA state.

19. A method for determining a CYP21A2 state, comprising: Receiving, under the control of a hardware processor, a first plurality of array reads generated from a sample obtained from a subject; Aligning the first plurality of array reads to a reference genome sequence to obtain a second plurality of array reads aligned to a CYP21A2 gene or a CYP21A1P pseudogene in the reference genome sequence; Determining the number of array reads of the second plurality of array reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference genome sequence; Determining a normalized number of the array reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene in the reference genome sequence; Given the normalized number of array reads aligned to the CYP21A2 gene or the CYP21A1P pseudogene, determining the total copy number of the CYP21A2 gene and the CYP21A1P pseudogene using a mixture of Gaussians model including a plurality of Gaussians each representing a different integer copy number; Phasing one or more haplotypes derived from the CYP21A2 gene or the CYP21A1P pseudogene in a region of the CYP21A2 gene or a corresponding region of the CYP21A1P pseudogene including a plurality of CYP21A2 / CYP21A1P discriminant bases, using the array reads of the second plurality of array reads aligned to the region or the corresponding region including the plurality of CYP21A2 / CYP21A1P discriminant bases; Determining the copy number of each of the one or more haplotypes using the total copy number of the CYP21A2 gene and the CYP21A1P gene and the number of array reads of the second plurality of array reads each including one or more of the plurality of CYP21A2 / CYP21A1P discriminant bases supporting the haplotype. ​ A method comprising determining the CYP21A2 status of the subject using one or more haplotypes derived from the CYP21A2 gene or the CYP21A1P pseudogene in the region of the CYP21A2 gene or the corresponding region of the CYP21A1P pseudogene, and / or the copy number of each of the one or more haplotypes.

20. A system for determining a recombinant mutant, comprising a non-transitory memory configured to store executable instructions and a first plurality of sequence reads generated from a sample obtained from a subject, a hardware processor in communication with the non-transitory memory, wherein the hardware processor, by the executable instructions, aligns the first plurality of sequence reads to a reference sequence to obtain a second plurality of sequence reads aligned to a gene or gene paralog in the reference sequence, or a region between them; given the number of sequence reads aligned to the gene or gene paralog, or a region between them, determines the total copy number of the gene and the gene paralog using a mixture of Gaussians including a plurality of Gaussians each representing a different integer copy number; phasing one or more haplotypes derived from the gene (including recombinant mutants of the gene) or the gene paralog, or a region of the gene or the corresponding region of the gene paralog, containing a plurality of gene / gene paralog identification bases, using the sequence reads of the second plurality of sequence reads aligned to the region or the corresponding region containing the plurality of gene / gene paralog identification bases; a system programmed to execute determining the copy number of each of the one or more haplotypes using the total copy number of the gene and the gene paralog and the number of sequence reads of the second plurality of sequence reads each containing one or more of the plurality of gene / gene paralog identification bases supporting the haplotype.

Images