Methods and Systems for Diagnosis from Whole Genome Sequencing Data

Abstract

Systems, devices, computer-readable media, and methods disclosed herein include for paralog genotyping such as determining the copy number of the survival motor neuron 1 gene and genotyping cytochrome P450 family 2 subfamily D member 6 gene using a Gaussian mixture model comprising a plurality of Gaussian functions each representing a different integer copy number.

Description

[0001] Cross - reference to related applications

[0002] This application claims the priority benefits of U.S. Provisional Patent Application No. 62 / 896,548, filed on September 5, 2019; U.S. Provisional Patent Application No. 62 / 908,555, filed on September 30, 2019; and U.S. Provisional Patent Application No. 63 / 006,651, filed on April 7, 2020. The entire content of each of the related applications is incorporated herein by reference. Background of the Invention Field of the Invention

[0003] The present disclosure generally relates to the field of paralog gene genotyping and, more particularly, to paralog gene genotyping using sequencing data. Background of the Invention

[0005] Genotyping is challenging. For example, spinal muscular atrophy is caused by the loss of function of the survival motor neuron 1 (SMN1) gene while the paralogous SMN2 gene is retained. Analyzing this region has been challenging because the sequences of SMN1 and its paralog SMN2 are nearly identical. Also, CYP2D6 is involved in the metabolism of 25% of all drugs. Genotyping CYP2D6 is challenging due to its high polymorphism, the presence of common structural variants (SVs), and its high sequence similarity to the pseudogene paralog CYP2D7 of the gene. Summary of the Invention

[0006] Disclosed herein are methods for determining the copy number of the survival motor neuron 1 (SMN1) gene. In some embodiments, a method for determining the copy number of the SMN1 gene is under the control of a processor, such as a hardware processor or a virtual processor, and includes: receiving sequence data that includes a plurality of sequence reads obtained from a sample of a subject and aligned to the SMN1 gene or the survival motor neuron 2 (SMN2) gene. The method may include: determining (i) a first number of sequence reads of the plurality of sequence reads that are aligned to a first SMN1 or SMN2 region that respectively includes at least one of exons 1 to 6 of the SMN1 gene or the SMN2 gene and (ii) a second number of sequence reads of the plurality of sequence reads that are aligned to a second SMN1 or SMN2 region that respectively includes at least one of exons 7 and 8 of the SMN1 gene or the SMN2 gene. The method may include: using (i) the length of the first SMN1 or SMN2 region and (ii) the length of the second SMN1 or SMN2 region to determine (i) a first normalized number of sequence reads that are aligned to the first SMN1 or SMN2 region and (ii) a second normalized number of sequence reads that are aligned to the second SMN1 or SMN2 region. The method may include: using a Gaussian mixture model that includes a plurality of Gaussian functions each representing a different integer copy number to determine (i) the copy number of the total survival motor neuron (SMN) genes that are respectively a full-length SMN1 gene, a full-length SMN2 gene, a truncated SMN1 gene, or a truncated SMN2 gene and (ii) the copy number of any full-length SMN genes that are respectively a full-length SMN1 gene or a full-length SMN2 gene, taking into account (i) the first normalized number of sequence reads that are aligned to the first SMN1 or SMN2 region and (ii) the second normalized number of sequence reads that are aligned to the second SMN1 or SMN2 region. The method may include: for a base among a plurality of SMN1 gene-specific bases associated with the full-length SMN1 gene, determining the most likely combination among a plurality of possible combinations of the possible copy number of the SMN1 gene and the possible copy number of the SMN2 gene that respectively include a total number of copies of any full-length SMN genes determined, taking into account (a) the number of sequence reads of the plurality of sequence reads that have a base supporting the SMN1 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads that have a base supporting the SMN2 gene-specific base corresponding to the SMN1 gene-specific base of the SMN2 gene. The method may include: using the most likely combination of the possible copy number of the SMN1 gene and the possible copy number of the SMN2 gene determined for the SMN1 gene-specific base to determine the copy number of the SMN1 gene.

[0007] In some embodiments, the sequencing data includes whole genome sequencing (WGS) data or short-read WGS data. In some embodiments, the subject is a neonatal subject, a pediatric subject, an adolescent subject, or an adult subject. The sample may comprise cellular or cell-free DNA.

[0008] In some embodiments, the sequence reads of the plurality of sequence reads are aligned with a first SMN1 or SMN2 region or a second SMN1 or SMN2 region, wherein the alignment quality score is about zero. The first SMN1 or SMN2 region may respectively comprise exons 1 to 6 of the SMN1 gene or the SMN2 gene and has a length of about 22.2 kb. The second SMN1 or SMN2 region may respectively comprise exons 7 and 8 of the SMN1 gene or the SMN2 gene and has a length of about 6 kb.

[0009] In some embodiments, determining (i) a first normalized number of sequence reads aligned to a first SMN1 or SMN2 region and (ii) a second normalized number of sequence reads aligned to a second region comprises: determining (i) the first normalized number of sequence reads aligned to the first SMN1 or SMN2 region and (ii) the second normalized number of sequence reads aligned to the second SMN1 or SMN2 region using, respectively, (i) the length of the first SMN1 or SMN2 region and (ii) the length of the second SMN1 or SMN2 region, and determining (iii) the depth of sequence reads in a region of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene in the sequence data. Determining (i) the first normalized number of sequence reads aligned to the first SMN1 or SMN2 region and (ii) the second normalized number of sequence reads aligned to the second SMN1 or SMN2 region may comprise: determining (i) the first SMN1 or SMN2 region length-normalized number of sequence reads aligned to the first SMN1 or SMN2 region and (ii) the second SMN1 or SMN2 region length-normalized number of sequence reads aligned to the second SMN1 or SMN2 region using, respectively, (i) the length of the first SMN1 or SMN2 region and (ii) the length of the second SMN1 or SMN2 region. Determining (i) the first normalized number of sequence reads aligned to the first SMN1 or SMN2 region and (ii) the second normalized number of sequence reads aligned to the second SMN1 or SMN2 region may comprise: using the depth of sequence reads in a region of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene, and determining (i) the first normalized depth of sequence reads aligned to the first SMN1 or SMN2 region and (ii) the second normalized depth of sequence reads aligned to the second SMN1 or SMN2 region based on, respectively, (i) the first SMN1 or SMN2 region length-normalized number and (ii) the second SMN1 or SMN2 region length-normalized number, wherein the first normalized number of sequence reads aligned to the first SMN1 or SMN2 region and the second normalized number of sequence reads aligned to the second SMN1 or SMN2 region are the first normalized depth and the second normalized depth, respectively.

[0010] In some embodiments, determining (i) a first normalized number of sequence reads aligned to a first SMN1 or SMN2 region and (ii) a second normalized number of sequence reads aligned to a second region includes: using (i) the GC content of the first SMN1 or SMN2 region and (ii) the GC content of the second SMN1 or SMN2 region to determine (i) the first normalized number of sequence reads aligned to the first SMN1 or SMN2 region and (ii) the second normalized number of sequence reads aligned to the second SMN1 or SMN2 region, and determining (iii) the depth of sequence reads in a region of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene in the sequence data, and determining (iv) the GC content of the region of the genome.

[0011] In some embodiments, the depth of the region includes the average depth or median depth of sequence reads in a region of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene in the sequencing data. The region may comprise approximately 3000 preselected regions of the subject's genome, each of length approximately 2 kb. In some embodiments, (i) the first normalized number of sequence reads aligned to the first SMN1 or SMN2 region and / or (ii) the second normalized number of sequence reads aligned to the second SMN1 or SMN2 region is from about 30 to about 40.

[0012] In some embodiments, the Gaussian mixture model includes a one-dimensional Gaussian mixture model. The plurality of Gaussian functions of the Gaussian mixture model may represent integer copy numbers from 0 to 10. The mean of each Gaussian function in the plurality of Gaussian functions may be the integer copy number represented by the Gaussian function.

[0013] In some embodiments, determining (i) the copy number of the total SMN gene and (ii) the copy number of any full-length SMN gene includes using a Gaussian mixture model and a first predetermined posterior probability threshold to determine (i) the copy number of the total SMN gene and (ii) the copy number of any full-length SMN gene, taking into account (i) the first normalized number of sequence reads aligned to the first SMN1 or SMN2 region and (ii) the second normalized number of sequence reads aligned to the second SMN1 or SMN2 region, respectively. The first predetermined posterior probability threshold may be 0.95.

[0014] In some embodiments, the method includes: using (i) the determined copy number of the total SMN gene and (ii) the determined copy number of the full-length SMN gene to determine the copy number of the truncated SMN gene. The copy number of the truncated SMN gene may be the difference between (i) the determined copy number of the total SMN gene and (ii) the determined copy number of the full-length SMN gene.

[0015] In some embodiments, the SMN1 gene-specific base is a splicing enhancer. The SMN1 gene-specific base can be the base at c.840 of the SMN1 gene. In some embodiments, considering (a) the number of sequence reads of the plurality of sequence reads having bases that support the SMN1 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads having bases that support the corresponding SMN2 gene-specific base, the most likely combination of the possible copy number of the SMN1 gene and the possible copy number of the SMN2 gene is associated with the highest posterior probability relative to other combinations in the plurality of combinations.

[0016] In some embodiments, determining the most likely combination of the possible copy number of the SMN1 gene and the possible combination of the SMN2 gene includes: determining, considering the ratio of (a) the number of sequence reads of the plurality of sequence reads having bases that support the SMN1 gene-specific base to (b) the number of sequence reads of the plurality of sequence reads having bases that support the SMN2 gene-specific base corresponding to the SMN1 gene-specific base, the most likely combination among the plurality of possible combinations of the possible copy number of the SMN1 gene and the possible copy number of the SMN2 gene, each including a total number of copies of any complete SMN gene determined. Determining the most likely combination of the possible copy number of the SMN1 gene and the possible combination of the SMN2 gene can include: determining (a) the number of sequence reads of the plurality of sequence reads having bases that support the SMN1 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads having bases that support the SMN2 gene-specific base corresponding to the SMN1 gene-specific base; determining the ratio of (a) the number of sequence reads of the plurality of sequence reads having bases that support the SMN1 gene-specific base to (b) the number of sequence reads of the plurality of sequence reads having bases that support the SMN2 gene-specific base corresponding to the SMN1 gene-specific base; and determining, based on the ratio of (a) the number of sequence reads of the plurality of sequence reads having bases that support the SMN1 gene-specific base to (b) the number of sequence reads of the plurality of sequence reads having bases that support the SMN2 gene-specific base corresponding to the SMN1 gene-specific base, the most likely combination among the plurality of possible combinations of the possible copy number of the SMN1 gene and the possible copy number of the SMN2 gene, each including a total number of copies of any complete SMN gene determined.

[0017] In some embodiments, determining the most likely combination of the possible copy numbers of the SMN1 gene and the possible combinations of the SMN2 gene includes, for each of the plurality of SMN1 gene-specific bases, determining the most likely combination among the plurality of possible combinations of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene that each include a total number of copies of any complete SMN gene determined, taking into account (a) the number of sequence reads of the plurality of sequence reads having a base that supports the SMN1 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads having a base that supports the SMN2 gene-specific base corresponding to the SMN1 gene-specific base of the SMN2 gene, that is associated with the highest posterior probability. Determining the copy number of the SMN1 gene may include determining the copy number of the SMN1 gene based on the possible copy number of the SMN1 gene in the most likely combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene determined for each of the plurality of SMN1 gene-specific bases.

[0018] In some embodiments, the SMN1 gene-specific base is in agreement with each of the plurality of SMN1 gene-specific bases other than the SMN1 gene-specific bases that exceed a predetermined consistency threshold. The consistency threshold may be 97%. The plurality of SMN1 gene-specific bases may include 8 SMN1 gene-specific bases. Each of the plurality of SMN1 gene-specific bases may be located on intron 6, exon 7, intron 7, or exon 8 of the SMN1 gene. If the subject is of the first race, the plurality of SMN1 gene-specific bases may be different, if the subject is of the second race, the plurality of SMN1 gene-specific bases may be different, and if the race of the subject is unknown, the plurality of SMN1 gene-specific bases may be different. The race of the subject may be unknown and the plurality of SMN1 gene-specific bases may not be race-specific. The race of the subject may be known and the plurality of SMN1 gene-specific bases may be specific to the race of the subject. In some embodiments, the method includes receiving race information of the subject. The method may include selecting the plurality of SMN1 gene-specific bases from the plurality of SMN1 gene-specific bases based on the received race information.

[0019] In some embodiments, determining the copy number of the SMN1 gene includes: determining the copy number of the SMN1 gene and the copy number of the SMN2 gene using the most likely combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene determined for each of the bases specific to the plurality of SMN1 genes. Determining the copy number may include: determining the copy number of the SMN1 gene using the most likely combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene determined for the bases specific to the SMN1 gene and a second predetermined posterior probability threshold of the combinations of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene. The second predetermined posterior probability threshold may be 0.6 or 0.8.

[0020] In some embodiments, most of the possible copy numbers of the determined SMN1 gene are consistent. The copy number of the determined SMN1 gene may be the consistent possible copy number of the SMN1 gene. The method may include: determining the possible combinations of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene that include the copy number of any complete SMN gene determined to be the total, taking into account (a) the number of sequence reads of the plurality of sequence reads having bases supporting any one of the bases of the plurality of SMN1 gene-specific bases and (b) the number of sequence reads of the plurality of sequence reads having bases supporting any one of the bases of the plurality of corresponding SMN2 gene-specific bases. The method may include: determining the possible copy number of the possible combination to be the consistent possible copy number of the SMN1 gene.

[0021] In some embodiments, determining the copy number of the SMN1 gene includes: determining that the copy number of the SMN1 gene is zero, one, or more than one. In some embodiments, the method includes: determining the spinal muscular atrophy (SMA) status of the subject based on the copy number of the SMN1 gene. The SMA status of the subject may include SMA, SMA carrier but not SMA, and not SMA carrier. In some embodiments, the method includes: using the number of sequence reads of the plurality of sequence reads aligned with g.27134 of the SMN1 gene and the bases of the sequence reads aligned with g.27134 of the SMN1 gene to determine that the subject is a silent SMA carrier.

[0022] In some embodiments, the method includes: determining a treatment recommendation for the subject based on the determined copy number of the SMN1 gene. The treatment recommendation may include administering Nusinersen and / or Zolgensma to the subject.

[0023] Disclosed herein are methods for genotyping the cytochrome P450 family 2 subfamily D member 6 (CYP2D6) gene. In some embodiments, the method for genotyping the CYP2D6 gene is under the control of a processor, such as a hardware processor or a virtual processor, and includes: receiving sequence data, the sequence data including a plurality of sequence reads obtained from a sample of a subject and aligned to the CYP2D6 gene or the cytochrome P450 family 2 subfamily D member 7 (CYP2D7) gene. The method may include: determining (i) a first number of sequence reads of the plurality of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene. The method may include: respectively using (i) the length of the CYP2D6 gene or the CYP2D7 gene to determine (i) a first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene. The method may include: using a Gaussian mixture model comprising a plurality of Gaussian functions each representing a different integer copy number to determine (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene, taking into account (i) the first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene. The method may include: for one of a plurality of CYP2D6 gene-specific bases, determining the most likely combination among a plurality of possible combinations of the possible copy number of the CYP2D6 gene and the possible copy number of the CYP2D7 gene, each including a total of the determined total copy number of the CYP2D6 gene and the CYP2D7 gene, taking into account (a) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D6 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D7 gene-specific base corresponding to the CYP2D6 gene-specific base. The method may include: using the most likely combination of the possible copy number of the CYP2D6 gene and the possible copy number of the CYP2D7 gene determined for the CYP2D6 gene-specific base to determine the alleles of the CYP2D6 gene that the subject has.

[0024] In some embodiments, the sequencing data includes whole genome sequencing (WGS) data or short-read WGS data. The subject may be a neonatal subject, a pediatric subject, an adolescent subject, or an adult subject. The sample may comprise cellular or cell-free DNA. The sample may comprise cellular or cell-free DNA.

[0025] In some embodiments, sequence reads of the plurality of sequence reads are aligned to the CYP2D6 gene or the CYP2D7 gene, wherein the alignment quality score is approximately zero. In some embodiments, determining (i) a first number of sequence reads of the plurality of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene comprises: determining (i) a first number of sequence reads of the plurality of sequence reads aligned to at least one exon or intron of the CYP2D6 gene or an exon or intron of the CYP2D7 gene.

[0026] In some embodiments, determining (i) a first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene comprises: using the length of (i) the CYP2D6 gene or the CYP2D7 gene, respectively, to determine (i) a first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, and determining (iii) the depth of sequence reads of a region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene in the sequence data. Determining (i) a first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene and (ii) a second normalized number of sequence reads aligned to a second region may comprise: using the length of (i) the CYP2D6 gene or the CYP2D7 gene, respectively, to determine (i) a first CYP2D6 gene or CYP2D7 gene length-normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene. Determining (i) a first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene and (ii) a second normalized number of sequence reads aligned to a second region may comprise: using the depth of sequence reads of a region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene to determine (i) a first normalized depth of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene according to the (i) CYP2D6 gene or CYP2D7 gene length-normalized number, and the first normalized depth of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene is the first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, respectively.

[0027] In some embodiments, determining a first normalized quantity of sequence reads aligned to (i) the CYP2D6 gene or the CYP2D7 gene includes: using (i) the GC content of the CYP2D6 gene or the CYP2D7 gene to determine the first normalized quantity of sequence reads aligned to (i) the CYP2D6 gene or the CYP2D7 gene, and determining (iii) the depth of sequence reads of a region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene in the sequence data, and (iv) determining the GC content of the region of the genome. The depth of the region may include the average depth or the median depth of sequence reads of a region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene in the sequencing data. The region may contain approximately 3000 preselected regions each of length approximately 2 kb and spanning the subject's genome. In some embodiments, the first normalized quantity of sequence reads aligned to (i) the CYP2D6 gene or the CYP2D7 gene and / or the second normalized quantity of sequence reads aligned to a second region is from about 30 to about 40.

[0028] In some embodiments, the Gaussian mixture model includes a one-dimensional Gaussian mixture model. The plurality of Gaussian functions of the Gaussian mixture model may represent integer copy numbers from 0 to 10. The mean of each of the plurality of Gaussian functions may be the integer copy number represented by the Gaussian function.

[0029] In some embodiments, determining the total copy number of (i) the CYP2D6 gene and the CYP2D7 gene includes: considering the first normalized quantity of sequence reads aligned to (i) the CYP2D6 gene or the CYP2D7 gene, using a Gaussian mixture model and a first predetermined posterior probability threshold to determine the total copy number of (i) the CYP2D6 gene and the CYP2D7 gene. The first predetermined posterior probability threshold may be 0.95.

[0030] In some embodiments, considering (a) the number of sequence reads of the plurality of sequence reads having bases supporting CYP2D6 gene-specific bases and (b) the number of sequence reads of the plurality of sequence reads having bases supporting corresponding CYP2D7 gene-specific bases, the most likely combination of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene is associated with the highest posterior probability relative to other combinations in the plurality of combinations.

[0031] In some embodiments, determining the most likely combination of the possible copy number of the CYP2D6 gene and the possible copy number of the CYP2D7 gene includes: determining the most likely combination among the multiple possible combinations of the possible copy number of the CYP2D6 gene and the possible copy number of the CYP2D7 gene, each including a total copy number equal to the determined total copy number of the CYP2D6 gene and the CYP2D7 gene, taking into account the ratio of (a) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D6 gene-specific bases to (b) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D7 gene-specific bases corresponding to the CYP2D6 gene-specific bases. Determining the most likely combination of the possible copy number of the CYP2D6 gene and the possible copy number may include: determining (a) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D6 gene-specific bases and (b) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D7 gene-specific bases corresponding to the CYP2D6 gene-specific bases; determining the ratio of (a) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D6 gene-specific bases to (b) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D7 gene-specific bases corresponding to the CYP2D6 gene-specific bases; and determining the most likely combination among the multiple possible combinations of the possible copy number of the CYP2D6 gene and the possible copy number of the CYP2D7 gene, each including a total copy number equal to the determined total copy number of the CYP2D6 gene and the CYP2D7 gene, taking into account the ratio of (a) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D6 gene-specific bases to (b) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D7 gene-specific bases corresponding to the CYP2D6 gene-specific bases.

[0032] In some embodiments, determining the alleles of the CYP2D6 gene that a subject has includes: determining one or more structural variants of the CYP2D6 gene that the subject has using the most likely combination of the likely copy number of the CYP2D6 gene and the likely copy number of the CYP2D7 gene determined for CYP2D6 gene-specific bases. In some embodiments, determining the most likely combination of the likely copy number of the CYP2D6 gene and the likely copy number of the CYP2D7 gene includes: for each of the plurality of CYP2D6 gene-specific bases, determining the likely copy number of the CYP2D6 gene and the likely copy number of the CYP2D7 gene in each of a plurality of possible combinations that together total the determined total copy number of the CYP2D6 gene and the CYP2D7 gene, taking into account (a) the number of sequence reads of the plurality of sequence reads having bases that support the CYP2D6 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads having bases that support the CYP2D7 gene-specific base corresponding to the CYP2D6 gene-specific base, and determining the most likely combination associated with the highest posterior probability. Determining the one or more structural variants of the CYP2D6 gene that the subject has may include: using the most likely combination of the likely copy number of the CYP2D6 gene and the likely copy number of the CYP2D7 gene determined for each of the plurality of CYP2D6 gene-specific bases to determine the one or more structural variants. In some embodiments, determining the one or more structural variants of the CYP2D6 gene that the subject has includes: determining the one or more structural variants of the CYP2D6 gene that the subject has based on the likely copy number of the CYP2D6 gene in the most likely combination determined for two or more different ones of the plurality of CYP2D6 gene-specific bases and the positions of the two or more CYP2D6 gene-specific bases.

[0033] In some embodiments, the CYP2D6 gene-specific bases are identical to each of the plurality of CYP2D6 gene-specific bases other than the CYP2D6 gene-specific bases that exceed a pre-determined identity threshold. The identity threshold can be 97%. The plurality of CYP2D6 gene-specific bases can include 118 CYP2D6 gene-specific bases. The plurality of CYP2D6 gene-specific bases can be different if the subject is of a first race, can be different if the subject is of a second race, and can be different if the subject is of an unknown race. The race of the subject may be unknown, and the plurality of CYP2D6 gene-specific bases may not be race-specific. The race of the subject may be known, and the plurality of CYP2D6 gene-specific bases may be specific to the race of the subject. In some embodiments, the method includes: receiving race information of the subject. The method can include: selecting the plurality of CYP2D6 gene-specific bases from the plurality of CYP2D6 gene-specific bases based on the received race information.

[0034] In some embodiments, the method includes: determining a second number of sequence reads of the plurality of sequence reads that align with a spacer region between the CYP2D7 gene and the repeat element REP7 downstream of the CYP2D7 gene. The method can include: using the length of the spacer region to determine a second normalized number of sequence reads that align with the spacer region. The method can include: using a Gaussian mixture model to determine the copy number of the spacer region in view of the second normalized number of sequence reads that align with the spacer region. Determining the one or more structural variants of the CYP2D6 gene that the subject has can include: using the most likely combination of the possible copy number of the CYP2D6 gene, the possible copy number of the CYP2D7 gene, and the copy number of the spacer region determined for the CYP2D6 gene-specific bases to determine the one or more structural variants of the CYP2D6 gene that the subject has. The one or more structural variants can include a CYP2D6 / CYP2D7 fusion allele having a spacer region and the repeat element REP7 downstream of the CYP2D6 / CYP2D7 fusion allele.

[0035] In some embodiments, the method includes: determining, using the received sequence data, one or more minor variants of the CYP2D6 gene that a subject has. In some embodiments, determining the one or more minor variants of the CYP2D6 gene that a subject has includes: for a minor variant position of the CYP2D6 gene that is associated with a minor variant allele of the CYP2D6 gene, determining, taking into account (a) the number of sequence reads having a base that supports the minor variant allele of the CYP2D6 gene at the minor variant position and (b) the number of sequence reads having a base that supports the reference allele of the CYP2D6 gene at the minor variant position, the most likely combination of the possible copy number of the minor variant allele of the CYP2D6 gene at the minor variant position and the possible copy number of the reference allele of the CYP2D6 gene at the minor variant position that totals the copy number of the CYP2D6 gene at the minor variant position, wherein the possible copy number of the minor variant allele of the CYP2D6 gene in the most likely combination at the minor variant position indicates the one or more minor variants of the CYP2D6 gene. In some embodiments, determining the one or more minor variants of the CYP2D6 gene that a subject has includes: for each minor variant position among a plurality of minor variant positions of the CYP2D6 gene, the minor variant position being associated with a minor variant allele of the CYP2D6 gene, determining, taking into account (a) the number of sequence reads having a base that supports the minor variant allele of the CYP2D6 gene at the minor variant position and (b) the number of sequence reads having a base that supports the reference allele of the CYP2D6 gene at the minor variant position, the most likely combination of the possible copy number of the minor variant allele of the CYP2D6 gene at the minor variant position and the possible copy number of the reference allele of the CYP2D6 gene at the minor variant position that totals the copy number of the CYP2D6 gene at the minor variant position, wherein the possible copy number of the minor variant alleles of the CYP2D6 gene in the most likely combination at the plurality of minor variant positions indicates the one or more minor variants of the CYP2D6 gene.

[0036] In some embodiments, the method includes: for a variant position of the CYP2D6 gene that is associated with a minor variant allele of the CYP2D6 gene, determining the most likely combination of the possible copy number of the minor variant allele of the CYP2D6 gene at the variant position and the possible copy number of the reference allele of the CYP2D6 gene at the variant position that together total the copy number of the CYP2D6 gene at the variant position, taking into account (a) the number of sequence reads that align with the CYP2D6 gene and overlap the variant position and have bases supporting the minor variant allele of the CYP2D6 gene at the variant position and (b) the number of sequence reads that align with the CYP2D6 gene and overlap the variant position and have bases supporting the reference allele of the CYP2D6 gene at the variant position; and using the possible copy number of the minor variant allele of the CYP2D6 gene in the determined most likely combination to determine one or more variants of the CYP2D6 gene. In some embodiments, the method includes: for each variant position among a plurality of variant positions of the CYP2D6 gene, where the variant position is associated with a minor variant allele of the CYP2D6 gene, determining the most likely combination of the possible copy number of the minor variant allele of the CYP2D6 gene at the variant position and the possible copy number of the reference allele of the CYP2D6 gene at the variant position that together total the copy number of the CYP2D6 gene at the variant position, taking into account (a) the number of sequence reads that align with the CYP2D6 gene and overlap the variant position and have bases supporting the minor variant allele of the CYP2D6 gene at the variant position and (b) the number of sequence reads that align with the CYP2D6 gene and overlap the variant position and have bases supporting the reference allele of the CYP2D6 gene at the variant position; and using the possible copy number of the minor variant allele of the CYP2D6 gene in the determined most likely combination at the plurality of variant positions to determine one or more variants of the CYP2D6 gene.

[0037] In some embodiments, the minor variant position is in the CYP2D6 / CYP2D7 homology region, and determining the most likely combination involves determining the most likely combination of the possible copy numbers of the minor variant alleles of the CYP2D6 gene at the minor variant position and the possible copy numbers of the reference alleles of the CYP2D6 gene at the minor variant position that together total the copy number of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads of bases supporting the minor variant allele of the CYP2D6 gene aligned to the CYP2D6 gene or the CYP2D7 gene and / or (b) the number of sequence reads of bases supporting the reference allele of the CYP2D6 gene at the minor variant position aligned to the CYP2D6 gene or the CYP2D7 gene. In some embodiments, the minor variant position is not in the CYP2D6 / CYP2D7 homology region, and determining the most likely combination involves determining the most likely combination of the possible copy numbers of the minor variant alleles of the CYP2D6 gene at the minor variant position and the possible copy numbers of the reference alleles of the CYP2D6 gene at the minor variant position that together total the copy number of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads of bases supporting the minor variant allele of the CYP2D6 gene aligned to the CYP2D6 gene and not to the CYP2D7 gene and / or (b) the number of sequence reads of bases supporting the reference allele of the CYP2D6 gene at the minor variant position aligned to the CYP2D6 gene and not to the CYP2D7 gene.

[0038] In some embodiments, the method includes determining the copy number of the CYP2D6 gene at the minor variant position. The copy number of the CYP2D6 gene at the minor variant position may include the copy number of the CYP2D6 gene. The copy number of the CYP2D6 gene at the minor variant position may include the copy number of the CYP2D6 gene of the determined most likely combination of the possible copy numbers of the CYP2D6 gene. The copy number of the CYP2D6 gene at the minor variant position may include the copy number of the CYP2D6 gene of the determined most likely combination and closest to the minor variant position of the possible copy numbers of the CYP2D6 gene. The copy number of the CYP2D6 gene at the minor variant position may include the copy number of the CYP2D6 gene at the 5' position or the 3' position of the minor variant position. In some embodiments, the method includes: (a) determining the number of sequence reads of bases supporting the minor variant allele of the CYP2D6 gene; and (b) determining the number of sequence reads of bases supporting the reference allele of the CYP2D6 gene.

[0039] In some embodiments, determining the alleles of the CYP2D6 gene that a subject has includes: determining the alleles of the CYP2D6 gene that the subject has (e.g., 2, 3, 4, 5 or more alleles). In some embodiments, determining the alleles of the CYP2D6 gene that a subject has includes: using the one or more structural variants of the determined CYP2D6 gene and / or the one or more minor variants of the determined CYP2D6 gene to determine the star alleles and / or haplotypes of the CYP2D6 gene that the subject has, optionally where the star alleles are associated with known functions.

[0040] In some embodiments, the method includes: using the determined alleles of the CYP2D6 gene to determine the level of CYP2D6 enzyme activity in a subject. The enzyme activity can be poor, moderate, normal or ultrametabolic. In some embodiments, the method includes determining a treatment dose recommendation and / or a treatment recommendation for a subject based on the alleles of the CYP2D6 gene that the subject has.

[0041] Disclosed herein are systems for paralog genotyping. In some embodiments, a system for paralog genotyping includes: a non-transitory memory configured to store executable instructions and sequence data, the sequence data including a plurality of sequence reads obtained from a sample of a subject and aligned to a first paralog or a second paralog. The system can include: a processor (such as a hardware processor or a virtual processor) in communication with the non-transitory memory, the processor programmed by the executable instructions to perform: using a Gaussian mixture model that includes a plurality of Gaussian functions each representing a different integer copy number to determine the copy number of a first type of paralog, taking into account (i) a first number of sequence reads aligned to a first region. The hardware processor is programmed by the executable instructions to perform: for one base among a plurality of first paralog-specific bases, determining the most likely combination among a plurality of possible combinations of the possible copy numbers of a first type of first paralog and the possible copy numbers of a first type of second paralog, each including a total of the determined copy number of the first type of paralog, taking into account (a) the number of sequence reads of the plurality of sequence reads having a base that supports the first paralog-specific base and (b) the number of sequence reads of the plurality of sequence reads having a base that supports a second paralog-specific base corresponding to the first paralog-specific base of the second paralog. The hardware processor is programmed by the executable instructions to perform: using the most likely combination of the possible copy numbers of the first paralog and the possible copy numbers of the second paralog determined for the first paralog-specific base to determine the copy number or allele of the first paralog. In some embodiments, the first paralog and the second paralog have at least 90% sequence identity.

[0042] In some embodiments, the hardware processor is programmed by executable instructions to perform: determining a first quantity of sequence reads of a plurality of sequence reads obtained from a subject in sequence data and aligned to a first region. The method may include: using the length of the first region to determine a first normalized quantity of sequence reads aligned to the first region, wherein determining the copy number of a first type of paralog includes: considering the first normalized quantity of sequence reads aligned to the first region and using a Gaussian mixture model to determine the copy number of the first type of paralog. The hardware processor may be programmed by executable instructions to perform: may include: receiving sequence data including the plurality of sequence reads aligned to the first region.

[0043] In some embodiments, the hardware processor is programmed by executable instructions to perform: considering a second quantity of sequence reads aligned to a second region and using Gaussian mixture to determine the copy number of one or more paralogs of a second type. Determining the copy number or allele of a first paralog may include: using the most likely combination of the possible copy numbers of the first paralog and the possible copy numbers of a second paralog determined for first paralog-specific bases and the copy number of the one or more paralogs of the second type to determine the copy number or allele of the first paralog. The method may include: determining the copy number of a third type of paralog from the copy number of the first type of paralog and the copy number of the second type of paralog. Determining the copy number or allele of a first paralog may include: using the most likely combination of the possible copy numbers of the first paralog and the possible copy numbers of a second paralog determined for first paralog-specific bases to determine the copy number or allele of the first paralog.

[0044] In some embodiments, the first paralog is the survival motor neuron 1 (SMN1) gene. The second paralog may be the survival motor neuron 2 (SMN2) gene. The first region may include at least exons 1 to 6 of the SMN1 gene and at least exons 1 to 6 of the SMN2 gene. The second region may include at least one of exons 7 and 8 of the SMN1 gene and at least one of exons 7 and 8 of the SMN2 gene. The first type of paralog may include the full-length SMN1 gene and the full-length SMN2 gene. The one or more paralogs of the second type may include the full-length SMN1 gene, the full-length SMN2 gene, a truncated SMN1 gene, or a truncated SMN2 gene. The copy number of the first paralog may include the copy number of the SMN1 gene.

[0045] In some embodiments, the first paralog is the cytochrome P450 family 2 subfamily D member 6 (CYP2D6) gene. The second paralog can be the cytochrome P450 family 2 subfamily D member 7 (CYP2D7) gene. The first region can include the CYP2D6 gene and the CYP2D7 gene. The second region can include the spacer region between the CYP2D7 gene and the repetitive element REP7 downstream of the CYP2D7 gene. The first type of paralog can include the CYP2D6 gene and the CYP2D7 gene. The second type of the one or more paralogs can include the CYP2D6 / CYP2D7 fusion allele with the spacer region and the repetitive element REP7 downstream of the CYP2D6 / CYP2D7 fusion allele. The copy number of the first paralog can include the allele of the CYP2D6 gene that the subject has, which is a minor variant or a structural variant of the CYP2D6 gene.

[0046] Embodiments disclosed herein include a system (e.g., a computing system) that includes a non-transitory memory configured to store executable instructions; and a processor (e.g., a hardware processor or a virtual processor) communicatively coupled to the non-transitory memory, the hardware processor programmed by the executable instructions to perform any of the methods disclosed herein. Embodiments disclosed herein include a device (e.g., an electronic device) that includes a non-transitory memory configured to store executable instructions; and a processor (e.g., a hardware processor or a virtual processor) communicatively coupled to the non-transitory memory, the hardware processor programmed by the executable instructions to perform any of the methods disclosed herein. Embodiments disclosed herein include a computer-readable medium that includes executable instructions that, when executed by a processor (e.g., a hardware processor or a virtual processor) of a system or a device, cause the hardware processor to perform any of the methods disclosed herein.

[0047] Details of one or more specific implementations of the subject matter described in this specification are set forth in the accompanying drawings and the following description. Other features, aspects, and advantages will become apparent from the description, the drawings, and the claims. Neither the summary nor the following detailed description is intended to limit or restrict the scope of the subject matter of the invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figures 1A to 1E Shows the reasons for SMA and SMN copy number calls according to one embodiment of the method disclosed herein.

[0049] Figures 2A to 2C Shows the population distribution of SMN1 / 2 copy numbers determined using one embodiment of the method disclosed herein.

[0050] Figure 3Shows the SMA identified in two trios of the next-generation children's project and verified using MLPA.

[0051] Figure 4 Shows that the population frequencies determined using an embodiment of the method disclosed herein are consistent with previous studies.

[0052] Figure 5 Non-limiting exemplary IGV snapshots showing that CYP2D6 is highly polymorphic and located downstream of CYP2D7 (a pseudogene paralog of CYP2D6).

[0053] Figure 6 Non-limiting exemplary schematic diagrams of CYP2D6 / 7 gene deletions, duplications, and fusion genes.

[0054] Figure 7 Non-limiting exemplary curves showing that the allele frequencies determined by the method are consistent with the PharmVar database from the Pharmacogene Variation (PharmVar) Consortium.

[0055] Figure 8 Flowchart of an exemplary method for showing the determination of the copy number of the survival motor neuron 1 (SMN1) gene using sequencing data.

[0056] Figure 9 Flowchart of an exemplary method for genotyping the cytochrome P450 family 2 subfamily D member 6 (CYP2D6) gene using sequencing data.

[0057] Figure 10 Flowchart of an exemplary method for genotyping paralogous genes using sequencing data.

[0058] Figure 11 Block diagram of an exemplary computing system configured to perform paralogous gene genotyping using sequencing data.

[0059] Figure 12A and Figure 12B Shows non-limiting exemplary curves illustrating common CNVs affecting the SMN1 / SMN2 locus. Figure 12AShows the depth profile across the SMN1 / SMN2 region. Samples with 2, 3, 4, and 5 total SMN1+SMN2 copy numbers are shown as dots respectively. For each CN category, the depths of 50 samples are summed. Each point represents the normalized depth value in a 100bp window. Read counts are calculated in each 100bp window, the reads for both SMN1 and SMN2 are summed, and normalized to the depth of the wild-type sample (CN = 4). SMN exons are represented as purple boxes. The two x-axes show the coordinates in SMN1 (bottom) and SMN2 (top). Figure 12B Shows the depth profile aggregated from 50 samples carrying deletions of exons 7 and 8, shown as dots. Read depths are calculated in the same way as Figure 12A described above.

[0060] Figure 13 Shows a non-limiting exemplary scatter plot of total SMN (SMN1+SMN2) copy number (x-axis, called from read depth in exons 1 to 6) and full SMN copy number (y-axis, called from read depth in exons 7 to 8).

[0061] Figures 14A to 14D Shows the distribution of SMN1 / SMN2 / SMN* copy numbers in the population. Figure 14A Is a non-limiting exemplary figure that shows the percentage of samples showing CN call concordance with c.840C>T at 16 SMN1-SMN2 base difference sites in African and non-African population groups. Site 13* is the c.840C>T splice variant site. The black horizontal line represents 85% concordance. Figure 14B Shows a non-limiting exemplary bar chart of SMN1, SMN2, and SMN* copy number distributions in five populations in the 1kGP and NIHR BioResource cohorts (values are shown in Table 15). Figure 14C Is a non-limiting exemplary curve graph of SMN1 CN versus total SMN2 CN (full SMN2+SMN*). Figure 14D Shows two trios where the SMA probands were detected by the caller and orthogonally confirmed in the NIHR BioResource cohort. The CN of each allele of SMN1, SMN2, and SMN* was phased and each member of the trio was labeled.

[0062] Figure 15 Shows non-limiting exemplary curve graphs, each of which shows the posterior probability distribution of simulated SMN1 CN using individual sites and SMN1:SMN2 CN combinations at different read depths.

[0063] Figure 16Shows a non - restrictive exemplary IGV snapshot of the SMN2 region in a sample with exon 7 - 8 deletion. The horizontal line connects two reads in pairs in the center - aligned track. The BLAT results of two split reads spanning the breakpoint are shown in the bottom track, which shows two fragments of the same read aligned to either side of the deletion breakpoint.

[0064] Figure 17 Shows a non - restrictive exemplary graph that shows the correlation between the raw SMN1 CN at a 15 - base difference near c840.C>T and the raw SMN1 CN at the c840.C>T locus. The raw SMN1 CN at each locus is calculated as the CN of the full SMN multiplied by the fraction of the read count supporting SMN1 in the read count supporting SMN1 + SMN2. The correlation coefficient is listed in the title of each graph.

[0065] Figure 18A and Figure 18B Shows a non - restrictive exemplary graph that shows the SMN1 / SMN2 haplotypes in samples with SMN1:2SMN2:0 and SMN1:2SMN2:1 in 1kGP. The y - axis shows the raw SMN1 CN as defined in Figure 16 The x - axis shows 16 loci, which are listed and explained in Table 8. Index #13 represents the c840.C>T locus. Samples with SMN1:2SMN2:0 are shown together in the upper - left panel. Samples with SMN1:2SMN2:1 are shown as 5 clusters. Figure 18A . Non - African. Figure 18B . African.

[0066] Figure 19 Shows a non - restrictive exemplary IGV snapshot showing a 1.9 - kb deletion of SMN1 in MB509.

[0067] Figure 20 Shows a non - restrictive exemplary graph that shows SMN1 / SMN2 / SMN*CN in 1kGP and the NIHR cohort.

[0068] Figure 21A and Figure 21B Shows the differences and no - calls in the validation samples.

[0069] Figure 22 Shows the CN calls derived from BWA and Isaac BAM.

[0070] Figure 23Non-limiting exemplary curves showing the quality of WGS data in the CYP2D6 / 7 region. The average mapping quality of 1kGP samples was plotted for each position in the CYP2D6 / 7 region. A median filter was applied in a 200bp window. Nine exons of REP6, REP7, and CYP2D6 / 7 are boxed in the left box (CYP2D6) and right box (CYP2D7). The two 2.8kb repeat regions downstream of CYP2D6 (REP6) and CYP2D7 (REP7) are identical and largely non-alignable. The dashed box indicates the spacer region between CYP2D7 and REP7. Two major homology regions within the gene are shaded.

[0071] Figure 24 Structural variants verified by PacBio CCS reads are shown. PacBio reads support deletions (*5), duplications, and fusions (*36, *68, and *13). Curves were generated using sv-viz2 (zotero.org / google-docs / ?xAunA6). For deletions and duplications, due to significant homology in the REP region, the exact position of the breakpoints within the REP is not available. The breakpoints in A and B are for illustrative purposes only.

[0072] Figure 25 Non-limiting exemplary curves showing the frequencies of CYP2D6 alleles in five ethnic groups for the ten most common haplotypes with altered CYP2D6 function. One haplotype (*2x2) has enhanced function, two haplotypes (*4 and *4+*68) are non-functional, and the remaining haplotypes have reduced function.

[0073] Figure 26 Shows that the CYP2D6 / CYP2D7 base difference sites have high variability in the population. The Y-axis shows the sample frequency in which the CN of the CYP2D6 base is called 2 in all samples with a total CYP2D6 + CYP2D7 CN of 4. The X-axis shows the genomic coordinates in hg38. CYP2D6 exons are plotted as gray boxes above the figure. The black horizontal line represents the 98% cutoff.

[0074] Figure 27 Shows the original CYP2D6 CN across the CYP2D6 / 7 differentiation site in an example with an SV. The original CYP2D6 CN was calculated as the total CYP2D6 + CYP2D7 CN multiplied by the ratio of CYP2D6 supporting reads to the total number of supporting reads for CYP2D6 and CYP2D7. The large diamond represents the copy number of the CYP2D6-derived gene (which can be a full CYP2D6 or a fusion gene ending with CYP2D6) at the end of the gene, calculated as the total CN of CYP2D6 + CYP2D7 minus the CN of the CYP2D7 spacer region (seeFigure 23 ) To detect SV, CYP2D6 CN was called at each locus, and a change in CYP2D6 CN within the gene indicated the presence of SV. For example, in HG01161, the CYP2D6 CN changed from 2 to 1 between exon 7 and exon 9, indicating a CYP2D7-CYP2D6 hybrid gene. In HG00553, the CYP2D6 CN changed from 2 to 3 between exon 1 and exon 2, indicating a CYP2D6-CYP2D7 hybrid gene.

[0075] Figure 28 PacBio data confirmed the *10D fusion in HG00421 is shown. A sample with *36 (HG00612) is shown in comparison. PacBio reads containing the fusion are those with shaded bases, which represent soft clips prepared by the aligner and derived from the CYP2D7 portion of the fusion. The fusion breakpoints are close to each other, but the breakpoint of *36 is upstream of the base differences in exon 9 (those within the black box), while the breakpoint of *10D is downstream, thus leaving the CYP2D6 gene intact.

[0076] Figure 29 PacBio data showed a false *61 (CYP2D6 / CYP2D7 hybrid) call by Aldy in HG02622. The expected genotype was *17 / *45, but Aldy called *61-like / *78 (both *61 and *78 are star alleles with SV). PacBio data showed no structural variants in this region (each read aligned completely, with no soft clips indicating unaligned parts).

[0077] Figure 30A and Figure 30B A novel *10+*36+*36+*83 haplotype in HG00597 is shown. Figure 30A . The depth profile is as Figure 27 shown, which shows that HG00597 has three copies of the *36-like fusion, all of which have breakpoints in the homologous region between exon 7 and exon 9. Figure 30B . An IGV screenshot of the PacBio data, which shows all reads containing the fusion, i.e., those aligned with soft clips. One copy of the fusion gene does not have g.42130692G>A, an SNP that is in *36 but not in *83, as shown in the region flanked by two black vertical lines. This copy is *83 and is different from the copy reported in PharmVar, which is a fusion gene with REP7 instead of REP6, otherwise the copy number in the region downstream of exon 9 would be 3 instead of Figure 30A 2 in

[0078] Figure 31A and Figure 31B Comparison was made between 1kGP and pharmGKB frequencies. Each point represents a haplotype with a frequency greater than or equal to 0.5% in 1kGP or pharmGKB. SV-related haplotypes are marked, including the two haplotypes with the largest deviations (*10+*36 in East Asians, *4+*68 in Europeans). Other haplotypes with deviation values are annotated (*2, *41, *34, *39, *2, and *29). A diagonal line is drawn for each subfigure. The correlation coefficients for each population are listed (*10+*36 was excluded in East Asians and *4+*68 was excluded in Europeans for calculation). Figure 31B Values in the low value range (<5%) are shown.

[0079] Figure 32 A non-limiting exemplary IGV snapshot showing de novo assembly of PacBio reads in HG00733 excluding the *68 fusion.

[0080] Throughout the figures, reference numerals may be repeated to indicate corresponding relationships between reference elements. The figures are provided to illustrate exemplary embodiments described herein and are not intended to limit the scope of the disclosure. Detailed Description

[0081] In the following detailed description, reference is made to the accompanying figures, which form a part of the detailed description. In the figures, like symbols generally identify like components unless the context dictates otherwise. The exemplary embodiments described in the detailed description, figures, 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 will be readily understood that the aspects of the present disclosure, as generally described and illustrated in the figures herein, may be arranged, substituted, combined, separated, and designed in a variety of different configurations, all of which are expressly contemplated herein and form a part of the present disclosure.

[0082] All patents, published patent applications, other publications, and sequences from GenBank, as well as other databases mentioned herein, are incorporated by reference in their entirety relative to the relevant art.

[0083] Methods disclosed herein include methods for determining the copy number of the survival motor neuron 1 (SMN1) gene and / or the survival motor neuron 2 (SMN2) gene. In some embodiments, a method for determining the copy number of the SMN1 gene and / or the SMN2 gene is under the control of a processor (such as a hardware processor or a virtual processor) and includes: receiving sequence data that includes a plurality of sequence reads obtained from a sample of a subject and aligned to the SMN1 gene or the SMN2 gene. The method can include: determining (i) a first number of sequence reads of the plurality of sequence reads that are aligned to a first SMN1 or SMN2 region that respectively includes at least one of exons 1 to 6 of the SMN1 gene or the SMN2 gene and (ii) a second number of sequence reads of the plurality of sequence reads that are aligned to a second SMN1 or SMN2 region that respectively includes at least one of exons 7 and 8 of the SMN1 gene or the SMN2 gene. The method can include: using (i) the length of the first SMN1 or SMN2 region and (ii) the length of the second SMN1 or SMN2 region to respectively determine (i) a first normalized number of sequence reads that are aligned to the first SMN1 or SMN2 region and (ii) a second normalized number of sequence reads that are aligned to the second SMN1 or SMN2 region. The method can include: using a Gaussian mixture model that includes a plurality of Gaussian functions each representing a different integer copy number to determine (i) the copy number of the total survival motor neuron (SMN) gene that is respectively a full-length SMN1 gene, a full-length SMN2 gene, a truncated SMN1 gene, or a truncated SMN2 gene and (ii) the copy number of any full-length SMN gene that is respectively a full-length SMN1 gene or a full-length SMN2 gene, taking into account (i) the first normalized number of sequence reads that are aligned to the first SMN1 or SMN2 region and (ii) the second normalized number of sequence reads that are aligned to the second SMN1 or SMN2 region. The method can include: for a base among a plurality of SMN1 gene-specific bases associated with the full-length SMN1 gene, determining the most likely combination among a plurality of possible combinations of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene that respectively include a total number of copies of any full-length SMN gene determined, taking into account (a) the number of sequence reads of the plurality of sequence reads that have a base supporting the SMN1 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads that have a base supporting the SMN2 gene-specific base corresponding to the SMN1 gene-specific base of the SMN2 gene. The method can include: using the most likely combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene determined for the SMN1 gene-specific base to determine the copy number of the SMN1 gene and / or the SMN2 gene.

[0084] Disclosed herein are methods for genotyping the cytochrome P450 family 2 subfamily D member 6 (CYP2D6) gene. In some embodiments, the method for genotyping the CYP2D6 gene is under the control of a processor, such as a hardware processor or a virtual processor, and includes: receiving sequence data including a plurality of sequence reads obtained from a sample of a subject and aligned to the CYP2D6 gene or the cytochrome P450 family 2 subfamily D member 7 (CYP2D7) gene. The method may include: determining (i) a first number of sequence reads of the plurality of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene. The method may include: using respectively (i) the length of the CYP2D6 gene or the CYP2D7 gene to determine (i) a first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene. The method may include: using a Gaussian mixture model comprising a plurality of Gaussian functions each representing a different integer copy number to determine (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene, taking into account (i) the first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene. The method may include: for one of a plurality of CYP2D6 gene-specific bases, determining the most likely combination among a plurality of possible combinations of the possible copy number of the CYP2D6 gene and the possible copy number of the CYP2D7 gene, each including a total of the determined total copy number of the CYP2D6 gene and the CYP2D7 gene, taking into account (a) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D6 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D7 gene-specific base corresponding to the CYP2D6 gene-specific base. The method may include: using the most likely combination of the possible copy number of the CYP2D6 gene and the possible copy number of the CYP2D7 gene determined for the CYP2D6 gene-specific base to determine the alleles of the CYP2D6 gene that the subject has.

[0085] Disclosed herein are methods for paralog genotyping. In some embodiments, the methods for paralog genotyping are under the control of a processor (such as a hardware processor or a virtual processor) and include: receiving sequence data including a plurality of sequence reads obtained from a sample of a subject and aligned to a first paralog or a second paralog. The method may include: using a Gaussian mixture model including a plurality of Gaussian functions each representing a different integer copy number to determine the copy number of a first type of paralog, considering (i) a first number of sequence reads aligned to a first region. The method may include: for a base among a plurality of first paralog-specific bases, determining the most likely combination among a plurality of possible combinations each including a possible copy number of a first type of first paralog and a possible copy number of a first type of second paralog that together total the determined copy number of the first type of paralog, considering (a) the number of sequence reads among the plurality of sequence reads having a base supporting the first paralog-specific base and (b) the number of sequence reads among the plurality of sequence reads having a base supporting a second paralog-specific base corresponding to the first paralog-specific base of the second paralog. The method may include: using the most likely combination of the possible copy numbers of the first paralog and the second paralog determined for the first paralog-specific base to determine the copy number or allele of the first paralog.

[0086] Embodiments disclosed herein include a system (e.g., a computing system) including a non-transitory memory configured to store executable instructions; and a processor (e.g., a hardware processor or a virtual processor) in communication with the non-transitory memory, the hardware processor programmed by the executable instructions to perform any method disclosed herein. Embodiments disclosed herein include a device (e.g., an electronic device) including a non-transitory memory configured to store executable instructions; and a processor (e.g., a hardware processor or a virtual processor) in communication with the non-transitory memory, the hardware processor programmed by the executable instructions to perform any method disclosed herein. Embodiments disclosed herein include a computer-readable medium including executable instructions that, when executed by a processor (e.g., a hardware processor or a virtual processor) of a system or device, cause the hardware processor to perform any method disclosed herein.

[0087] Spinal muscular atrophy diagnosis and carrier screening from whole-genome sequencing data

[0088] Spinal muscular atrophy (SMA) is characterized by voluntary muscle weakness and is the leading genetic cause of early childhood death, with an incidence of 1 in 6,000 to 10,000 live births and a carrier frequency of 1:40 to 1:801,2. SMA is caused by mutations in the SMN1 (survival motor neuron 1) gene ( Figure 1A ). The duplicated gene SMN2 differs from SMN1 by only a few base pairs, and one of them (the c.840C>T splicing variant in exon 7) has a functional outcome. By disrupting a splicing enhancer, the c.840C>T mutation results in increased skipping of exon 7 and a reduction in full-length transcripts in SMN23 ( Figures 1B to 1D ). The genomic region undergoes unequal crossing-over and gene conversion, resulting in variable copy numbers of SMN1 and SMN2 ( Figure 1B ). Due to the high incidence and disease severity, population-wide SMA screening is recommended, and the key to this screening is to determine the copy number of SMN1 for SMA diagnosis and carrier testing. Additionally, the copy number of SMN2 defines the severity of SMA and is important for clinical classification and prognosis.

[0089] Conventional SMA carrier testing uses PCR-based methods such as multiplex ligation-dependent probe amplification (MLPA), quantitative PCR (qPCR), and digital PCR. These methods mainly target the c.840C>T locus. Incorporating SMA screening into high-throughput NGS-based testing that allows profiling of large numbers of genes or even the entire genome may be advantageous. The almost complete sequence identity between SMN1 and SMN2 makes variant calling challenging for standard GSS-based methods.

[0090] Disclosed herein is an SMN copy number caller based on a bioinformatics method that utilizes whole-genome sequencing (WGS) data to determine the copy numbers of SMN1 and SMN2 ( Figure 1E)。The method may include calling the copy number of SMN1+SMN2 in two regions (exons 1 to 6 and exons 7 to 8) by adding the reads in SMN1 and SMN2. The method may include using the read counts at fixed base differences to distinguish SMN1 from SMN2. In some embodiments, the method does not include realigning the aligned sequences to a modified reference. The method is the first SMN copy number calling tool capable of identifying both patients and carriers with SMA from WGS data. Some embodiments of the method are not limited to exons 7 and 8 and do not primarily focus on c.840C>T. The method employs a whole-genome approach and provides the most comprehensive set of calls, including the copy numbers of full-length SMN1 and SMN2, as well as truncated forms of SMN lacking exons 7 and 8. The method can be easily applied to any WGS data and will be a valuable tool for SMA diagnosis and carrier screening to incorporate into high-throughput population-wide WGS screening.

[0091] Figures 1A to 1E Shows the reasons for SMA SMN copy number calling according to an embodiment of the bioinformatics method disclosed herein. Table 1 shows the discrimination of SMN1 from SMN2 based on fixed single nucleotide polymorphisms (SNPs) according to an embodiment of the method. SMN1 copy number calling is performed at 16 sites near c.840C>T. Nine sites with a high percentage of identity to c.840C>T are selected for a combined call of the SMN1 copy number. Figures 2A to 2C And Table 2 shows the overall distribution of the determined SMN1 / 2 copy numbers. When there are fewer copies of SMN2, more copies of SMN1 are observed, indicating that gene conversion is the mechanism of CN variability between SMN1 and SMN2. Table 3 shows the validation of the copy number calls determined using the bioinformatics method relative to those determined using digital PCR. The validation for digital PCR shows 100% concordance in SMN1 CN and 98% concordance in SMN2 CN. Figure 3 Shows SMA identified in two trios of the Next Generation Children's Project and verified using MLPA. Figure 4 And Table 4 shows that the population frequencies determined using the bioinformatics method are consistent with previous studies.

[0092] Table 1. Distinguishing SMN1 from SMN2 based on fixed single nucleotide polymorphisms (SNPs) 。

[0093]

[0094] Table 2. Population distribution of SMN1 / 2 copy numbers

[0095]

[0096] Table 3. Validation of copy number calls determined using bioinformatics methods relative to those determined using digital PCRValidation .

[0097]

[0098] Table 4. Population frequencies determined using bioinformatics methods are consistent with previous studies .

[0099]

[0100] a Hendrickson et al. “Differences in SMN1 allele frequencies among ethnic groups within North America”, J Med Genet., Vol. 46, No. 9, 2009: pp. 641-644. doi: 10.1136 / jmg.2009.066969.

[0101] b Sugarman et al. “Pan-ethnic carrier screening and prenatal diagnosis for spinal muscular atrophy: clinical laboratory analysis of >72 400 specimens”, Eur J Hum Genet., Vol. 20, No. 1, 2012: pp. 27-32. doi: 10.1038 / ejhg.2011.134.

[0102] * African Americans

[0103] Characterizing medically actionable variants in 2,500 publicly available high-depth genomes of diverse ancestries

[0104] Whole-genome sequencing (WGS) data for population-scale is increasingly available. For example, public sequence data of >2,500 samples from the 1000 Genomes Project (1kGP), such as high-depth (>30x) WGS data, are available. This has greatly improved the clinical interpretation of simple single nucleotide variants (SNVs) and insertions / deletions (indels). However, many medically important regions and variants, such as triplet repeats and paralogs, are not included in WGS-based databases because annotating these regions and variants requires specialized bioinformatics methods. For this reason, population-level characterization of known clinical variants is needed to maximize the impact of population sequencing experiments. In some embodiments, the methods disclosed herein address three drawbacks of standard secondary analysis pipelines: 1) detection and carrier screening for spinal muscular atrophy (SMA), 2) CYP2D6 genotyping for pharmacogenetic applications, and 3) detection of triplet repeat expansions. The method can be targeted to call the copy number of SMN1 / 2, CYP2D6 star alleles, and repeat expansions in the 1kGP population and quantify the differences between subsets. The frequency distributions of subsets and vertical validation of these methods using validation data generated from high-quality long reads are described herein.

[0105] CYP2D6

[0106] CYP2D6 is an important drug-metabolizing enzyme that is highly polymorphic ( Figure 5 ). CYP2D6 has high sequence similarity with its pseudogene paralog (CYP2D7). Genotyping CYP2D6 using WGS is challenging due to the common gene conversion between CYP2D6 and CYP2D7 (hereinafter referred to as CYP2D6 / 7), common SVs (gene deletions, duplications, and CYP2D6 / 7 fusion genes; see Figure 6 for illustration), and sequence similarity between CYP2D / 7, which results in ambiguous read alignments to either gene ( Figure 5 ). Disclosed herein is a CYP2D6 caller based on bioinformatics methods that can call (e.g., unambiguously call) haplotypes targeted to star alleles with known functions (e.g., all star alleles). In some embodiments, the method includes the following actions

[0107] 1. Call the total copy number of CYP2D6 + CYP2D7.

[0108] 2. Call CNV / hybrids based on copy number calls across CYP2D6 / CYP2D7 discrimination sites.

[0109] 3. Call 56 SNPs / indels from BAM (or another file containing sequence reads).

[0110] - Use copy number information.

[0111] - Count reads at the CYP2D6 and CYP2D7 positions in the homologous region.

[0112] 4. Call star alleles and haplotypes based on all called variants.

[0113] Table 5 shows the validation results of CYP2D6 star allele calls performed by this method. The CYP2D6 star allele calls performed on 92 out of 96 samples by this method were consistent with the GeT-RM consensus sequence calls from multiple platforms. This method outperformed callers such as Aldy (which called 89 out of 96 samples for CYP2D6 star alleles, consistent with the GeT-RM consensus sequence) and Stargazer (which called 83 out of 96 samples for CYP2D6 star alleles, consistent with the GeT-RM consensus sequence).

[0114] Table 5. CYP2D6 caller validation .

[0115] Sample CYP2D6 call GeT-RM consensus sequence Aldy Stargazer NA24008 *1 / *4+*68 *1 / *4 *1 / *4+*68 *1 / *4+*68 NA21781 *2×2 / *4+*68 *2×2 / *68+*2 *2×2 / *4+*68 *2×2 / *4+*68 NA23874 *4 / *4+*68 *4 / *4 No call *4 / *4+*68 NA18565 *10 / *10+*36 *10 / *36x2 *10 / *10+*36 *10 / *10+*36

[0116] Figure 7 Shows that the allele frequencies determined by this method are consistent with the PharmVar database from the Pharmacogene Variation (PharmVar) Consortium.

[0117] Determining copy number of the survival motor neuron 1 gene using sequencing data

[0118] Figure 8 Is a flowchart of an exemplary method 800 for determining the copy number of the survival motor neuron 1 gene using sequencing data such as whole genome sequencing data. Method 800 may be included 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, Figure 11 The computing system 1100 shown and described in more detail below can execute a set of executable program instructions to implement method 800. When method 800 is initiated, the executable program instructions can be loaded into a memory such as RAM and executed by one or more processors of the computing system 1100. Although method 800 is described with respect to Figure 11 The computing system 1100 shown, this description is merely exemplary and not intended to be limiting. In some embodiments, method 800 or portions thereof can be executed serially or in parallel by multiple computing systems.

[0119] After method 800 starts at block 804, method 800 proceeds to block 808, where the computing system (such as the reference Figure 11The described computing system 1100 determines: (i) a first quantity of sequence reads of the plurality of sequence reads aligned with at least one of exons 1 to 6 of a survival motor neuron 1 (SMN1) gene or a survival motor neuron 2 (SMN2) gene, respectively, and (ii) a second quantity of sequence reads of the plurality of sequence reads aligned with at least one of exons 7 and 8 of the SMN1 gene or the SMN2 gene, respectively. The first quantity of sequence reads aligned with the first SMN1 or SMN2 region (or the second quantity of sequence reads aligned with the second SMN1 or SMN2 region) can be or be about, for example, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 or higher.

[0120] At least one of exons 1 to 6 of the SMN1 gene may include exon 1, exon 2, exon 3, exon 4, exon 5, and / or exon 6 of the SMN1 gene. At least one of exons 1 to 6 of the SMN2 gene may include exon 1, exon 2, exon 3, exon 4, exon 5, and / or exon 6 of the SMN2 gene. The first SMN1 or SMN2 region may respectively contain exons 1 to 6 of the SMN1 gene or the SMN2 gene, and the length may be about 22.2 kb. The second SMN1 or SMN2 region may respectively contain exons 7 and 8 of the SMN1 gene or the SMN2 gene, and the length may be about 6 kb.

[0121] In some embodiments, the computing system receives sequence data that includes a plurality of sequence reads obtained from a sample of a subject and aligned with the SMN1 gene or the SMN2 gene. The sequencing data may include whole genome sequencing (WGS) data or short read WGS data. In some embodiments, the subject is a neonatal subject, a pediatric subject, an adolescent subject, or an adult subject. The sample may contain cellular or cell-free DNA.

[0122] In some embodiments, the sequence reads of the plurality of sequence reads are aligned with the first SMN1 or SMN2 region or the second SMN1 or SMN2 region, wherein the alignment quality score is about zero. The alignment quality can be or be about, for example, 0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10 or higher (on a scale of 0 to 1 of the alignment score).

[0123] Method 800 proceeds from block 808 to block 812, where the computing system determines (i) a first normalized quantity of sequence reads aligned to a first SMN1 or SMN2 region and (ii) a second normalized quantity of sequence reads aligned to a second SMN1 or SMN2 region using (i) the length of the first SMN1 or SMN2 region and (ii) the length of the second SMN1 or SMN2 region, respectively. The first normalized quantity of sequence reads aligned to the first SMN1 or SMN2 region (or the second normalized quantity of sequence reads aligned to the second SMN1 or SMN2 region) can be or be approximately, for example, 1, 2, 3, 4, 5, 6, 7, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or higher. The length of the first SMN1 or SMN2 region can be or be approximately, for example, 3 kb, 6 kb, 9 kb, 12 kb, 15 kb, 18 kb, 21 kb, 22.2 kb, 24 kb or a longer length. The length of the second SMN1 or SMN2 region can be or be approximately, for example, 3 kb, 6 kb or a longer length.

[0124] In some embodiments, to determine (i) a first normalized quantity of sequence reads aligned to a first SMN1 or SMN2 region and (ii) a second normalized quantity of sequence reads aligned to a second region, the computing system can use (i) the length of the first SMN1 or SMN2 region and (ii) the length of the second SMN1 or SMN2 region to determine (i) a first normalized quantity of sequence reads aligned to the first SMN1 or SMN2 region and (ii) a second normalized quantity of sequence reads aligned to the second SMN1 or SMN2 region, and to determine (iii) the depth of sequence reads of a region of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene in the sequence data. The depth of sequence reads of a region of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene in the sequence data can be or be approximately, for example, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or higher.

[0125] To determine (i) a first normalized quantity of sequence reads aligned to a first SMN1 or SMN2 region and (ii) a second normalized quantity of sequence reads aligned to a second SMN1 or SMN2 region, a computing system uses (i) the length of the first SMN1 or SMN2 region and (ii) the length of the second SMN1 or SMN2 region to determine (i) a first SMN1 or SMN2 region length-normalized quantity of sequence reads aligned to the first SMN1 or SMN2 region and (ii) a second SMN1 or SMN2 region length-normalized quantity of sequence reads aligned to the second SMN1 or SMN2 region. The computing system can use the depth of sequence reads of a region of the genome of a subject other than the locus containing the SMN1 gene and the SMN2 gene to determine (i) a first normalized depth of sequence reads aligned to the first SMN1 or SMN2 region and (ii) a second normalized depth of sequence reads aligned to the second SMN1 or SMN2 region, based on (i) the first SMN1 or SMN2 region length-normalized quantity and (ii) the second SMN1 or SMN2 region length-normalized quantity, respectively. The first normalized quantity of sequence reads aligned to the first SMN1 or SMN2 region and the second normalized quantity of sequence reads aligned to the second SMN1 or SMN2 region can be the first normalized depth and the second normalized depth, respectively.

[0126] In some embodiments, to determine (i) a first normalized number of sequence reads aligned to a first SMN1 or SMN2 region and (ii) a second normalized number of sequence reads aligned to a second region, a computing system may use (i) the GC content of the first SMN1 or SMN2 region and (ii) the GC content of the second SMN1 or SMN2 region, respectively, to determine (i) the first normalized number of sequence reads aligned to the first SMN1 or SMN2 region and (ii) the second normalized number of sequence reads aligned to the second SMN1 or SMN2 region, and to determine (iii) the depth of sequence reads of a region of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene in the sequence data, and to determine (iv) the GC content of the region of the genome. The GC content of the first SMN1 or SMN2 region (or the GC content of the second SMN1 or SMN2 region) may be or be about, for example, 40%, 41%, 42%, 43%, 44%, 45%, 46%, 47%, 48%, 49%, 50%, 51%, 52%, 53%, 54%, 55%, 56%, 57%, 58%, 59% or 60%. The depth of sequence reads of a region of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene in the sequence data may be or be about, for example, 3, 4, 5, 10, 20, 30, 40, 50, 100 or higher. The GC content of a region of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene in the sequence data may be or be about, for example, 40%, 41%, 42%, 43%, 44%, 45%, 46%, 47%, 48%, 49%, 50%, 51%, 52%, 53%, 54%, 55%, 56%, 57%, 58%, 59% or 60%.

[0127] In some embodiments, the depth of the region includes the average depth of sequence reads of a region of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene in the sequencing data. The depth of the region may include the median depth of sequence reads of a region of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene in the sequencing data. The depth of the region may be or be about, for example, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or higher. The region may contain about 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000 or more preselected regions each having a length of about 0.5 kb, 1 kb, 1.5 kb, 2 kb, 2.5 kb or 3 kb across the subject's genome. For example, the region may contain about 3000 preselected regions each having a length of about 2 kb and each spanning the subject's genome.

[0128] In some embodiments, the first normalized number of sequence reads aligned to the first SMN1 or SMN2 region (or the second normalized number of sequence reads aligned to the second SMN1 or SMN2 region) is or is about 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, or higher. For example, (i) the first normalized number of sequence reads aligned to the first SMN1 or SMN2 region and / or (ii) the second normalized number of sequence reads aligned to the second SMN1 or SMN2 region is about 30 to about 40.

[0129] Method 800 proceeds from block 812 to block 816, where, taking into account (i) the first normalized number of sequence reads aligned to the first SMN1 or SMN2 region and (ii) the second normalized number of sequence reads aligned to the second SMN1 or SMN2 region, respectively, the computational system uses a Gaussian mixture model that includes a plurality of Gaussian functions each representing a different integer copy number to determine (i) the copy number of the total survival motor neuron (SMN) gene and (ii) the copy number of any full-length SMN gene. The total survival motor neuron gene can include a full-length SMN1 gene, a full-length SMN2 gene, a truncated SMN1 gene, and / or a truncated SMN2 gene. Any full-length SMN gene can include a full-length SMN1 gene and / or a full-length SMN2 gene. The copy number of the total SMN gene (or any gene of the present disclosure) can be or about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or higher. The copy number of any full-length SMN gene (or any gene of the present disclosure) can be or about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or higher.

[0130] In some embodiments, the Gaussian mixture model includes a one-dimensional Gaussian mixture model. The plurality of Gaussian functions of the Gaussian mixture model can represent integer copy numbers, such as 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 plurality of Gaussian functions of the Gaussian mixture model can represent integer copy numbers from 0 to 10. The mean of each of the plurality of Gaussian functions (e.g., 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or greater) can be the integer copy number represented by the Gaussian function (e.g., a copy number of 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or higher). The standard deviation of the Gaussian function can be or about, for example, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, or higher.

[0131] In some embodiments, to determine (i) the copy number of the total SMN gene and (ii) the copy number of any full-length SMN gene, the computational system may consider, respectively, (i) the first normalized number of sequence reads aligned to the first SMN1 or SMN2 region and (ii) the second normalized number of sequence reads aligned to the second SMN1 or SMN2 region, and use a Gaussian mixture model and a first predetermined posterior probability threshold to determine (i) the copy number of the total SMN gene and (ii) the copy number of any full-length SMN gene. The first predetermined posterior probability threshold (or any predetermined posterior probability threshold of the present disclosure) may be or about, for example, 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 higher. For example, the first predetermined posterior probability threshold may be 0.95.

[0132] Method 800 proceeds from block 816 to block 820, where, for one base among a plurality of SMN1 gene-specific bases (also referred to herein as SMN discriminative bases) associated with a full-length SMN1 gene, the computational system determines the most likely combination among a plurality of possible combinations that each include a possible copy number of SMN1 genes and a possible copy number of SMN2 genes that together total the determined copy number of any full-length SMN genes, taking into account (a) the number of sequence reads of the plurality of sequence reads that have a base supporting the SMN1 gene-specific base (e.g., the unnormalized or normalized number of sequence reads) and (b) the number of sequence reads of the plurality of sequence reads that have a base supporting the SMN2 gene-specific base corresponding to the SMN1 gene-specific base of the SMN2 gene. The possible copy number of SMN1 genes may be or about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or higher. The possible copy number of SMN2 genes may be or about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or higher.

[0133] In some embodiments, considering (a) the number of sequence reads of the plurality of sequence reads having bases supporting SMN1 gene-specific bases and (b) the number of sequence reads of the plurality of sequence reads having bases supporting corresponding SMN2 gene-specific bases, the most likely combination of the possible copy number of the SMN1 gene and the possible copy number of the SMN2 gene is associated with the highest posterior probability relative to other combinations in the plurality of combinations. The highest posterior probability (or any probability of the present disclosure) can be or be about, for example, 60%, 61%, 62%, 63%, 64%, 65%, 66%, 67%, 68%, 69%, 70%, 71%, 72%, 73%, 74%, 75%, 76%, 77%, 78%, 79%, 80%, 81%, 82%, 83%, 84%, 85%, 86%, 87%, 88%, 89%, 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, 99% or higher. The difference in the posterior probability (or any probability of the present disclosure) can be or be about, for example, 1%, 2%, 3%, 4%, 5%, 6%, 7%, 8%, 9%, 10%, 11%, 12%, 13%, 14%, 15%, 16%, 17%, 18%, 19%, 20%, 21%, 22%, 23%, 24%, 25%, 26%, 27%, 28%, 29%, 30% or higher.

[0134] In some embodiments, to determine the most likely combination of the possible copy numbers of the SMN1 gene and the possible combinations of the SMN2 gene, the computing system can determine the most likely combination among the multiple possible combinations of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene, each including a total copy number of any complete SMN gene determined, taking into account the ratio of (a) the number of sequence reads of the multiple sequence reads having bases supporting the SMN1 gene-specific bases to (b) the number of sequence reads of the multiple sequence reads having bases supporting the SMN2 gene-specific bases corresponding to the SMN1 gene-specific bases. To determine the most likely combination of the possible copy numbers of the SMN1 gene and the possible combinations of the SMN2 gene, the computing system can determine (a) the number of sequence reads of the multiple sequence reads having bases supporting the SMN1 gene-specific bases and (b) the number of sequence reads of the multiple sequence reads having bases supporting the SMN2 gene-specific bases corresponding to the SMN1 gene-specific bases. The computing system can determine the ratio of (a) the number of sequence reads of the multiple sequence reads having bases supporting the SMN1 gene-specific bases to (b) the number of sequence reads of the multiple sequence reads having bases supporting the SMN2 gene-specific bases corresponding to the SMN1 gene-specific bases. Based on the ratio of (a) the number of sequence reads of the multiple sequence reads having bases supporting the SMN1 gene-specific bases to (b) the number of sequence reads of the multiple sequence reads having bases supporting the SMN2 gene-specific bases corresponding to the SMN1 gene-specific bases, the computing system can determine the most likely combination among the multiple possible combinations of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene, each including a total copy number of any complete SMN gene determined.

[0135] In some embodiments, to determine the most likely combination of the possible copy numbers of the SMN1 gene and the possible combinations of the SMN2 gene, for each of the plurality of SMN1 gene-specific bases, the computing system determines the most likely combination associated with the highest posterior probability among the plurality of possible combinations of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene, each including a total number of copy numbers of the SMN1 gene and the SMN2 gene that is the copy number of any complete SMN gene determined, taking into account (a) the number of sequence reads of the plurality of sequence reads having bases that support the SMN1 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads having bases that support the SMN2 gene-specific base corresponding to the SMN1 gene-specific base of the SMN2 gene. The number of sequence reads aligned to the SMN1 gene-specific base (or SMN2 gene-specific base) can be or be about, for example, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or higher. To determine the copy number of the SMN1 gene, the computing system can determine the copy number of the SMN1 gene based on the possible copy number of the SMN1 gene in the most likely combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene determined for each of the plurality of SMN1 gene-specific bases.

[0136] In some embodiments, the SMN1 gene-specific base is a splicing enhancer. The SMN1 gene-specific base can be the base at c.840 of the SMN1 gene. In some embodiments, the SMN1 gene-specific base has identity with each of the plurality of SMN1 gene-specific bases other than the SMN1 gene-specific base that exceeds a predetermined identity threshold. The predetermined identity threshold (or any threshold of the present disclosure) can be or be about, for example, 80%, 81%, 82%, 83%, 84%, 85%, 86%, 87%, 88%, 89%, 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, 99% or higher. For example, the identity threshold can be 97%. The plurality of SMN1 gene-specific bases can comprise or comprise about 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21 or more SMN1 gene-specific bases. For example, the plurality of SMN1 gene-specific bases can include 8 SMN1 gene-specific bases. Each of the plurality of SMN1 gene-specific bases can be located on intron 6, exon 7, intron 7 or exon 8 of the SMN1 gene.

[0137] If the subject is of a first race (or ethnic group), the plurality of SMN1 gene-specific bases may be different, if the subject is of a second race (or ethnic group), the plurality of SMN1 gene-specific bases may be different, and if the subject is of an unknown race, the plurality of SMN1 gene-specific bases may be different. The race can be, for example, Caucasian, African, African American, American Indian, Alaska Native, Asian, South Asian, East Asian, Native Hawaiian, Pacific Islander, or a combination thereof. The race (or ethnic group) of the subject may be unknown, and the plurality of SMN1 gene-specific bases may not be race-specific (or ethnic group-nonspecific). The race (or ethnic group) of the subject may be known, and the plurality of SMN1 gene-specific bases may be specific to the subject's race (or ethnic group). In some embodiments, the computing system may receive race (or ethnic group) information of the subject. The computing system may select the plurality of SMN1 gene-specific bases from the plurality of SMN1 gene-specific bases based on the received race (or ethnic group) information.

[0138] Method 800 proceeds from block 820 to block 824, where the computing system determines the copy number of the SMN1 gene using the most likely combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene determined for the SMN1 gene-specific bases. Alternatively or in addition, the computing system determines the copy number of the SMN2 gene using the most likely combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene determined for the SMN1 gene-specific bases.

[0139] In some embodiments, to determine the copy number of the SMN1 gene, the computing system may use the most likely combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene determined for each base among the plurality of SMN1 gene-specific bases to determine the copy number of the SMN1 gene and the copy number of the SMN2 gene. To determine the copy number, the computing system may use the most likely combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene determined for the SMN1 gene-specific bases and a second predetermined posterior probability threshold of the combinations of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene to determine the copy number of the SMN1 gene. The second predetermined posterior probability threshold (or any predetermined posterior probability threshold of the present disclosure) may be or be about, for example, 0.50, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.60, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 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 higher. For example, the second predetermined posterior probability threshold may be 0.6 or 0.8.

[0140] In some embodiments, most of the possible copy numbers of the determined SMN1 gene are consistent. The copy number of the determined SMN1 gene may be the consistent possible copy number of the SMN1 gene. The computing system may determine the possible combinations of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene that include a total of the copy numbers of any complete SMN genes determined, considering (a) the number of sequence reads of the plurality of sequence reads having bases that support any base among the SMN1 gene-specific bases and (b) the number of sequence reads of the plurality of sequence reads having bases that support any one of the plurality of corresponding SMN2 gene-specific bases. The computing system may determine the possible copy number of the possible combination to be the consistent possible copy number of the SMN1 gene.

[0141] In some embodiments, to determine the copy number of the SMN1 gene, a computing system may determine that the copy number of the SMN1 gene is zero, one, or more than one. In some embodiments, the computing system may determine the spinal muscular atrophy (SMA) status of a subject based on the copy number of the SMN1 gene. The SMA status of the subject may include SMA, SMA carrier but not SMA, and not an SMA carrier. In some embodiments, the computing system may use the number of sequence reads of the plurality of sequence reads aligned to g.27134 of the SMN1 gene and the bases of the sequence reads aligned to g.27134 of the SMN1 gene to determine that the subject is a silent SMA carrier.

[0142] For example, the computing system at block 820 may identify the ratio of read overlap positions of SMN1 and SMN2, where the genes have different base sequences. For positions where SMN1 differs from SMN2, the computing system may extract the overlapping reads based on SMN1 or SMN2. From these reads, the computing system may count the number of SMN1 - specific bases and the number of SMN2 - specific bases. The computing system may determine the fraction of SMN1 or SMN2 reads. The computing system may calculate the CN of SMN1 and SMN2 at positions where SMN1 differs from SMN2. The computing system may combine the full - length CN with the ratio of SMN1 to SMN2 to call the CN of SMN1 and the CN of SMN2. The computing system at block 824 may combine the CNs from multiple fixed differences between SMN1 and SMN2 to obtain the accurate CN of SMN1 and the accurate CN of SMN2.

[0143] Calling SMA / non-SMA or carrier / non-carrier In some embodiments, the computing system may use (i) the determined copy number of the total SMN gene and (ii) the determined copy number of the full - length SMN gene to determine the copy number of the truncated SMN gene. The copy number of the truncated SMN gene may be the difference between (i) the determined copy number of the total SMN gene and (ii) the determined copy number of the full - length SMN gene.

[0144] Treatment In some embodiments, the computing system may determine treatment recommendations for a subject based on the determined copy number of the SMN1 gene. The treatment recommendations may include administering Nusinersen and / or Zolgensma to the subject.

[0145] Method 800 ends at block 828.

[0146] Genotyping the cytochrome P450 family 2 subfamily D member 6 gene using sequencing data

[0147] Figure 9A flowchart showing an exemplary method 900 for genotyping the cytochrome P450 family 2 subfamily D member 6 gene using sequencing data such as whole genome sequencing data. Method 900 may be included in a set of executable program instructions stored on a computer-readable medium of a computing system such as one or more disk drives. For example, Figure 11 The computing system 1100 shown and described in more detail below may execute a set of executable program instructions to implement method 900. When method 900 is initiated, the executable program instructions may be loaded into a memory such as RAM and executed by one or more processors of the computing system 1100. Although method 900 is described with respect to Figure 11 the computing system 1100 shown, this description is merely exemplary and not intended to be limiting. In some embodiments, method 900 or portions thereof may be executed serially or in parallel by multiple computing systems.

[0148] The number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene (e.g., the unnormalized or normalized number of sequence reads) may be used to determine the total copy number (CN) of the CYP2D6 gene and the CYP2D7 gene using a Gaussian mixture model. The total CN of the CYP2D6 gene and the CYP2D7 gene may be used to determine the CN of CYP2D6 at various CYP2D6 / CYP2D7 discriminatory bases (also referred to herein as CYP2D6 gene-specific bases) by iterating over all possible combinations of CYP2D6 CN and CYP2D7 CN at the CYP2D6 / CYP2D7 discriminatory bases. The CYP2D6 CN at various CYP2D6 / CYP2D7 discriminatory bases may be used to call structural variants. For example, at each CYP2D6 / CYP2D7 discriminatory base (also referred to herein as a CYP2D6 gene-specific base), the number of chromosomes carrying the CYP2D6 gene and the number of chromosomes carrying the CYP2D7 gene may be called by combining the total CN of the CYP2D6 gene and the CYP2D7 gene with the read counts supporting each of the gene-specific bases. Based on the total CN of the schedule, all possible combinations of CYP2D6 CN and CYP2D7 CN may be iterated to derive the combination that produces the highest posterior probability for the observed number of reads supporting CYP2D6 and CYP2D7. Structural variants may be called by identifying bases where the CN of the CYP2D6 gene has changed.

[0149] One or more minor variants can be determined. For each minor variant position of the minor variants, the minor variants can be determined by iterating over all possible combinations of the variant allele CN and the reference (non-variant) allele CN to determine the most likely variant allele CN using sequence reads that overlap the minor variant positions in the CYP2D6 or CYP2D7 gene. For example, if there are a total of three CYP2D6 gene copies and there are 10 reads supporting the variant allele and 20 reads supporting the reference allele, the variant allele CN can be determined to be one, i.e., there is one CYP2D6 gene copy carrying the minor variant. For example, minor variants that define a star allele can be searched for in the sequence data (e.g., in a BAM file). The minor variants of interest can be divided into those that belong to the CYP2D6 / CYP2D7 homologous regions and those that do not. For the former, variant reads that align with the CYP2D6 gene or the CYP2D7 gene and overlap each minor variant position of the CYP2D6 gene of interest or the corresponding position in the CYP2D7 gene can be searched for. For the latter, reads that align with the CYP2D6 gene and overlap the minor variant positions of the CYP2D6 gene of interest can be searched for. The CN called in this region can also be considered during the minor variant call. The called structural variants and minor variants can be matched to the definition of the star allele to name the star alleles, and these star alleles can be further grouped into haplotypes.

[0150] After method 900 starts at block 904, method 900 proceeds to block 908, where the computing system (e.g., the computing system 1100 described in the reference Figure 11 determines (i) a first quantity of sequence reads of a plurality of sequence reads that align with the cytochrome P450 family 2 subfamily D member 6 (CYP2D6) gene or the cytochrome P450 family 2 subfamily D member 7 (CYP2D7) gene. The first quantity of sequence reads that align with the first CYP2D6 gene or CYP2D7 gene (or any gene of the present disclosure) can be or can be about, for example, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 or higher.

[0151] A computing system can receive sequence data that includes multiple sequence reads obtained from a subject and aligned to the CYP2D6 gene or the CYP2D7 gene. In some embodiments, the sequencing data includes whole genome sequencing (WGS) data or short-read WGS data. The subject can be a neonatal subject, a pediatric subject, an adolescent subject, or an adult subject. The sample can include cellular or cell-free DNA. The sample can include cellular or cell-free DNA.

[0152] In some embodiments, the sequence reads of the multiple sequence reads are aligned to the CYP2D6 gene or the CYP2D7 gene, where the alignment quality score is about zero. The alignment quality can be or be about, for example, 0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, or higher (on a scale of 0 to 1 of the alignment score).

[0153] In some embodiments, to determine (i) a first quantity of the sequence reads of the multiple sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, the computing system can determine (i) a first quantity of the sequence reads of the multiple sequence reads aligned to at least one exon or intron of the CYP2D6 gene (e.g., one of exons 1 to 9 or one of introns 1 to 8 of the CYP2D6 gene) and / or at least one exon or intron of the CYP2D7 gene (e.g., one of exons 1 to 9 or one of introns 1 to 8 of the CYP2D7 gene).

[0154] Method 900 proceeds from block 908 to block 912, where the computing system determines (i) a first normalized quantity of the sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene using (i) the length of the CYP2D6 gene or the CYP2D7 gene, respectively. The first normalized quantity of the sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene (or any gene of the present disclosure) can be or be about, for example, 1, 2, 3, 4, 5, 6, 7, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, or higher. The length of the CYP2D6 gene can be or be about, for example, 4.4 kb. The length of the CYP2D7 gene can be or be about, for example, 4.9 kb.

[0155] In some embodiments, to determine (i) a first normalized quantity of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, a computing system may use (i) the length of the CYP2D6 gene or the CYP2D7 gene, respectively, to determine (i) the first normalized quantity of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, and to determine (iii) the depth of sequence reads in a region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene. The depth of sequence reads in a region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene (or any gene of the present invention) may be or may be about, for example, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or higher.

[0156] To determine (i) a first normalized quantity of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene and (ii) a second normalized quantity of sequence reads aligned to a second region, a computing system may use (i) the length of the CYP2D6 gene or the CYP2D7 gene, respectively, to determine (i) a first CYP2D6 gene or CYP2D7 gene length-normalized quantity of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene. The computing system may use the depth of sequence reads in a region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene to determine a first normalized depth of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene based on (i) the CYP2D6 gene or CYP2D7 gene length-normalized quantity. The first normalized depth of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene may be the first normalized quantity of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, respectively.

[0157] In some embodiments, to determine (i) a first normalized quantity of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, a computing system can use (i) the GC content of the CYP2D6 gene or the CYP2D7 gene to determine (i) the first normalized quantity of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, and to determine (iii) the depth of sequence reads in a region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene in the sequence data, and (iv) to determine the GC content of the region of the genome. The GC content of the CYP2D6 gene or the CYP2D7 gene (or any gene of the present disclosure) can be or be about, for example, 40%, 41%, 42%, 43%, 44%, 45%, 46%, 47%, 48%, 49%, 50%, 51%, 52%, 53%, 54%, 55%, 56%, 57%, 58%, 59%, or 60%. The depth of sequence reads in a region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene in the sequence data can be or be about, for example, 3, 4, 5, 10, 20, 30, 40, 50, 100, or higher. The GC content of a region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene (or any gene of the present invention) in the sequence data can be or be about, for example, 40%, 41%, 42%, 43%, 44%, 45%, 46%, 47%, 48%, 49%, 50%, 51%, 52%, 53%, 54%, 55%, 56%, 57%, 58%, 59%, or 60%.

[0158] The depth of the region can include the average depth of sequence reads in a region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene in the sequencing data. The depth of the region can include the median depth of sequence reads in a region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene in the sequencing data. The depth of the region can be or be about 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, or higher. The region can contain about 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, or more preselected regions each having a length spanning about 0.5 kb, 1 kb, 1.5 kb, 2 kb, 2.5 kb, or 3 kb of the subject's genome. For example, the region can contain about 3000 preselected regions each having a length of about 2 kb and each spanning the subject's genome.

[0159] In some embodiments, (i) the first normalized quantity of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene and / or (ii) the second normalized quantity of sequence reads aligned to the second region is or is about 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or higher. For example, (i) the first normalized quantity of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene and / or (ii) the second normalized quantity of sequence reads aligned to the second region is about 30 to about 40.

[0160] Method 900 proceeds from block 912 to block 916, where, taking into account (i) the first normalized quantity of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, the computational system uses a Gaussian mixture model that includes a plurality of Gaussian functions each representing a different integer copy number to determine (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene. The total copy number of the CYP2D6 gene and the CYP2D7 gene (or any gene of the present invention) can be or can be about 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or higher.

[0161] In some embodiments, the Gaussian mixture model includes a one-dimensional Gaussian mixture model. The plurality of Gaussian functions of the Gaussian mixture model can represent integer copy numbers, such as 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 plurality of Gaussian functions of the Gaussian mixture model can represent integer copy numbers from 0 to 10. The mean of each of the plurality of Gaussian functions (e.g., 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or greater) can be the integer copy number represented by the Gaussian function (e.g., a copy number of 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or higher). The standard deviation of the Gaussian function can be or can be about, for example, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 or higher.

[0162] In some embodiments, to determine (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene, the computing system may use a Gaussian mixture model and a first predetermined posterior probability threshold to determine (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene, taking into account (i) the first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene. The first predetermined posterior probability threshold (or any predetermined posterior probability threshold of the present disclosure) can be or about, for example, 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 higher. For example, the first predetermined posterior probability threshold can be 0.95.

[0163] Method 900 proceeds from block 916 to block 920, where for one of a plurality of CYP2D6 gene-specific bases (also referred to herein as CYP2D6 / CYP2D7 discriminatory bases), the computing system determines the most likely combination among a plurality of possible combinations that each include a possible copy number of the CYP2D6 gene and a possible copy number of the CYP2D7 gene that together total the determined total copy number of the CYP2D6 gene and the CYP2D7 gene, taking into account (a) the number of sequence reads of the plurality of sequence reads that have bases supporting the CYP2D6 gene-specific base (e.g., the unnormalized or normalized number of sequence reads) and (b) the number of sequence reads of the plurality of sequence reads that have bases supporting the CYP2D7 gene-specific base corresponding to the CYP2D6 gene-specific base (e.g., the unnormalized or normalized number of sequence reads). The possible copy number of the CYP2D6 gene can be or about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or higher. The possible copy number of the CYP2D7 gene can be or about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or higher.

[0164] In some embodiments, considering (a) the number of sequence reads of the plurality of sequence reads having bases supporting CYP2D6 gene-specific bases and (b) the number of sequence reads of the plurality of sequence reads having bases supporting the corresponding CYP2D7 gene-specific bases, the most likely combination of the possible copy number of the CYP2D6 gene and the possible copy number of the CYP2D7 gene is associated with the highest posterior probability relative to other combinations in the plurality of combinations. The highest posterior probability (or any probability of the present disclosure) can be or be about, for example, 60%, 61%, 62%, 63%, 64%, 65%, 66%, 67%, 68%, 69%, 70%, 71%, 72%, 73%, 74%, 75%, 76%, 77%, 78%, 79%, 80%, 81%, 82%, 83%, 84%, 85%, 86%, 87%, 88%, 89%, 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, 99% or higher. The difference in the posterior probability (or any probability of the present disclosure) can be or be about, for example, 1%, 2%, 3%, 4%, 5%, 6%, 7%, 8%, 9%, 10%, 11%, 12%, 13%, 14%, 15%, 16%, 17%, 18%, 19%, 20%, 21%, 22%, 23%, 24%, 25%, 26%, 27%, 28%, 29%, 30% or higher.

[0165] In some embodiments, to determine the possible copy numbers of the CYP2D6 gene and the most likely combination of the possible copy numbers, the computing system may determine the most likely combination among the multiple possible combinations of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene, each including a total copy number equal to the determined total copy numbers of the CYP2D6 gene and the CYP2D7 gene, taking into account the ratio of (a) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D6 gene-specific bases to (b) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D7 gene-specific bases corresponding to the CYP2D6 gene-specific bases. To determine the possible copy numbers of the CYP2D6 gene and the most likely combination of the possible copy numbers, the computing system may determine (a) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D6 gene-specific bases and (b) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D7 gene-specific bases corresponding to the CYP2D6 gene-specific bases. The computing system may determine the ratio of (a) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D6 gene-specific bases to (b) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D7 gene-specific bases corresponding to the CYP2D6 gene-specific bases. Taking into account the ratio of (a) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D6 gene-specific bases to (b) the number of sequence reads of the multiple sequence reads having bases supporting the CYP2D7 gene-specific bases corresponding to the CYP2D6 gene-specific bases, the computing system may determine the most likely combination among the multiple possible combinations of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene, each including a total copy number equal to the determined total copy numbers of the CYP2D6 gene and the CYP2D7 gene.

[0166] In some embodiments, to determine the most likely combination of the possible copy numbers of the CYP2D6 gene and the possible combinations of the CYP2D7 gene, for each of the plurality of CYP2D6 gene-specific bases, the computing system determines the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene in a plurality of possible combinations, each including a total copy number of the determined CYP2D6 gene and CYP2D7 gene, taking into account (a) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D6 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D7 gene-specific base corresponding to the CYP2D6 gene-specific base of the CYP2D7 gene. The number of sequence reads aligned to the SMN1 gene-specific base (or SMN2 gene-specific base) can be or be about, for example, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or higher. To determine the alleles of the CYP2D6 gene that the subject has, the computing system can use the most likely combination of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene determined for each of the plurality of CYP2D6 gene-specific bases to determine whether the alleles of the CYP2D6 gene that the subject has are minor variants or structural variants of the CYP2D6 gene, or neither.

[0167] In some embodiments, the CYP2D6 gene-specific base is consistent with each of the plurality of CYP2D6 gene-specific bases other than the CYP2D6 gene-specific base that exceeds a pre-determined consistency threshold. The consistency threshold (or any threshold of the present disclosure) can be or be about, for example, 80%, 81%, 82%, 83%, 84%, 85%, 86%, 87%, 88%, 89%, 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, 99% or higher. For example, the pre-determined consistency threshold can be 97%. The plurality of CYP2D6 gene-specific bases can comprise or comprise about, for example, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 118, 120, 130, 140, 150, 160, 170 or more CYP2D6 gene-specific bases. For example, the plurality of CYP2D6 gene-specific bases can include 118 CYP2D6 gene-specific bases.

[0168] Method 900 proceeds from block 920 to block 924, where the computing system determines one or more resulting variants of the CYP2D6 gene that the subject has using the most likely combination of the likely copy number of the CYP2D6 gene and the likely copy number of the CYP2D7 gene determined for CYP2D6 gene-specific bases. For example, the computing system may identify the ratio of read overlap positions of CYP2D6 and CYP2D7, where the genes have different base sequences. For positions where CYP2D6 differs from CYP2D7, the computing system may extract overlapping reads based on CYP2D6 or CYP2D7. From these reads, the computing system may count the number of CYP2D6-specific bases and the number of CYP2D7-specific bases. The computing system may determine the fraction of CYP2D6 or CYP2D7 reads. The computing system may calculate the CN of CYP2D6 and CYP2D7 at positions where CYP2D6 differs from CYP2D7. The computing system may combine the total CN of CYP2D6 and CYP2D7 with the ratio of CYP2D6 to CYP2D7 to call the CN of CYP2D6 and the CN of CYP2D7. The computing system may use the CN of CYP2D6 and CYP2D7 for small variant calling under one or more fixed differences between CYP2D6 and CYP2D7. The computing system may perform structural variant calling by combining the CN of CYP2D6 and CYP2D7 with multiple fixed differences between CYP2D6 and CYP2D7 to determine the presence of a transition between the CNs at CYP2D6 and CYP2D7, which defines the type of structural variant in the sample.

[0169] Fusion genes containing REP. In some embodiments, the computing system may determine a second number of sequence reads of the plurality of sequence reads that align to a spacer region between the CYP2D7 gene and a repeat element REP7 downstream of the CYP2D7 gene. The second number of sequence reads of the plurality of sequence reads that align to the spacer region between the CYP2D7 gene and the repeat element REP7 downstream of the CYP2D7 gene may be or may be about, for example, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 or higher. The computing system may use (ii) the length of the spacer region to determine (ii) a second normalized number of sequence reads that align to the spacer region. The second normalized number of sequence reads that align to the spacer region may be or may be about, for example, 1, 2, 3, 4, 5, 6, 7, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or higher. The length of the spacer region may be or may be about, for example, 1.5 kb. Taking into account (ii) the second normalized number of sequence reads that align to the spacer region, the computing system may use a Gaussian mixture model to determine (ii) the copy number of the spacer region. The copy number of the spacer region may be or may be about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or higher. To determine the alleles of the CYP2D6 gene that a subject has, the computing system may use a combination of the possible copy number of the CYP2D6 gene, the possible copy number of the CYP2D7 gene, and the copy number of the spacer region determined for CYP2D6 gene-specific bases to determine whether the alleles of the CYP2D6 gene that the subject has are small variants or structural variants of the CYP2D6 gene, or neither. The structural variant may include a CYP2D6 / CYP2D7 fusion allele having the spacer region and a repeat element REP7 downstream of the CYP2D6 / CYP2D7 fusion allele.

[0170] Method 900 proceeds from block 924 to block 928, where for a variant position of the CYP2D6 gene that is associated with a minor variant allele of the CYP2D6 gene, taking into account (a) the number of sequence reads (e.g., unnormalized or normalized number of sequence reads) that align with the CYP2D6 gene, overlap with the variant position, and have bases that support the minor variant allele of the CYP2D6 gene at the variant position and (b) the number of sequence reads (e.g., unnormalized or normalized number of sequence reads) that align with the CYP2D6 gene, overlap with the variant position, and have bases that support the reference allele of the CYP2D6 gene at the variant position, a computing system can determine the most likely combination of the possible copy number of the minor variant allele of the CYP2D6 gene at the variant position and the possible copy number of the reference allele of the CYP2D6 gene at the variant position that together total the copy number of the CYP2D6 gene at the variant position. The possible copy number of the minor variant allele of the CYP2D6 gene in the most likely combination at the variant position can indicate the one or more minor variants of the CYP2D6 gene.

[0171] For each variant position among multiple variant positions of the CYP2D6 gene, where the variant position is associated with a minor variant allele of the CYP2D6 gene, taking into account (a) the number of sequence reads (e.g., unnormalized or normalized number of sequence reads) that align with the CYP2D6 gene, overlap with the variant position, and have bases that support the minor variant allele of the CYP2D6 gene at the variant position and (b) the number of sequence reads (e.g., unnormalized or normalized number of sequence reads) that align with the CYP2D6 gene, overlap with the variant position, and have bases that support the reference allele of the CYP2D6 gene at the variant position, a computing system can determine the most likely combination of the possible copy number of the minor variant allele of the CYP2D6 gene at the variant position and the possible copy number of the reference allele of the CYP2D6 gene at the variant position that together total the copy number of the CYP2D6 gene at the variant position. The possible copy number of the minor variant allele of the CYP2D6 gene in the most likely combination at the multiple variant positions can indicate the one or more minor variants of the CYP2D6 gene.

[0172] In some embodiments, a computing system may determine the copy number of the CYP2D6 gene at a minor variant position. The copy number of the CYP2D6 gene at a minor variant position may include the copy number of the CYP2D6 gene. The copy number of the CYP2D6 gene at a minor variant position may include the copy number of the CYP2D6 gene that is the determined most likely combination of the possible copy numbers of the CYP2D6 gene. The copy number of the CYP2D6 gene at a minor variant position may include the copy number of the CYP2D6 gene that is the determined most likely combination and closest to the minor variant position of the possible copy numbers of the CYP2D6 gene. The copy number of the CYP2D6 gene at a minor variant position may include the copy number of the CYP2D6 gene at the 5' position or the 3' position of the minor variant position.

[0173] In some embodiments, a computing system may (a) determine the number of sequence reads (e.g., the unnormalized or normalized number of sequence reads) of bases that support a minor variant allele of the CYP2D6 gene. The computing system may (b) determine the number of sequence reads (e.g., the unnormalized or normalized number of sequence reads) of bases that support a reference allele of the CYP2D6 gene.

[0174] Method 900 proceeds from block 928 to block 932, where the computing system determines one or more minor variants of the CYP2D6 gene using the determined most likely combination of the possible copy numbers of the minor variant alleles of the CYP2D6 gene. The computing system may determine one or more minor variants of the CYP2D6 gene using the possible copy numbers of the minor variant alleles of the CYP2D6 gene at the plurality of minor variant positions in the determined most likely combination.

[0175] In some embodiments, the minor variant position is in the CYP2D6 / CYP2D7 homology region. To determine the most likely combination, the computing system can determine the most likely combination including determining the most likely combination of the possible copy number of the minor variant allele at the minor variant position of the CYP2D6 gene and the possible copy number of the reference allele at the minor variant position of the CYP2D6 gene that together total the copy number of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads having bases supporting the minor variant allele at the minor variant position of the CYP2D6 gene aligned to the CYP2D6 gene or the CYP2D7 gene and / or (b) the number of sequence reads having bases supporting the reference allele at the minor variant position of the CYP2D6 gene aligned to the CYP2D6 gene or the CYP2D7 gene. In some embodiments, the minor variant position is not in the CYP2D6 / CYP2D7 homology region. To determine the most likely combination, the computing system can determine the most likely combination including determining the most likely combination of the possible copy number of the minor variant allele at the minor variant position of the CYP2D6 gene and the possible copy number of the reference allele at the minor variant position of the CYP2D6 gene that together total the copy number of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads having bases supporting the minor variant allele at the minor variant position of the CYP2D6 gene aligned to the CYP2D6 gene (and not to the CYP2D7 gene) and / or (b) the number of sequence reads having bases supporting the reference allele at the minor variant position of the CYP2D6 gene aligned to the CYP2D6 gene (and not to the CYP2D7 gene).

[0176] For example, the computing system can first determine the SV (structural variant, such as a deletion or duplication) pattern and breakpoint based on the CN of paralog-specific bases. In addition or alternatively, the computing system can then call a predefined set of minor variants (which are variants specific to the gene of interest, such as CYP2D6, and are a set of variants different from the bases distinguishing paralogs) based on the read alignment, total CN, and (sometimes) SV pattern and breakpoint determined in the first step. Since the alignment is not always accurate, the computing system can extract the bases of interest from the reads aligned to either paralog.

[0177] Method 900 proceeds from block 932 to block 936, where the computing system uses the one or more structural variants of the identified CYP2D6 gene and / or the one or more minor variants of the identified CYP2D6 gene to determine the star alleles and / or haplotypes of the CYP2D6 gene that the subject has. Star alleles can be associated with known functions. The star alleles and / or haplotypes of the CYP2D6 gene can include, for example, CYP2D6*1, *2, *3, *4, *5, *6, *7, *9, *10, *11, *13, *14, *15, *17, *21, *22, *28, *29, *31, *33, *34, *35, *36, *37, *38, *39, *40, *41, *43, *45, *46, *47, *49, *52, *54, *56, *57, *59, *64, *65, *68, *71, *72, *82, *84, *86, *94, *95, *99, *100, *101, *106, *108, *111, *112, *113, or combinations thereof.

[0178] Enzyme activity . In some embodiments, the computing system can use the identified alleles of the CYP2D6 gene to determine the level of CYP2D6 enzyme activity in the subject. The enzyme activity can be poor, moderate, normal, or ultrametabolic. The computing system can determine a treatment dose recommendation and / or a treatment recommendation for the subject based on the one or more minor variants and / or the one or more structural variants.

[0179] Method 900 ends at block 940.

[0180] Performing paralog genotyping using sequencing data

[0181] Figure 10 A flowchart showing an exemplary method 1000 for using sequencing data such as whole-genome sequencing data for paralog genotyping. Method 1000 can be included 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, Figure 11 the computing system 1100 shown and described in more detail below can execute a set of executable program instructions to implement method 1000. When method 1000 is initiated, the executable program instructions can be loaded into a memory such as RAM and executed by one or more processors of the computing system 1100. Although method 1000 is described with respect to Figure 11 the computing system 1100 shown, this description is merely exemplary and not intended to be limiting. In some embodiments, method 1000 or portions thereof can be executed serially or in parallel by multiple computing systems.

[0182] After method 1000 begins at block 1004, method 1000 proceeds to block 1008, where a computing system (e.g., computing system 1100 described with reference to Figure 11 receives sequence data that includes a sample obtained from a subject and multiple sequence reads aligned to a first paralog or a second paralog. Techniques for generating the sequence reads include sequencing-by-synthesis using, for example, MINISEQ, MISEQ, NEXTSEQ, HISEQ, and NOVASEQ sequencing instruments from Illumina, Inc., San Diego, CA.

[0183] Method 1000 proceeds from block 1008 to block 1012, where, taking into account (i) a first quantity of sequence reads aligned to a first region, the computing system uses a Gaussian mixture model that includes multiple Gaussian functions each representing a different integer copy number to determine the copy number of a first type of paralog. The copy number of a first type (or any type disclosed herein) of paralog can be or be about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or higher.

[0184] The multiple Gaussian functions of the Gaussian mixture model can represent integer copy numbers, such as 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 multiple Gaussian functions of the Gaussian mixture model can represent integer copy numbers from 0 to 10. The mean of each of the multiple Gaussian functions (e.g., 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or greater) can be the integer copy number represented by the Gaussian function (e.g., 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or higher copy number). The standard deviation of the Gaussian function can be or be about, for example, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, or higher.

[0185] In some embodiments, a computing system may determine (i) a first number of sequence reads of a plurality of sequence reads in sequence data that are from a sample obtained from a subject and aligned to a first region. The first number of sequence reads of the plurality of sequence reads in sequence data that are from a sample obtained from a subject and aligned to the first region (or any region of the present disclosure) may be or may be about, for example, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, or higher. The computing system may use (i) the length of the first region to determine (i) a first normalized number of sequence reads aligned to the first region. The first normalized number of sequence reads aligned to the first region (or any region of the present disclosure) may be or may be about, for example, 1, 2, 3, 4, 5, 6, 7, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, or higher. The length of the first region may be or may be about, for example, 1 kb, 2 kb, 3 kb, 4 kb, 5 kb, 6 kb, 7 kb, 8 kb, 9 kb, 10 kb, 11 kb, 12 kb, 13 kb, 14 kb, 15 kb, 16 kb, 17 kb, 18 kb, 19 kb, 20 kb, 21 kb, 22 kb, 23 kb, 24 kb, 25 kb, 26 kb, 27 kb, 28 kb, 29 kb, 30 kb, or higher. To determine the copy number of a first type of paralog, the computing system may use a Gaussian mixture model to determine the copy number of the first type of paralog, taking into account (i) the first normalized number of sequence reads aligned to the first region.

[0186] In some embodiments, the computing system can use Gaussian mixture to determine the copy number of one or more paralogs of a second type, given (ii) the second number of sequence reads aligned to the second region. The copy number of the one or more paralogs of the second type (or any type of the present disclosure) can be or can be about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or higher. To determine the copy number or allele of a first paralog, the computing system can use the most likely combination of the possible copy numbers of the first paralog and the possible copy numbers of the second paralog determined for the first paralog-specific bases, along with the copy number of the one or more paralogs of the second type, to determine the copy number or allele of the first paralog. The computing system can determine the copy number of a third type of paralog from the copy numbers of the first type of paralog and the second type of paralog. The copy number of the third type of paralog (or any type of the present disclosure) can be or can be about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or higher. To determine the copy number of the first paralog, the computing system can use the most likely combination of the possible copy numbers of the first paralog and the possible copy numbers of the second paralog determined for the first paralog-specific bases to determine the copy number or allele of the first paralog.

[0187] Methods for aligning sequence reads to a reference genomic sequence can utilize aligners such as the Burrows-Wheeler aligner (BWA) and iSAAC. Other alignment methods include 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&NovoalignCS, NextGENe, Omixon, PALMapper, Partek, PASS, PerM, PRIMEX, QPalma, RazerS, REAL, cREAL, RMAP, rNA, RT 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.

[0188] Method 1000 proceeds from block 1012 to block 1016, where, for one of the plurality of first paralog-specific bases, the computing system determines the most likely combination among a plurality of possible combinations that each include a possible copy number of a first type of first paralog and a possible copy number of a first type of second paralog that together total the determined copy number of the first type of paralog, taking into account (a) the number of sequence reads of the plurality of sequence reads having bases that support the first paralog-specific base and (b) the number of sequence reads of the plurality of sequence reads having bases that support a second paralog-specific base corresponding to the first paralog-specific base for the second paralog.

[0189] Method 1000 proceeds from block 1016 to block 1020, where the computing system uses the most likely combination of the possible copy numbers of the first paralog and the second paralog determined for the first paralog-specific base to determine the copy number or allele of the first paralog.

[0190] In some embodiments, the first paralog is the survival motor neuron 1 (SMN1) gene. The second paralog can be the survival motor neuron 2 (SMN2) gene. The first region can include at least exons 1 to 6 of the SMN1 gene and at least exons 1 to 6 of the SMN2 gene. The second region can include at least one of exons 7 and 8 of the SMN1 gene and at least one of exons 7 and 8 of the SMN2 gene. The first type of paralog can include the complete SMN1 gene and the complete SMN2 gene. The second type of the one or more paralogs can include the complete SMN1 gene, the complete SMN2 gene, a truncated SMN1 gene, or a truncated SMN2 gene. The copy number of the first paralog can include the copy number of the SMN1 gene. The computing system can determine to achieve the reference Figure 8 the copy number of the SMN1 gene for the method 800 (or a portion thereof) described.

[0191] In some embodiments, the first paralog is the cytochrome P450 family 2 subfamily D member 6 (CYP2D6) gene. The second paralog can be the cytochrome P450 family 2 subfamily D member 7 (CYP2D7) gene. The first region can include the CYP2D6 gene and the CYP2D7 gene. The second region can include the CYP2D7 gene and the spacer region between the CYP2D7 gene and the repeat element REP7 downstream. The first type of paralog can include the CYP2D6 gene and the CYP2D7 gene. The second type of the one or more paralogs can include the CYP2D6 / CYP2D7 fusion allele with the spacer region and the repeat element REP7 downstream of the CYP2D6 / CYP2D7 fusion allele. The allele of the first paralog can include the allele of the CYP2D6 gene that the subject has, which is a minor variant or a structural variant of the CYP2D6 gene. The computing system can determine to achieve the reference Figure 9 the allele of the CYP2D6 gene for the method 900 (or a portion thereof) described.

[0192] In various embodiments, the first paralog and the second paralog can be different. Examples of the first paralog and the second paralog include, but are not limited to, the SMN1 gene and the SMN2 gene; the CYP2D6 gene and the CYP2D7 gene; the double homeobox 4 (DUX4) gene, the DUX4c gene, the double homeobox 4-like 2 (DUX4L2) gene, the double homeobox 4-like 3 (DUX4L3) gene, the double homeobox 4-like 4 (DUX4L4) gene, the double homeobox 4-like 5 (DUX4L5) gene, the double homeobox 4-like 6 (DUX4L6) gene, the double homeobox 4-like 7 (DUX4L7) gene, and the double homeobox 2 (DUX2) gene; and the ribosomal protein S17 (RpS17) gene and the RpS17-like (RpS17L) gene. In some embodiments, the computing system can determine the copy number or alleles of the first paralog that implement the method 800 (or a portion thereof) and / or the method 900 (or a portion thereof) referred to in Figure 8 and / or the method 900 (or a portion thereof) referred to in Figure 9 the copy number or alleles of the first paralog.

[0193] In some embodiments, the first paralog and the second paralog have or have approximately 80%, 81%, 82%, 83%, 84%, 85%, 86%, 87%, 88%, 89%, 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, 99% or higher sequence identity. For example, the first paralog and the second paralog have at least 90% sequence identity.

[0194] Method 1000 ends at block 1024.

[0195] Execution environment

[0196] In Figure 11 is depicted the general architecture of an example computing device 1100 configured for paralog genotyping. Figure 11 The depicted general architecture of the computing device 1100 includes an arrangement of computer hardware and software components. The computing device 1100 can include more than Figure 11Those additional (or fewer) elements shown. However, to provide an enabling disclosure, it is not necessary to show all such generally conventional elements. As shown, computing device 1100 includes a processing unit 1110, a network interface 1120, a computer-readable media drive 1130, an input / output device interface 1140, a display 1150, and an input device 1160, all of which may communicate with each other via a communication bus. The network interface 1120 may provide a connection to one or more networks or computing systems. Thus, the processing unit 1110 may receive information and instructions from other computing systems or services via the network. The processing unit 1110 may also communicate with a memory 1170 via the input / output device interface 1140 and also provide output information for an optional display 1150. The input / output device interface 1140 may also accept input from an optional input device 1160 (such as a keyboard, mouse, digital pen, microphone, touch screen, gesture recognition system, voice recognition system, game pad, accelerometer, gyroscope, or other input device).

[0197] Memory 1170 may contain computer program instructions (grouped as modules or components in some embodiments) for the processing unit 1110 to execute to implement one or more embodiments. Memory 1170 generally includes RAM, ROM, and / or other persistent, auxiliary, or non-transitory computer-readable media. Memory 1170 may store an operating system 1172 that provides computer program instructions for use by the processing unit 1110 in the general management and operation of the computing device 1100. Memory 1170 may also include computer program instructions and other information for implementing aspects of the present disclosure.

[0198] For example, in one embodiment, memory 1170 includes a paralog genotyping module 1174 for genotyping one or more paralogs using sequencing data, such as the method 1000 described in the reference Figure 10 Alternatively or in addition, the paralog genotyping module 1174 may be or may include a module for determining the SMN1 copy number using sequencing data, such as the method 800 described in the reference Figure 8 Alternatively or in addition, the paralog genotyping module 1174 may be or may include a module for genotyping the CYP2D6 gene using sequencing data, such as the method 900 described in the reference Figure 9 In addition, memory 1170 may include a data repository 1190 and / or communicate with one or more other data repositories that store sequencing data and / or the results of genotyping one or more paralogs.

[0199] Examples

[0200] Some aspects of the embodiments discussed above are further disclosed in detail in the following examples, which are not intended to limit the scope of the present disclosure in any way.

[0201] Example 1

[0202] Spinal muscular atrophy diagnosis and carrier screening from whole-genome sequencing data

[0203] Spinal muscular atrophy (SMA), caused by loss-of-function of the SMN1 gene but retention of the paralogous SMN2 gene, is a leading genetic cause of early childhood death. Due to the nearly identical sequences of SMN1 and its paralog SMN2, analysis of this region by next-generation sequencing (NGS)-based pipelines has been challenging. The American College of Medical Genetics recommends preconception population-wide SMA screening of potential parents to quantify the copy number (CN) of SMN1.

[0204] This example describes a bioinformatics method for accurately identifying the CN of SMN1 and SMN2 using whole-genome sequencing (WGS) data. The method uses the read depth between SMN1 and SMN2 and eight informative reference genome differences to calculate the CN of SMN1 and SMN2.

[0205] Characterize the status of SMN1 / 2 in a total of 12,747 short-read whole genomes sequenced to high depth (>30x) in five ethnic groups. Among these samples, a total of 251 (1317) samples with complete gene loss (gain) of SMN1 and 6241 (374) samples with loss (gain) of SMN2 were identified. A pan-ethnic carrier frequency of 2% was calculated, which is consistent with previous studies. Additionally, the CN calls were validated, and all (48 / 48) SMN1 and 98% (47 / 48) SMN2 CN calls were consistent with those measured by digital PCR.

[0206] This WGS-based SMN copy number calling method can be used to identify carriers and affected SMA status, enabling the provision of SMA testing as a comprehensive test in neonatal care and also providing an accurate screening tool for carrier status in large-scale WGS sequencing projects.

[0207] Introduction

[0208] With the latest advancements in next-generation sequencing (NGS), it is now possible to profile large numbers of genes or even entire genomes at high throughput within clinically relevant timeframes. Driven by these advancements, large-scale population sequencing efforts are underway in many countries, and the testing for rare genetic diseases, including carrier status, will be one of the major driving factors. Spinal muscular atrophy (SMA) is an autosomal recessive neuromuscular disease characterized by the loss of α motor neurons, leading to severe muscle weakness and atrophy at birth or shortly after birth. SMA is the leading genetic cause of death in infants after cystic fibrosis. The incidence of SMA is 1 in every 6,000 to 10,000 live births, and the carrier frequency in different ethnic groups is 1:40 to 80. SMA is classified into four clinical types based on the age of onset and disease severity: infants who are too weak to sit unaided (type I), those who can sit weakly but cannot stand (type II), ambulatory patients with weaker legs than arms (type III), and a rather benign adult-onset SMA (type IV). Since two early treatments (Nusinersen and Zolgensma) that have received FDA approval for improving SMA are available, early detection of SMA can be crucial for long-term quality of life.

[0209] The SMN region includes two paralogous genes: SMN1 and SMN2. SMN2 is located 875 kb from SMN1 on 5q and was formed by a duplication of the ancestral gene unique to the human lineage. The genomic region around SMN1 / 2 undergoes unequal crossing-over and gene conversion, resulting in variable copy numbers (CN) of SMN1 and SMN2. SMN2 has greater than 99.9% sequence identity with SMN1, and one of the base differences (c.840C>T in exon 7) has important functional consequences. By disrupting a splicing enhancer, c.840T promotes the skipping of exon 7, resulting in the instability and incomplete functionality of the vast majority of SMN2-derived transcripts (70%-85%, depending on the tissue). Approximately 95% of SMA cases result from the biallelic loss of the functional c.840C nucleotide caused by the deletion of SMN1 or gene conversion (c.840T) to SMN2. In the remaining 5% of SMA cases, patients have other pathogenic variants in trans to the c.840C-deleted allele in SMN1. SMN2 can produce small amounts of functional protein, and the copy number of SMN2 in an individual modifies disease severity and is highly correlated with the clinical types described above.

[0210] Due to the high incidence and disease severity of SMA, the American College of Medical Genetics recommends population-wide SMA screening. The utility of population-wide carrier screening has been demonstrated in pilot studies. Screening for SMA includes: 1) determining the copy number of SMN1 for SMA diagnosis and carrier testing, and 2) determining the copy number of SMN2 for clinical classification and prognosis. Traditionally, SMA testing and carrier testing have been performed using polymerase chain reaction (PCR)-based assays such as quantitative PCR (qPCR), multiplex ligation-dependent probe amplification (MLPA), and digital PCR. These methods primarily determine the copy number of SMN1 based on the different c.840C>T locus between SMN1 and SMN2. This example demonstrates that WGS can achieve or exceed the performance of these tests and suggests that current and future precision medicine research initiatives can utilize genomic data for population-level screening.

[0211] Due to the nearly complete sequence identity between SMN1 and SMN2, replicating current SMA testing protocols poses problems for high-throughput WGS. Additionally, frequent gene conversion between SMN1 and SMN2 is thought to result in the generation of hybrid genes. These challenges require bioinformatics methods that can overcome the difficulties in this region. Two NGS-based tests for SMA carrier detection have been reported. Larson et al. (Validation of a high resolution NGS method for detecting spinal muscular atrophy carriers among phase 3 participants in the 1000 Genomes Project, BMC Med Genet., 2015, Vol. 16: p. 100, doi: 10.1186 / s12881-015-0246-2) used a Bayesian hierarchical model to calculate the probability that the fraction of SMN1-derived reads at three-base differences between SMN1 and SMN2 is equal to or less than 1 / 3. The method disclosed by Larson can test for SMA; however, since this method does not perform copy number calling, it is not an ideal solution for screening carriers. Instead, Feng et al. (The next generation of population-based spinal muscular atrophy carrier screening: comprehensive pan-ethnic SMN1 copy-number and sequence variant analysis by massively parallel sequencing, Genet Med Off J Am Coll Med Genet., 2017, Vol. 19 No. 8: pp. 936-944, doi: 10.1038 / gim.2016.215) described a copy number caller for both SMN1 and SMN2 based on targeted sequencing data that closely mimics current qPCR methods. Feng's method was designed for targeted sequencing and thus requires specialized normalization, which limits the method to one assay at one locus. This method derives the total copy number of SMN (including both SMN1 and SMN2) from the read coverage in exon 7 and calculates the SMN1:SMN2 ratio based on the number of reads supporting SMN1 and SMN2 at the c.840C>T locus. Using the total coverage and the SMN1:SMN2 ratio, this method derives the absolute copy numbers of SMN1 and SMN2.Since this method relies on only a single locus, it is unreliable for WGS data where the depth variability for each locus can be very high.

[0212] Compared to targeted sequencing, WGS provides a more uniform coverage across the genome and offers a less biased method for detecting copy number variants (CNVs). Additionally, WGS provides an opportunity to comprehensively analyze the population variant spectrum in the SMN region, which is poorly understood at the sequence level. This example describes a novel method for detecting the CNs of both SMN1 and SMN2 using WGS data. While most conventional assays only test for the deletion of c.840C as representative of the deletion of "exon 7 deletion", this example describes a method that can more comprehensively characterize the variability of the region, including: 1) DNA deletions, including whole gene deletions / duplications and partial deletions of the region encompassing exons 7 and 8; and 2) small variant detection, including the g.27134T>G SNP associated with SMA "silent" carriers (two copies of SMN1 on the same haplotype). The accuracy of this method was demonstrated by comparing the calls using digital PCR with the WGS-based calls of the example. A 100% (48 / 48) concordance for SMN1 and 98% (47 / 48) concordance for SMN2 were shown. Additionally, this method was applied to 2,504 unrelated samples from the 1000 Genomes Project and 10,243 unrelated samples from the NIHR BioResource Project to report the population distribution of SMN1 and SMN2 copy numbers. The carrier frequency of SMA determined using the method of this example was consistent with the carrier frequency of SMA reported in previous PCR-based studies. In addition to demonstrating the accuracy of this method for quantifying variants in the SMN region, this example also highlights the importance of using diverse ethnic groups when developing novel informatics methods to address clinically relevant difficult regions of the genome.

[0213] Materials and methods

[0214] Samples and data processing

[0215] Samples validated using digital PCR were collected from the Neuromuscular Disease Research Laboratory (Nemours Alfred I. duPont Hospital for Children) and generated from cell lines as previously described. The cohort contained 29 SMA samples (14 type I SMA, 1 type I / II SMA, 10 type II SMA, 3 type III SMA, and 1 SMA with unknown clinical grade), 6 non-SMA neuromuscular disease samples (including hereditary sensory and autonomic neuropathy type 3, myotonic dystrophy type I, distal hereditary motor neuropathy type I, and Charcot-Marie-Tooth peripheral neuropathy type IA), and 13 normal samples. SMA testing and carrier testing were performed using TruSeq DNA PCR-free sample preparation, where 150 bp paired-end reads were sequenced on an Illumina (San Diego, CA) HiSeq X instrument. Read alignment was performed using the genomic build GRCh37.

[0216] For the population study, 13,343 individuals were queried from the NIHR BioResource Rare Diseases project (EGAS00001001012), which performed WGS on individuals with rare diseases and their close relatives. An additional cohort of individuals (n = 840) from the Next Generation Children's project (EGAD00001004357) was also studied, which performed diagnostic trio WGS on patients from neonatal and pediatric intensive care units in the UK and their parents. WGS in these studies was performed using the Illumina TruSeq DNA PCR-free sample preparation kit, where 100 bp or 125 bp paired-end reads were sequenced on an Illumina HiSeq 2500, or 150 bp paired-end reads were sequenced on a HiSeq X instrument. Read alignment was performed using the genomic build GRCh37. When performing population analysis, related individuals and individuals of unknown origin were excluded, leaving 10,243 unrelated individuals.

[0217] For 1000 Genomes Project (1kGP) data, WGS BAM files were downloaded from ncbi.nlm.nih.gov / bioproject / PRJEB31736 / . These BAM files were generated by sequencing 2 × 150 bp reads from a PCR-free library on an Illumina NovaSeq 6000 instrument with an average sequencing depth of at least 30X and aligning them to the human reference sequence hs38DH using BWA-MEM v0.7.15 (with an average genome coverage of greater than 30X).

[0218] SMN copy number analysis by orthogonal methods

[0219] For validation samples, SMN1 and SMN2 CN were measured using allele-specific exon 7 probes on a QuantStudio 3D digital PCR system (Life Technologies, Carlsbad, CA) as previously described. SMN1 and SMN2 copy numbers were normalized to the copy number of RPPH1 (RNase P). Standard MLPA (SALSA MLPA P060 SMA carrier probe mix, MRC-Holland) was used to confirm SMA samples detected in the Next Generation Children's Project.

[0220] Copy number calls for full-length and truncated SMN

[0221] The SMN1 and SMN2 loci are affected by two common CNVs, namely whole-gene CNV and partial gene deletions of exons 7 and 8 (see the results of this example). The truncated form of SMN with partial deletions of exons 7 and 8 was named SMN*. This method calls the copy numbers of the full SMN1+SMN2 (hereinafter referred to as SMN) and truncated SMN (SMN*) genes using the following steps.

[0222] Identifying from SMN1 and SMN2 Reads are mapped and counted: Read counts are calculated directly from the BAM file of the WGS alignment using all reads mapped to SMN1 or SMN2 (including reads with a mapping quality of zero). In many cases, reads will align to these regions with a mapping quality of 0 because the sequences between these two regions are identical. These two genes only share sequences with each other and not with other regions of the genome. The total SMN (SMN1, SMN2, and SMN*) CN is calculated using the read counts in the 22.2 kb region including exons 1 to 6, and the CN of full SMN (SMN1 and SMN2) is calculated using the read counts in the 6 kb region including exons 7 and 8.

[0223] Calculating the normalized depth of the SMN region The read counts for the above two regions are each normalized by the region length and further normalized by dividing by the median depth of 3000 preselected 2 kb regions across the genome.

[0224] Converting the normalized depth to copy number The normalized depth values across populations are modeled using a one-dimensional mixture of 11 Gaussians, where the constrained means centered on each integer copy number value represent the copy number state in the range from 0 to 10. The copy numbers of total SMN and full SMN are called from the Gaussian mixture model (GMM) with a posterior probability threshold of 0.95.

[0225] Calculating the CN of full-length and truncated SMN:Define the full SMN CN as the CN of the 6.3 kb region covering exons 7 and 8. The copy number of truncated SMN (SMN*) was derived by subtracting the full SMN CN from the total SMN CN calculated from the 22.2 kb region including exons 1 - 6.

[0226] Genotyping the copy number of alleles at individual bases

[0227] Call the number of chromosomes carrying SMN1 and SMN2 bases by combining the total SMN CN with the read counts supporting each base in the gene - specific bases. Based on the called copy number of full SMN at each position, the method iterates through all possible combinations of SMN1 and SMN2 copy numbers and derives the combination that produces the highest posterior probability for the observed number of reads supporting SMN1 and SMN2. In addition to calling the CN of bases specific to SMN1 or SMN2, this method can also be applied to variant positions to identify the copy number of SNPs known to be specific to one of the two genes (e.g., g27134T>G), as described below.

[0228] Copy numbers of SMN1 and SMN2

[0229] For 16 positions (localized in intron 6 to exon 8) that differ between SMN1 and SMN2 in the reference genome, test whether these sites are truly fixed in the population by comparing the CN call of the SMN1 allele at these positions with the CN call of the splice - variant base SMN1 c.840C. Based on the consistency with the splice - variant base, eight positions where the SMN1 base is fixed or nearly fixed in the population were identified, including c.840C>T (see the results section of this example, Figure 14A ). The remaining sites may be polymorphic in the population and may not be reliable for CN calling.

[0230] To make a final CN call, the method requires: 1) the SMN1 CN calls to be consistent across at least 5 out of 8 sites with a posterior probability cutoff of 0.8, or 2) at least 5 out of 8 sites (posterior probability > 0.6) to be consistent with the CN call derived from reads overlapping all 8 sites (posterior probability > 0.9). Otherwise, no call is made for both SMN1 CN and SMN2 CN. SMA samples are identified as having zero copies of full SMN1, and carrier samples are identified as having one copy of full SMN1.

[0231] At higher CN values, the read depth is expected to have greater variability, resulting in lower confidence (lower posterior probability) in CN calls for individual loci and greater inconsistency between loci. Thus, non-calls are more likely to occur in samples with high SMN1 / SMN2 CN (i.e., both values greater than or equal to 2) (see Figure 15 ). However, it is still possible to confidently determine whether the SMN1 CN is 0 (SMA) or 1 (carrier) in such samples, which allows for the calling of SMA / non-SMA or carrier / non-carrier. When the SMN1 copy number is a non-call, the sample is called "non-SMA" if at least seven of the SMN1 CN calls are confidently greater than zero. Similarly, the sample is called "non-carrier" if at least seven of the SMN1 CN calls are confidently greater than 1. Additionally, when the SMN1 CN is a non-call, the deletion of the c.840C allele, which indicates SMA, is directly tested. This is done by testing whether the number of reads supporting the SMN1 base (c.840C) is more likely to be derived from zero or one copy of SMN1.

[0232] Results

[0233] Common CNVs affecting the SMN1 / SMN2 locus

[0234] The genes SMN1 and SMN2 are located in a ~2 Mb region of the reference genome with a large number of complex segmental and inverted segmental duplications. While existing methods (e.g., PCR-based methods) mainly focus on the c.840C>T locus, this example illustrates a copy number method based on sequencing data from the whole gene. The SMN1 copy number is defined as the number of SMN genes carrying the c.840C allele, and the SMN2 copy number is defined as the number of SMN genes with the c.840T allele. Sequence analysis was performed using high-depth (>30x) WGS data from 2,504 samples of the 1000 Genomes Project (1kGP) and 10,243 unrelated samples from the NIHR BioResource project (see the methods of this example).

[0235] To develop a CN calling strategy, two common CNVs that result in DNA deletions were first characterized. The primary CNV evaluated involved the entire SMN1 / SMN2 gene region. The read depth across a ~30 kb homologous region encompassing the SMN1 and SMN2 genes was examined. Figure 12AShows the normalized read depth in a 100 bp sliding window in samples with different SMN1 + SMN2 CNs across the region (representing both SMN1 and SMN2). The depth profiles indicate that the entire region is either deleted or duplicated in these samples. Due to extensive sequence homology both within and outside this region, the exact breakpoints of this CNV are expected to vary among samples and can only be resolved with high resolution using long-read sequencing. For SMA testing, the analysis was limited to the region (∼30 kb) containing the SMN gene (SMN1 or SMN2).

[0236] In addition to whole-gene CNVs, a 6.3 kb partial gene deletion encompassing both exons 7 and 8 was also detected ( Figure 12B 、 Figure 16 ). The sequences at the breakpoints are identical between SMN1 and SMN2, so the deletion occurred at chr5:70244114 - 70250420 in SMN1 or chr5:69368689 - 69375000 in SMN2 ( Figure 16 , hg19). However, approximately 500 bp downstream of the breakpoint defining the end of the deletion, there are three base differences between the SMN1 and SMN2 loci (70250881A>69375425C, 70250981A>69375525G, 70250991A>69375535G). In samples containing this deletion, 245 read pairs from 237 samples were identified, where one read spanned the breakpoint and the other read spanned at least two of the three SMN-discriminating bases. Analysis of these read pairs indicated that 100% were consistent with a deletion occurring on the SMN2 sequence background. This truncated form of SMN2 was named "SMN*", and since both exon 7 and exon 8 are deleted, SMN* most likely has limited or no biological function. Therefore, SMN* is an important variant that any SMN CN caller should consider.

[0237] Figure 12A and Figure 12B Shows non-limiting exemplary curves illustrating common CNVs affecting the SMN1 / SMN2 loci. Figure 12A Shows depth profiles across the SMN1 / SMN2 region. Samples with 2, 3, 4, and 5 total SMN1 + SMN2 copy numbers are shown as dots respectively. For each CN category, the depths of 50 samples were summed. Each point represents the normalized depth value in a 100 bp window. Read counts were calculated in each 100 bp window, the reads for both SMN1 and SMN2 were summed, and normalized to the depth of wild-type samples (CN = 4). The SMN exons are represented as purple boxes. The two x-axes show the coordinates in SMN1 (bottom) and SMN2 (top).Figure 12B A depth profile aggregated from 50 samples carrying exon 7 and 8 deletions is shown as dots. Read depth is calculated in the same manner as in Figure 12A .

[0238] After searching for abnormal read pairs, no other common CNVs were found in the SMN region. Combining this information, CNs of the SMN gene were called by dividing the gene into two regions (a 6.3 kb region containing exons 7-8 and a 22.2 kb region containing exons 1-6) to specifically identify the numbers of the full and truncated forms. The CNs of these two regions were calculated based on the read depth as described in the method section of this example. The calculated CN from the exon 7-8 region provides the number of full SMN genes. Samples with SMN* have a higher CN call from the exon 1-6 region compared to the CN call from the exon 7-8 region, and the difference between them represents the CN of SMN*. Figure 13 The calculation results for a cohort of 12,747 samples are shown, where 2,144 cases of SMN* were identified, including 140 samples with two copies of SMN* and one sample with three copies of SMN*.

[0239] Figure 13 A non-limiting exemplary scatter plot of the total SMN (SMN1 + SMN2) copy number (x-axis, called by read depth in exons 1 to 6) and the full SMN copy number (y-axis, called by read depth in exons 7 to 8) is shown.

[0240] Distinguishing SMN1 CN from SMN2 CN

[0241] After calculating the total copy number of the SMN gene, SMN1 and SMN2 were distinguished as follows. Since c.840C>T is the most important functional difference between SMN1 and SMN2, the absolute copy numbers of these two genes can theoretically be derived using the ratio between the numbers of reads supporting SMN1 and SMN2 at this locus. However, for WGS datasets, the read depth at a diploid position is usually 30 - 40X, and sometimes it does not provide sufficient power to clearly distinguish different CN states (see Figure 15 ). Therefore, additional base differences near c.840C>T were utilized so that when making CN calls, the information at these sites can be combined with c.840C>T. Since it is necessary to distinguish full SMN1 from SMN2, variants occurring within the 6.3 kb deletion were considered. SNPs and short tandem repeats (STRs) in homopolymers, which are more likely to be error-prone, were excluded, resulting in 16 base differences between SMN1 and SMN2 (Table 8).

[0242] For these 16 base differences, the CNs of SMN1 and SMN2 alleles were called independently (see the method section of this example), and the CN call at each position was compared with the CN call at the splice variant site ( Figure 14A , Figure 17 ). There was a significant difference in the call consistency between African and non - African populations ( Figure 14A ). For non - African samples, there were 13 sites with high (>85%) CN consistency with the splice site. In contrast, for African samples, only seven sites had high CN consistency with the splice site, and the consistency values were lower than those in non - African populations. This is consistent with the intragenic variation at many of these positions and the higher frequency of these non - reference alleles in African populations. Splice variants and seven positions with high consistency with the splice variant were selected in African and non - African populations to make CN calls for SMN1 and SMN2. By restricting to two CN states (SMN1 = CN2 and SMN2 = CN0 or SMN1 = CN2 and SMN2 = CN1) that allow easy identification of hybrid alleles, it was possible to estimate the allele frequencies at these sites on the SMN1 and SMN2 genes (Table 9, Figure 18A and Figure 18B ). Based on this analysis, in these eight positions, at most 0.5% of the SMN1 gene was estimated to contain SMN2 alleles. In contrast, at most 0.9% of the SMN2 gene was estimated to carry SMN1 alleles. These observations may be the result of gene conversion, or many of these eight sites are polymorphic in the population. Most of these hybrid alleles are from African populations (Table 9).

[0243] Figures 14A to 14D shows the distribution of SMN1 / SMN2 / SMN* copy numbers in the population. Figure 14A is a non - restrictive exemplary figure that shows the percentage of samples showing CN call consistency with c.840C>T at 16 SMN1 - SMN2 base - difference sites in African and non - African populations. Site 13* is the c.840C>T splice variant site. The black horizontal line represents 85% consistency. Figure 14B shows a non - restrictive exemplary bar graph of the SMN1, SMN2, and SMN* copy number distributions in five populations in the 1kGP and NIHR BioResource cohorts (values are shown in Table 15). Figure 14C is a non - restrictive exemplary curve of SMN1 CN versus total SMN2 CN (complete SMN2+SMN*). Figure 14DTwo trios are shown where the SMA proband was detected by the caller and orthogonally confirmed in the NIHR BioResource cohort. The CN of each allele of SMN1, SMN2, and SMN* was phased and each member of the trio was genotyped.

[0244] Introducing more base differences improved the ability to distinguish SMN1 from SMN2. However, because these sites are not truly invariant in their respective genes and CN calls at individual sites can be in error, the likelihood that a call in an individual call deviates from the true copy number state increases. To make a final call, the SMN1 CN calls need to agree with each other at 5 or more of the 8 sites (for a full description of the CN calling rules, see the Methods section of this example). At a posterior probability cutoff of 0.8, most samples had agreeing calls at at least 5 of the eight sites, and only 1.4% of samples had fewer than 5 agreeing sites (Table 10). In 80% of these samples, a confident CN call was made based on a second consensus rule (which needed to agree with the CN call made by summing the reads at all 8 sites). Due to lower posterior probabilities rather than discrepant calls, the "disagreeing" sites were more often non-calls, and only 15.3% of these sites were confident calls that disagreed with the consensus at other sites. Similarly, most of the disagreements came from the African population (Table 10). Using fewer sites for the majority rule resulted in a greater number of non-calls and incorrect calls compared to using eight sites (Table 11).

[0245] Validation of the SMN copy number caller

[0246] To test the method, 48 samples with known SMN1 and SMN2 CN (including 29 SMA probands, 6 SMA carriers, and 13 samples with SMN1 CN greater than 1) were sequenced. In all 48 cases, the SMN1 CN calls agreed with the digital PCR results, and in 47 of the 48 cases, the SMN2 CN calls agreed (97.9%) (Tables 6A and 6B). In this single discrepant case (MB509), the method called an SMN2 CN of 3, while digital PCR showed an SMN2 CN of 2 (Table 12). Upon careful examination, a 1884bp deletion in SMN1 in this sample was found (chr5:70247145-70249029, hg19) ( Figure 19)。The deletion is small (it does not significantly change the depth in the 6 kb region used to determine the full SMN CN) and has not been reported previously (nor found in population data), so the method was not designed to detect it. Therefore, the sample was correctly identified as SMA, but the SMN2 CN was overestimated by one. The deletion is consistent with the CN calls made at the 8 SMN1 - SMN2 differential sites, where the first 2 sites are not in the deletion and are called at SMN1 CN = 1, and the last 6 sites are in the deletion and are called at SMN1 CN = 0.

[0247] The concordance of SMN1 / SMN2 / SMN* CN calls was analyzed in 258 trios from the Next Generation Children's Project cohort (see the Methods section of this Example). There were no Mendelian errors in any of the calls (Table 13).

[0248] Table 6A. Validation performed on samples with known SMN1 / SMN2 CN 。

[0249]

[0250] Table 6B. Validation on samples with known SMN1 / SMN2 copy number (CN) 。

[0251]

[0252] Copy numbers of SMN1, SMN2, and SMN* by population

[0253] Given the high accuracy demonstrated by validating the digital PCR results, the method was applied to the high - depth (>30x) WGS data of 12,747 unrelated samples from the 1kGP and NIHR BioResource projects (Table 14). The CN distributions were analyzed by population (Europeans, Africans, East Asians, South Asians, and Admixed Americans consisting of Caucasians, Mexican - Americans, Peruvians, and Puerto Ricans). Figure 14B A bar graph showing the number of individuals with different CNs of full SMN1, full SMN2, and SMN* is presented. The distributions between the 1kGP samples and the NIHR BioResource samples are similar ( Figure 20 ). Generally, individuals have more copies of SMN1 than SMN2. The most common combinations of SMN1 / SMN2 copy numbers are 2 / 2 (44.9%) and 2 / 1 (33.4%). Except for Africans, who exhibit higher variability in SMN1 and SMN2 CNs, the variability in SMN1 copy number is much lower than that in SMN2 copy number. In contrast, 54.7% of Africans have three or more copies of SMN1, which is more than twice the copy number observed in any of the other four populations ( Figure 14B 、Table 7). There is an inverse relationship between the copy number of SMN1 and the copy number of SMN2, where the CN of SMN2 decreases as the CN of SMN1 increases (Figure 14C , correlation coefficient -0.344, p-value < 2.2e-16). This observation is consistent with the mechanism of gene conversion occurring between SMN1 and SMN2. The observed higher SMN1 CN relative to SMN2 CN may be the result of SMN2-to-SMN1 conversion or bias in selection against low SMN1 CN. Africans have significantly lower SMN2 CN than other groups.

[0254] The number of SMA carriers identified in the populations is summarized in Tables 7 and 15. Among 12,683 individuals with confident SMN1 / SMN2 CN calls, Europeans have the highest carrier frequency of 2.2%, followed by admixed Americans (2.05%), East Asians (1.35%), and South Asians (1.67%). Africans have the lowest carrier frequency (0.44%). The CN frequency distribution observed in this example is consistent with previous studies of the SMN1 / SMN2 CN distribution in the general population. In addition, the frequency of exon 7-8 deletions (SMN*) was determined in the populations: 21.2% of Europeans and 11.5% of admixed Americans have at least one copy of SMN*, while the frequencies are lower for South Asians (3.35%), Africans (1.1%), and East Asians (0.34%).

[0255] In the Next Generation Children's Project cohort (see the Methods section of this example), SMA was identified in two neonatal probands from trio analysis with independent confirmation. In addition, the phasing of SMN1 CN, SMN2 CN, and SMN* CN was performed for each trio member ( Figure 14D ).

[0256] Simulation of single-site CN calling

[0257] Based on the Poisson distribution, the number of reads at a single locus with median depths of 30X, 35X, and 40X in the sample was simulated, and reads supporting SMN1 were sampled based on the binomial model for all possible combinations of SMN1 CN and SMN2 CN, where the total SMN CN is between 2 and 6. The posterior probability of the simulated SMN1 CN was derived based on the number of reads supporting SMN1 and SMN2 (see the Methods section of this example). When at least one of the SMN1 CN or SMN2 CN values is low (less than or equal to 1), the posterior probability is high (greater than 0.9) ( Figure 16)。When both values are greater than 2, i.e., in the SMN1:SMN2 combinations of 2:2, 2:3, 2:4, 3:2, 3:3, and 4:2, the posterior probability tends to be lower and drops below 0.9. This is because when the expected CN is higher, the variability of the read depth is greater. Therefore, in these scenarios, using a single locus to make SMN1 and SMN2 CN calls may not be accurate enough.

[0258] Differences in validation samples

[0259] There is a sample MB509 that shows a discrepancy between our CN call and the digital PCR result. Upon further examination, it was found that this sample has two copies of SMN2 and one copy of SMN1, with an 1884 bp deletion (chr5: 70247145 - 70249029, hg19, Figure 20 )。Although the read alignment in the SMN1 / 2 region is not always accurate, careful analysis of the split reads indicates that these reads or their mates overlap with bases specific to SMN1. Without intending to be bound by theory, it is assumed that the deletion is correctly located on SMN1. The deletion is small (it does not significantly change the depth in the 6.3 kb region used to determine the full SMN CN) and has not been reported previously (nor found in the 1kGP samples, so it is a very rare variant), and thus the method was not designed to detect deletions. Therefore, the total copy number of SMN1 + SMN2 called by this method is 3. The deletion is consistent with the CN calls made at 8 SMN1 - SMN2 differential sites, where the first 2 sites are not in the deletion and are called at SMN1 CN = 1, and the last 6 sites are in the deletion and are called at SMN1 CN = 0 ( Figure 21A )。Based on the majority rule, the SMN1 copy number called by this method is 0, correctly identifying the sample as SMA. The SMN2 copy number is calculated as the total copy number minus the SMN1 copy number, so the SMN2 copy number called by this method is 3, overestimating by 1.

[0260] Four other samples, MB231, MB367, MB383, and LP2101748, have discrepancies between the CN calls made and the results from digital PCR or MLPA. The read counts and normalized depth values (read counts divided by the haploid sample depth) at 8 base - difference sites support our CN calls ( Figure 21A ),and the discrepancies may be caused by errors in the orthogonal methods. In two samples, the genomic sequencing (GS) calls and the digital PCR calls differ by a factor of two (MB231: GS - 0,2, PCR - 0,4 and MB383: GS - 3,1, PCR - 6,2). There may be normalization issues in digital PCR, resulting in an overestimation of the copy number by a factor of two.

[0261] When comparing the CN calls made in 1109 1kGP samples with the MLPA results, one sample with a no call for SMN2Δ7-8 due to a low posterior probability of the total SMN CN was excluded, as well as three samples with no calls for SMN1 and SMN2 CN due to inconsistent CN calls at 8 selected sites that did not meet the common sequence rule. Figure 21B )

[0262] Detection of "silent" carriers

[0263] The g.27134T>G SNP may be associated with the 2+0 SMA silent carrier state, where one chromosome carries two copies of SMN1 (either by SMN1 duplication or gene conversion of SMN2 to SMN1), while the other chromosome has no copy of SMN1. The method of this example can also detect the presence of this SNP and thus can be used to screen potential silent carriers. This SNP is closely associated with the double-copy SMN1 allele in Africans, where 84.5% of individuals with three copies of SMN1 and 92.6% of individuals with four copies of SMN1 have the g.27134T>G SNP (Table 7). Calling this SNP greatly improves the carrier detection rate in Africans because Africans have a higher frequency of alleles carrying two copies of SMN1 (Tables 17 and 18). However, 33% of individuals with two copies of SMN1 also have the g.27134T>G SNP, indicating that a large portion of the single-copy SMN1 alleles also carry this SNP. Calculate the maximum likelihood estimates of the percentages of single- and double-copy SMN1 alleles carrying g.27134T>G (Table 17) and the maximum likelihood estimates of the residual risks of the combination of CN and SNP calls (Table 18). The calculated estimates are similar to previous studies, but there is considerable variability among all of these estimates. This variability may be caused by population variability, such as Africans (this example) compared to African Americans (previous studies) and Northern Europeans (overrepresented in this example) relative to the more diverse Caucasians sampled (previous studies).

[0264] Table 7. SMN1 CN and g.27134T>G frequency by population 。

[0265]

[0266] Comparison between two aligners, BWA and Isaac

[0267] The method of this embodiment allows for the analysis of reads in both SMN1 and SMN2, and thus is insensitive to how the aligner differentiates between the two genes. Therefore, using different aligners should produce similar results. The BAM data analyzed in this embodiment was generated using two different aligners: BWA for 1kGP data and various versions of Isaac for the remaining data. The consistent SMN1 / 2 CN distributions between 1kGP and NIHR (Table 19, Figure 20 ) samples indicate that our method is insensitive to the aligner. Additionally, the consistency of the method was tested by aligning 117 samples (including 5 SMA samples and 3 carriers) with both BMA and Isaac. Using the method of this embodiment, all 117 samples produced exactly the same calls (SMN1 / SMN2 / SMN2Δ7-8 CN), and the normalized depths for both exons 1-6 and exons 7-8 were nearly identical (Pearson's r > 0.999, Figure 22 ).

[0268] Comparison between carrier calls in this study and those by Larson et al.

[0269] Comparing the carrier calls (N = 37) made in the 1kGP samples of this embodiment with the carrier calls (N = 36) reported by Larson et al., 26 overlapping calls were found (Table 15). Assuming the calls made by the method of this example are correct, Larson et al. made 10 false positive (FP) and 11 false negative (FN) calls. Larson et al. identified carriers by determining whether the fraction of reads supporting SMN1 was less than or equal to 1 / 3. This study used low-depth sequencing data that was expected to result in some errors, but more importantly, their method is error-prone without calling the total copy number. For example, a sample with one copy of SMN1 and one copy of SMN2 would be called a non-carrier (SMN1 fraction of 1 / 2), and a sample with two copies of SMN1 and four copies of SMN2 would be called a carrier (SMN1 fraction of 1 / 3), resulting in false positives and false negatives (Table 16).

[0270] Additional figures and tables

[0271] Figure 15 Non-limiting exemplary curves are shown, each curve showing the posterior probability distribution of simulated SMN1 CN using a single locus at different read depths and SMN1:SMN2 CN combinations.

[0272] Figure 16Shows a non - restrictive exemplary IGV snapshot of the SMN2 region in a sample with exon 7 - 8 deletion. The horizontal lines connect two reads in pairs in the center - aligned track. The BLAT results of two split reads spanning the breakpoint are shown in the bottom track, which shows two segments of the same read aligned to either side of the deletion breakpoint.

[0273] Figure 17 Shows a non - restrictive exemplary graph that shows the correlation between the raw SMN1 CN at a 15 - base difference near c840.C>T and the raw SMN1 CN at the c840.C>T locus. The raw SMN1 CN at each locus is calculated as the CN of full SMN multiplied by the fraction of read counts supporting SMN1 in the read counts supporting SMN1 + SMN2. The correlation coefficient is listed in the title of each graph.

[0274] Figure 18A and Figure 18B Shows a non - restrictive exemplary graph that shows the SMN1 / SMN2 haplotypes in samples with SMN1:2SMN2:0 and SMN1:2SMN2:1 in 1kGP. The y - axis shows the raw SMN1 CN as defined in Figure 16 The x - axis shows 16 loci, which are indexed and explained in Table 8. Index #13 represents the c840.C>T locus. Samples with SMN1:2SMN2:0 are shown together in the upper - left panel. Samples with SMN1:2SMN2:1 are shown as 5 clusters. Figure 18A . Non - African. Figure 18B . African.

[0275] Figure 19 Shows a non - restrictive exemplary IGV snapshot showing a 1.9 kb deletion of SMN1 in MB509.

[0276] Figure 20 Shows a non - restrictive exemplary graph that shows SMN1 / SMN2 / SMN*CN in 1kGP and the NIHR cohort.

[0277] Figure 21A and Figure 21B Shows the differences and no - calls in the validation samples. Figure 21AShows five samples with differences between GS calls and digital PCR or MLPA results. The x-axis shows 16 loci, the indices of which are listed and explained in Table 8. Index #13 represents the c840.C>T locus. The left y-axis of the bar graph shows the read counts supporting SMN1 and SMN2. The right y-axis of the bar graph shows the normalized read depths of SMN1 and SMN2 (representative of copy number, read count divided by haploid depth). The title of each subplot shows the GS and digital PCR / MLPA calls for SMN1 and SMN2 of each sample, separated by commas. Figure 21B Shows three 1kGP validation samples in which the SMN caller makes no calls for SMN1 and SMN2 CN due to inconsistencies between SMN1 / SMN2 base difference loci. The eight loci for the common sequence rules used in the method are #7-8 and #10-15. The y-axis shows the raw SMN1 CN as Figure 17 defined.

[0278] Figure 22 Shows CN calls derived from BWA and Isaac BAM.

[0279] Table 8. Genomic coordinates of base differences between SMN1 and SMN2 。

[0280]

[0281] Table 9. Frequencies of SMN1 haplotypes with SMN2 alleles and SMN2 haplotypes with SMN1 alleles in two simple CN states (SMN1 = CN2 and SMN2 = CN0 or SMN1 = CN2 and SMN2 = CN1). The numbers in parentheses represent the haplotypes contributed by the African population.

[0282]

[0283] Table 10. Number of samples with different numbers of concordant sites at 8 SNP loci. The numbers in parentheses indicate the number of samples contributed by African populations 。

[0284]

[0285] * Calls were made in these samples based on the second majority rule (see Methods).

[0286] Table 11. Number of no calls due to discordance and differential calls made using a reduced number of loci 。

[0287]

[0288] Table 12. Validation samples

[0289]

[0290]

[0291]

[0292] Table 13. SMN1, SMN2, and SMN* CN calls for 258 trios in the Next Generation Children Project cohort 。

[0293]

[0294]

[0295] Table 14. Number of samples by population in the 1kGP and NIHR BioResource cohorts 。

[0296]

[0297]

[0298] Table 16. Comparison of carrier calls made in this example and by Larson et al. in 1kGP samples 。

[0299]

[0300]

[0301] Table 17. Maximum likelihood estimates of the percentages of single and double copies of the SMN1 allele carrying g.27134T>G estimates 。

[0302]

[0303] *The NIHR BioResource cohort represents the majority of the European population analyzed in this example due to its large sample size, including Nordic samples that carry a lower frequency of the g.27134T>G SNP than the more diverse European samples from the 1000 Genomes Project.

[0304]

[0305] Table 19. SMN1 / SMN2 / SMN2Δ7-8 CN in the 1kGP and NIHR cohorts

[0306]

[0307] Discussion

[0308] Due to the high sequence homology between SMN1 and SMN2, the SMN region has been difficult to resolve using both short-read and long-read sequencing, and to date, this important region has been excluded from standard WGS analysis. This example demonstrates a method for independently resolving the copy numbers (CNs) of SMN1 and SMN2 using short-read WGS data, thus filling an important gap in SMA diagnosis and carrier screening to enable precision medicine research initiatives. Accurate measurement of SMN1 and SMN2 CN is important not only for the diagnosis of SMA but also as a prognostic indicator and basis for treatment selection. SMN2 CN has been used as a criterion in many SMA clinical trials, including Nusinersen and Zolgensma.

[0309] As proof of the method, CN calls were made for SMN1 and SMN2 using sequencing data from 12,747 samples covering five different subpopulations. The following samples were identified: 251 samples with a complete gene loss (less than two copies) of SMN1 and 1317 samples with a complete gene gain (more than two copies); 6241 samples with a complete gene loss of SMN2 and 1274 samples with a complete gene gain; 2144 samples carrying one or more copies of the truncated form SMN*. The driving role of deletions, duplications, or gene conversions in CN changes in this region could not be accurately quantified. Evidence supporting all three mechanisms included: 1) 3853 samples with a total (SMN1 + SMN2) CN < 4 (deletion), 2) 670 samples with a total CN > 4 (duplication), and 3) a strong inverse correlation between SMN1 and SMN2 CN (gene conversion, Figure 14C ). Additionally, carrier frequencies between 1:42 and 1:101 were identified according to the ancestral population (Table 7). CN frequencies varied widely by population, and the results for each population in this example were consistent with previous population studies. While this consistency provided qualitative support for the accuracy of the method, the accuracy of the method was directly evaluated by comparing the CN calls made by the method with the results of digital PCR. In this direct comparison, all (48 / 48) SMN1 and 98% (47 / 48) SMN2 CN calls made by the method were consistent with the digital PCR-based results. One inconsistency was due to a 2-kb deletion not targeted by the method, and importantly, the method correctly identified the SMA status of the sample.

[0310] In this example, the CN calls were optimized to be applicable to individuals of any ancestry and thus restricted SMN1 / 2 differentiation to functionally important splice variants highly consistent with splice variants in all populations plus seven loci ( Figure 14A) By quantifying the concordance between all the differences in the reference differences and the splice variants, this method can identify changes in these fixed differences that, if not properly accounted for (e.g., removed from our analysis), may lead to errors in our CN calls. Not accounting for fixed differences when analyzing Africans will be particularly problematic because these Africans have more diverse haplotypes. Population genetics studies (e.g., including the use of long-read sequencing) can help to more directly profile haplotype diversity between populations and identify new variant loci that can further improve the accuracy of SMN1 / SMN2 differentiation.

[0311] One type of "silent" carrier occurs when an individual has two copies of the SMN1 gene but they are both on the same haplotype. The SNP (g.24134T>G) has been used to identify individuals at increased risk of being a carrier when SMN1 CN is 2, but the risk associated with this SNP can vary widely between studies and populations (Table 17). When an individual has only one copy of SMN1, the individual can be definitively identified as a carrier, but this variant only indicates a 2% to 8% chance of being a carrier when SMN1 CN is 2. Using WGS, it is possible to classify the different variants that occur for different CN combinations of SMN1 and SMN2 and to identify additional markers that can be used to improve our ability to identify these "silent" carriers. In addition, the loss of the current c.840C>T splice variant accounts for approximately 95% of SMA cases, and the remaining cases include other pathogenic variants. These other pathogenic variants represent another type of "silent" carrier. This method can directly genotype these other pathogenic variants as part of the testing process, thus further improving the ability to detect SMA carriers and cases.

[0312] Although there are difficult regions in the genome where normal WGS pipelines do not deliver variant calls, this example demonstrates the ability to apply WGS paired with targeted bioinformatics methods to address one such difficult region. This targeted strategy (WGS + specialized bioinformatics) can be applied to many difficult variants, such as the repeat expansions and CYP2D6 disclosed herein. Traditionally, it has been cost-effective to perform all known genetic tests and carrier screenings on each individual, and thus information such as carrier rates and family history has been used to identify candidates for specific genetic tests. However, this process means that many people without a family history who could benefit from knowing their SMA status generally do not have access to this data. Once WGS analysis can accurately detect all SNVs and CNVs in all clinically relevant genes, a more universal population-wide genetic testing strategy can be achieved through a single test. Improving WGS to make it an economic alternative to current genetic testing will help to promote the integration of more genetic testing and carrier screening into WGS, thus making genetic testing more widely available to the general population. WGS provides a valuable opportunity to assess genetic variation across the entire genome, and leveraging more targeted bioinformatics solutions developed for difficult regions will help to bring the prospect of personalized medicine closer to reality.

[0313] Example 2

[0314] Accurate CYP2D6 genotyping using whole-genome sequencing data

[0315] This example and Appendix A describe the genotyping of CYP2D6 using whole-genome sequencing data. The content of Appendix A is incorporated herein by reference in its entirety.

[0316] CYP2D6 is involved in the metabolism of 25% of all drugs and is a key target for personalized medicine. Genotyping CYP2D6 is challenging due to its high polymorphism, the presence of common structural variants (SVs), and its high sequence similarity to the pseudogene paralog CYP2D7 of the gene. This article discloses bioinformatics methods that can accurately genotype CYP2D6 using whole-genome sequencing (WGS) data, also referred to herein as Cyrius. In 138 samples with GeT-RM consensus sequence calls and 50 additional samples with sequencing data from Pacific Biosciences of California, Inc. (Menlo Park, CA), also known as PacBio, this method (97.9% concordant with truth) had superior performance compared to other methods (85.6% to 88.8%). A specific differentiator of this method is the ability to call structural variant star alleles. This method correctly identified 97.2% (70 / 72) of structural variant star alleles, compared to 77.8% - 88.9% (56 / 72 and 64 / 72) of structural variant star alleles identified by other methods. Applying this method to 2504 samples from the 1000 Genomes Project (1kGP), the frequency of CYP2D6 star alleles estimated to involve SVs was 32.2% higher than previously reported for some populations. This example provides benchmarking results for the largest validation dataset to date. In some embodiments, this method is a useful tool for pharmacogenetic applications of WGS. This method can help bring the prospect of precision medicine closer to reality.

[0317] Introduction

[0318] Individuals vary significantly in their response to a large number of clinically prescribed drugs. A powerful factor contributing to this differential drug response is the genetic makeup of drug-metabolizing genes. Precision medicine requires genotyping of drug genes to enable personalized treatment. Cytochrome P450 2D6 (CYP2D6) is one of the most important drug-metabolizing genes and is involved in the metabolism of 25% of drugs. The CYP2D6 gene is highly polymorphic, with 106 star alleles (Pharmvar.org / gene / CYP2D6) defined by the Pharmacogene Variation (PharmVar) Consortium. CYP2D6 star alleles are CYP2D6 gene copies defined by a combination of small variants such as single nucleotide variants (SNVs) and insertions / deletions (indels) and structural variants (SVs), and correspond to different levels of CYP2D6 enzyme activity, such as poor, intermediate, normal, or ultra-rapid metabolizers.

[0319] Genotyping of CYP2D6 is challenged by the presence of the non-functional paralog CYP2D7, which is located upstream of CYP2D6 and has 94% sequence similarity with several nearly identical regions. Deletions and duplications of CYP2D6, as well as fusions between CYP2D6 and its pseudogene paralog CYP2D7, are common. Traditionally, CYP2D6 genotyping has been performed using array- or polymerase chain reaction (PCR)-based methods such as TaqMan assays, droplet digital PCR (ddPCR), and long-range PCR. These assays differ in the number of star alleles (variants) they target, resulting in differences in genotyping results among different assays. Common limitations of these methods are: 1) the wild-type allele *1 is usually the default call when the targeted variant is not detected, or 2) parental alleles such as *2 are assigned when variants defining true star alleles are not tested. These assays are low-throughput and generally difficult to detect structural variants.

[0320] Whole genome profiling may be possible at high throughput over clinically relevant time periods by next-generation sequencing (NGS). Large-scale population sequencing efforts have been carried out, and pharmacogenetic testing can be a desired goal. CYP2D6 genotyping using NGS is particularly challenging due to the common gene conversion between CYP2D6 and CYP2D7 (hereinafter referred to as CYP2D6 / 7), common SVs (gene deletions, duplications, and CYP2D6 / 7 fusion genes), and sequence similarity between CYP2D / 7, which results in ambiguous read mapping for these two genes. Some existing callers cannot detect complex structural variants and have been shown to have low performance. Other existing callers, such as Aldy (Numanagic et al., Allelic decomposition and exact genotyping of highly polymorphic and structurally variant genes, Nat Commun., 2018, Vol. 9, No. 1: pp. 1-11, Doi: 10.1038 / s41467-018-03273-1) and Stargazer (Lee et al., Stargazer: a software tool for calling star alleles from next-generation sequencing data using CYP2D6 as a model, GenetMed., 2019, Vol. 21, No. 2: p. 361, Doi: 10.1038 / s41436-018-0054-0), rely on exact read mapping of sequence reads to CYP2D6 in order to detect SVs based on depth and derive haplotype configurations based on observed small variants and SVs. However, accurate read mapping of sequence reads to CYP2D6 is often not possible at many positions across the gene because the sequence is highly similar to CYP2D7 or even indistinguishable due to gene conversion. Thus, the depth pattern can be ambiguous, and the caller can make false positive / negative small variant calls. Some callers do not support hg38, so many studies will need to be realigned to hg37 to use these tools.

[0321] A set of reference samples provided by the CDC Genetic Testing Reference Materials Program (GeT-RM; Gaedigk et al., Characterization of Reference Materials for Genetic Testing of CYP2D6 Alleles: A GeT-RM Collaborative Project, J Mol Diagn JMD, August 2019, Doi: 10.1016 / j.jmoldx.2019.06.007) enables the genotyping accuracy of newly developed methods to be evaluated, where the consensus genotypes of major pharmacogenetic genes are derived using multiple genotyping platforms. GeT-RM covers 43 of the 106 CYP2D6 star alleles. Additionally, many of the single-marker methods for these consensus genotypes can be error-prone, leading to conflicts between methods. High-quality long reads are available to provide a complete picture of CYP2D6 to improve the validation of complex variants and haplotypes. This paper discloses Cyrius, a WGS-based CYP2D6 genotyping method that overcomes the challenges of CYP2D6 and CYP2D7 (referred to herein as CYP2D6 / 7). Cyrius has better genotyping accuracy than Aldy and Stargazer in 138 GeT-RM reference samples and 50 samples with whole-genome PacBio sequencing data, covering 41 of the 106 known star alleles. The method was applied to high-depth sequence data from 2504 unrelated samples from the 1000 Genomes Project (1kGP) to report the distribution of star alleles in five ethnic groups. This analysis demonstrated differences in frequencies in PharmGKB, highlighting the potential errors associated with combined limited star allele calls using multiple techniques designed to identify specific subgroups of known star alleles. This analysis expands the current understanding of CYP2D6 gene diversity, particularly for complex star alleles with SVs.

[0322] Materials and methods

[0323] Samples

[0324] Analyze the following samples: Whole-genome sequencing (WGS) data of 138 GeT-RM reference samples (including 96 samples that were genotyped in the initial GeT-RM study and updated in the latest GeT-RM version) and 42 additional samples newly added in the latest GeT-RM version. For the first batch of 96 samples, WGS was performed using TruSeq DNA PCR-free sample preparation, with 150-bp paired-end reads sequenced on an Illumina, Inc. (San Diego, CA) HiSeq X instrument. Read alignment was performed using the genome build GRCh37. Sequence data of 70 of these samples were downloaded from ebi.ac.uk / ena / data / view / PRJEB19931. WGS data of the second batch of 42 samples were downloaded from the NYGC as part of the 1000 Genomes Project (see below).

[0325] For population studies, 1000 Genomes Project (1kGP) data were used, where WGS BAMs of 2,504 samples were downloaded from ncbi.nlm.nih.gov / bioproject / PRJEB31736 / . These BAM files were generated by sequencing 2 × 150-bp reads from a PCR-free library on an Illumina NovaSeq 6000 instrument and aligning them to the human reference sequence hs38D. WGS data of 70 GeT-RM samples were downloaded from ebi.ac.uk / ena / data / view / PRJEB19931.

[0326] PacBio sequencing

[0327] gDNA samples were purchased from the Coriell Institute for Medical Research (Coriell, NJ, USA). The quality of gDNA samples was evaluated by Nanodrop (ThermoFisher, MA, USA). The A280 / A260 ratio needed to be in the range of 1.8 - 2.0, and the A260 / 230 ratio ≥ 2.0. The molecular weight of gDNA was evaluated by a femtosecond pulse system (Agilent CA, USA). Most DNA fragment sizes should be > 40 kb. If the quality of gDNA samples from Coriell was below the protocol requirements, fresh DNA was extracted from B-lymphocytes (Coriell, NJ, USA) using a Qiagen DNA extraction kit (Qiagen, CA, USA).

[0328] According to the manufacturer's instructions (Covaris, MA, USA), 10 μg of gDNA was fragmented into 15 kb using a Covaris g-Tube. DNA was purified using 0.45× AMPure XP beads (Beckman Coulter, IN, USA) according to the manufacturer's instructions. The size of the sheared DNA was confirmed by a femtosecond pulse system (Agilent, CA, USA).

[0329] Construct the library according to the protocol of PacBio "Preparing HiFi Library Using SMRTbell Template Preparation Kit 1.0" or "Preparing HiFi Library Using SMRTbell Express Template Preparation Kit 2.0" (PacBio, CA, USA). Select a library size of 15 - 20 kb using a SageElf instrument (SageScience, MA, USA) with 0.75% agarose. Quality control of all libraries was performed using Qubit (Life Technologies, CA, USA) and femtosecond pulse (Agilent, CA, USA). library" or "Preparing HiFi library" of the protocol (PacBio, CA, USA). Use a SageElf instrument (SageScience, MA, USA) with 0.75% agarose to select a library size of 15 - 20 kb. Quality control of all libraries was performed using Qubit (Life Technologies, CA, USA) and femtosecond pulse (Agilent, CA, USA).

[0330] The PacBio Sequel II sequencing platform was used for sequencing. 20× coverage WGS data were generally obtained from 2-3 SMRT cells (Pacific Biosciences, CA, USA). CYP2D6 genotyping method

[0331] The method Cyrius described in this example first calls the sum of the copy numbers (CN) of CYP2D6 / 7 in a method similar to that described in Example 1. Using all reads mapped to CYP2D6 or CYP2D7 (including reads with a mapping quality of zero), the read count was directly calculated from the BAM file of the WGS alignment to account for regions with high sequence homology. The sum of the read counts was normalized by the region length. Then, GC correction was performed on 3000 preselected 2 kb regions across the genome. These 3000 normalized regions were randomly selected from the genome for stable coverage of cross-population samples to infer the sequencing depth and capture GC bias. The normalized depth values of the entire population were modeled using a one-dimensional mixture of 11 Gaussians, where the constrained means centered on each integer CN value represent the CN states in the range of 0 to 10. The CN of CYP2D6 + CYP2D7 was called from the Gaussian mixture model (GMM), where the posterior probability threshold was 0.95. The same method was used to call the CN of the 1.5 kb spacer region between repeat REP7 and CYP2D7 to infer the CN of the REP7-containing fusion gene ( Figure 23 ).

[0332] Figure 23A non-limiting exemplary graph showing the quality of WGS data in the CYP2D6 / 7 region. The average mapping quality of 1kGP samples is plotted for each position in the CYP2D6 / 7 region. A median filter is applied in a 200bp window. Nine exons of REP6, REP7, and CYP2D6 / 7 are boxed on the left (CYP2D6) and right (CYP2D7). The two 2.8kb repeat regions downstream of CYP2D6 (REP6) and CYP2D7 (REP7) are identical and largely non-alignable. The dashed box indicates the spacer region between CYP2D7 and REP7. The two major homology regions within the gene are shaded.

[0333] The method identified 118 CYP2D6 / CYP2D7 discriminatory bases (see additional information in this example, Figure 26 ). At each of these discriminatory base positions, Cyrius calls the number of chromosomes carrying CYP2D6 and the number of chromosomes carrying CYP2D7 by combining the total CYP2D6 + CYP2D7 CN with the read counts supporting each base in the gene-specific bases. Based on the total CN called, Cyrius iterates through all possible combinations of CYP2D6 CNs and CYP2D7 CNs and derives the combination that produces the highest posterior probability for the observed number of reads supporting CYP2D6 and CYP2D7. When the CN of CYP2D6 changes within the gene, gene fusions are called by identifying bases ( Figure 27 ).

[0334] Cyrius resolved read alignments to identify small variants defining star alleles. Variants of interest were divided into variants belonging to the CYP2D6 / CYP2D7 homology region (i.e., Figure 23 the low mapping quality region on ) and variants occurring in the unique region of CYP2D6. For the former, Cyrius looked for variant reads in CYP2D6 and their corresponding sites in CYP2D7. For the latter, Cyrius used reads aligned to CYP2D6. The CN called in this region was also considered during small variant calling. For example, in a sample where the *68 repeat fusion has been identified, one haplotype should have a full copy of CYP2D6 plus a copy of *68, while the other haplotype should have a full copy of CYP2D6, so the CYP2D6 CN should be at position 3 upstream of exon 2 and position 2 downstream of exon 2.

[0335] Finally, Cyrius matches the called structural and minor variants with the definition of star alleles (downloaded and parsed from PharmVar, pharmvar.org / gene / CYP2D6, last accessed in March 2019) to call star alleles, and further groups star alleles into haplotypes when, for example, there are more than two copies of CYP2D6. For this prior, information defining the exact haplotypes is included, e.g., *68 is on the same haplotype as *4, *36 is on the same haplotype as *10). These priors are made based on the tandem arrangement patterns described in PharmVar and are also supported by our real data (12 / 12 for *68 and 25 / 25 for *36). An option is available to match only the called structural and minor variants with star alleles of known function.

[0336] Of the 131 star alleles defined in PharmVar (last accessed in March 2020), 25 star alleles are still awaiting full screening, so this example excludes these alleles and focuses mainly on the 106 screened star alleles (another option is provided in Cyrius to include those unscreened star alleles). Among these 106 star alleles, four star alleles are removed from our target list, and none of these four star alleles are in GeT-RM. The removed star alleles include *61 and *63 (both of unknown function), which are CYP2D6 / 7 hybrid genes, very similar to *36, with the fusion breakpoint slightly upstream. Since it is impossible to distinguish the exon 7-exon 8 region between CYP2D6 / 7 ( Figure 26 ), these two star alleles cannot be distinguished from *36 and will be called *36 by Cyrius. Additionally, *27 (normal function) and *32 (unknown function) are removed; *27 and *32 share g.42126938C>T, which is a gene conversion variant in a highly homologous region (variant reads will align perfectly with CYP2D7). By counting the reads supporting CYP2D6 and CYP2D7 at a single locus, it may be difficult to accurately distinguish 1 copy of CYP2D6 and 3 copies of CYP2D7 from their respective 2 copies. Therefore, *27 will be called *1, and *32 will be called *41.

[0337] Validation according to the truth of GeT-RM and long-read values

[0338] When comparing CYP2D6 calls made by Cyrius, Aldy, and Stargazer with the consensus genotypes provided by GeT-RM, genotypes are considered to match as long as all star alleles in the star allele of the true genotype are present, even if the haplotype assignments are different. Examples of this situation occur in several samples where GeT-RM lists *1 / *10+*36+*36 but is called by Aldy as *1+*36 / *10+*36.

[0339] When validating genotype calls for PacBio data, PacBio reads covering the entire CYP2D6 gene are analyzed to identify small variants that define known star alleles. The long (∼10 kb) reads allow these variants to be fully phased into haplotypes, and these haplotypes are matched against the star allele table to determine which star allele each read represents. Reads carrying structural variants are identified by aligning the reads to a set of reference alleles that are constructed to represent known structural variants (*5 / *13 / *36 / *68 / repeat fusions).

[0340] Running Aldy and Stargazer

[0341] Run Aldy v2.2.5 using the command "aldy genotype-p illumina-g CYP2D6".

[0342] Use VDR as a control gene and run Stargazer v1.0.7 with GDF and VCF files as input to genotype CYP2D6.

[0343] Since Aldy and Stargazer only support GRCh37, for 1kGP samples that were initially aligned to hs38DH, use Isaac to realign to GRCh37.

[0344] Results

[0345] Validation and performance comparison

[0346] Cyrius, Aldy, and Stargazer made CYP2D6 calls for 188 samples for which high-quality ground truth was obtained. These 188 samples included 138 GeT-RM samples, and 50 fact samples with PacBio whole-genome sequencing were compared (Tables 20, 21). PacBio CCS data allowed the breakpoints of common and rare structural variants in this region to be located and visualized ( Figure 24), and thus serve as a valuable resource for studying complex star alleles and confirming the phasing of star allele variants. In the case of short reads, these samples with SVs show different depth signals that allow accurate calling of SVs ( Figure 27 )

[0347] Table 20. Summary of benchmark test results against the truth .

[0348]

[0349] * After seeing three inconsistent samples, Cyrius was improved, and then Cyrius was able to accurately call 187 out of 188 of these samples.

[0350] Table 21. Results of Cyrius / Aldy / Stargazer for GeT-RM and PacBio truth .

[0351]

[0352]

[0353]

[0354]

[0355]

[0356]

[0357]

[0358]

[0359] By comparing with GeT-RM samples, three samples were found in which the calls of all three callers were either consistent or inconsistent with the GeT-RM consensus sequence. Whole-genome PacBio sequencing confirmed that the calls of the three callers were correct and that the GeT-RM consensus sequence should be updated ( Figure 24 )

[0360] Figure 24Shows structural variants validated by PacBio CCS reads. PacBio reads support deletions (*5), duplications, and fusions (*36, *68, and *13). Curves were generated using sv-viz2 (zotero.org / google-docs / ?xAunA6). For deletions and duplications, due to the identical sequences in the REP6 / 7 region, the exact positions of breakpoints within REP6 / 7 are not available. The breakpoints in A and B are for illustrative purposes only. The genotypes of the samples in subfigures A to E are *2 / *5, *17 / *2x2, *10 / *10+*36, *29 / *4+*68, and *1 / *2+*13, respectively.

[0361] Cyrius initially made four discordant calls from the truth set GeT-RM, showing a sensitivity of 97.9%. Among these discrepancies was sample NA19908 (defined as *1 / *46 by GeT-RM), where Cyrius called *1 / *46 and *43 / *45 as two possible haplotypes. Both of these star allele combinations produce the same set of variants. Neither phasing of reads nor population frequency analysis could rule out either genotype combination. Genotyping results from various assays of the GeT-RM consensus sequence from which this sample was generated also showed discrepancies between *1 / *46 and *43 / *45, highlighting the difficulty of these combinations (Table 22). Future sequencing of more samples for either haplotype may help identify new variants that distinguish between the two.

[0362] Table 22. GeT-RM results for sample NA19908 。

[0363]

[0364] In the remaining three samples where Cyrius was discordant with the truth set, the errors were determined and Cyrius was improved to call the correct genotype. First, in NA23275 (*1 / *40), an 18bp insertion defining *40 was initially missed because reads containing the insertion were typically not aligned as insertions but as soft clips. When looking for variants, the caller was improved to consider soft clips. Second, in HG03225 (*5 / *56), CYP2D7-derived reads aligned to CYP2D6, preventing the variant defined by *56 from being called. The caller was improved to be more sensitive to variant reads in this region. Finally, in HG00421 (*10x2 / *2), a fusion was mis-called as *36, as was the case with two other callers. A more detailed examination of this sample with PacBio data showed a different fusion, *10D, where the fusion breakpoint was located downstream of exon 9 ( Figure 28)。This fusion has the same function (reduced function) as *10, while *36 is non-functional due to CYP2D7-derived exon 9. The caller was improved to be able to call *10D. Although these three samples were considered mis-called in this example, the improvements made to Cyrius after seeing these three samples allowed accurate calling of 187 out of 188 samples, highlighting how more real data and more population data can identify limitations of callers that can be improved for subsequent samples.

[0365] In contrast, when compared to these samples, the sensitivity of two other CYP2D6 callers was less than 90%. The sensitivity of Aldy was 88.8%. Specifically, Aldy over-called CYP2D6 / CYP2D7 fusions such as *61, *63, *78, and *83 (calling 8 out of 21 inconsistent samples, Table 21). Figure 29 The PacBio data in can demonstrate that the fusions called by Aldy are incorrect. Stargazer had a sensitivity of 85.6% and was most error-prone when SVs were present. The sensitivity in samples with SVs was only 77.8%, and 16 out of 27 inconsistent calls were from samples with structural variants. Notably, Stargazer mis-called NA19317 (*5 / *5) as *2 / *2, completely missing the double deletion. Stargazer could not genotype two samples with the *13 fusion (Table 21). In addition, Stargazer showed a high error rate in the *36 fusion (7 mis-calls out of a total of 25 samples with *36). Specifically, Stargazer mis-called all 5 samples where there was more than one copy of *36 in a single haplotype.

[0366] The 188 validation samples used in this example together confirmed the accuracy of CYP2D6 calls in 48 different haplotypes (Table 23), including 41 star alleles and several common and rare SV structures such as duplications, *2+*13, *4+*68, *10+*36, *10+*36+*36, and *10+*36+*36+*83 (novel haplotypes not previously reported, see Figure 30A and Figure 30B ). These 41 star alleles tested in the validation data represent 38.7% of the 106 screened star alleles currently listed in PharmVar and 53.4% of those alleles with known function (31 out of 58). They overlap 96.4% with the Cyrius haplotypes called from 1kGP samples (Table 23, also see the next section).

[0367] Table 23. Haplotypes validated in this example and their frequencies in the 1kGP .

[0368]

[0369]

[0370]

[0371]

[0372] CYP2D6 haplotype frequencies in five ethnic groups

[0373] Given the high accuracy demonstrated in the previous section, in addition to the validation samples, Cyrius was used to study CYP2D6 in the global population. The haplotype distributions of populations (Europeans, Asians, East Asians, South Asians, and admixed Americans consisting of Caucasians, Mexican Americans, Peruvians, and Puerto Ricans) in 2,504 1kGP samples were analyzed ( Figure 25 , Table 23). Cyrius made an unambiguous diplotype call in 2,445 (97.6%) of the 2,504 samples that called 46 different star alleles, and 41 of the star alleles overlapped with those already included in the validation data. These 41 previously validated star allele calls accounted for 96.5% of all star allele calls in the 1kGP samples (Table 23).

[0374] Figure 25 As a non-limiting exemplary graph, it shows the frequencies of CYP2D6 alleles of the ten most common haplotypes with altered CYP2D6 function in five ethnic groups. One haplotype (*2x2) has enhanced function, two haplotypes (*4 and *4+*68) have no function, and the remaining haplotypes have reduced function.

[0375] In 59 samples where Cyrius did not make an unambiguous diplotype call, 10 samples had ambiguous SV calls, 30 samples had variant calls that did not match any of the known star alleles, four samples had the same ambiguity between *1 / *46 and *43 / *45 as described in the validation sample NA19908 above, and 15 samples had unambiguous star allele calls that Cyrius could not unambiguously phase into a diplotype.

[0376] In most cases, the haplotype frequencies were consistent with pharmGKB ( Figure 31A and Figure 31B, Table 24). For example, Africans have high frequencies of *17 (∼20%) and *29 (∼10%), South Asians have a high frequency of *41 (∼12%), Europeans have a high frequency of *4 (18% to 20%, including *4 + *68), and East Asians have a high frequency of *10 (40% to 50%, including *10 + *36). As Cyrius improves in sensitivity to structural variants, it is possible to provide a more comprehensive picture of the frequencies of structural variants across populations. Among them, the haplotype *10 + *36 containing the fusion is very common in East Asians (>30%, reported in PharmGKB as 1% to 2%, Figure 31A and Figure 31B ), and another haplotype *4 + *68 containing the fusion is also quite common in Europeans (>5%, data not available in PharmGKB, Figure 31A and Figure 31B ). In summary, in East Asians, Europeans, Americans, Africans, and South Asians, the frequencies of haplotypes involving SVs are estimated to be 32.2%, 5.57%, 1.47%, 1.34%, and 0.45% higher, respectively, than those reported in PharmGKB (the total frequencies reported in PharmGKB are 7.48%, 5.33%, 5.17%, 9.9%, and 6.19%, respectively).

[0377] There are several other haplotypes reported to have lower frequencies than PharmGKB ( Figure 31A and Figure 31B), thus highlighting the difficulty of combining data from multiple studies using different techniques. These include *2 in Africans and South Asians. Since *2 is the default assignment, its frequency can be overestimated in PharmGKB if some other star alleles are not tested. A lower frequency of *41 was determined in Africans. According to PharmGKB, *41 is not always determined by the SNPs that define it in various studies, leading to an overestimation of the frequency of *41, especially for Africans. The much higher frequency of *29 in South Asians in PharmGKB (estimated at 6% vs. 0% in this example) is caused by an error in PharmGKB: 0.2% was incorrectly included in PharmGKB as 20% in the article by Sistonen et al. (CYP2D6 worldwide genetic variation shows high frequency of altered activity variants and no continental structure. Pharmacogenet Genomics, 2007, Vol. 17, No. 2: pp. 93–101, doi: 10.1097 / 01.fpc.0000239974.69464.f2). In Europeans, the frequencies of *34 and *39 were estimated to be much lower. *34 and *39 are each defined by one of two variants that define *2, so both of these variants should have been tested in any study mapping CYP2D6.*34 and *39 were reported in only 3 out of 91 studies on Europeans at >1% in PharmGKB. Among them, Wesmiller et al. (The Association of CYP2D6 Genotype and Postoperative Nausea and Vomiting in Orthopedic Trauma Patients. Biol Res Nurs., 2013, Vol. 15, No. 4, pp. 382 - 389, doi:10.1177 / 1099800412449181) reported only *39 with a limited sample size (N = 112), Kapedanovska Nestorovska (Distribution of the most Common Genetic Variants Associated with a Variable Drug Response in the Population of the Republic of Macedonia, Balk J Med Genet BJMG, 2014, Vol. 17, No. 2, pp. 5 - 14, doi:10.2478 / bjmg - 2014 - 0069) reported both *34 and *39, and was for a specific country Macedonia and also had a small sample size (N = 184), and Del Tredici et al. (Frequency of CYP2D6 Alleles Including Structural Variants in the United States, Front Pharmacol, 2018, Vol. 9, doi:10.3389 / fphar.2018.00305) did not report *34 or *39, but PharmGKB may have wrongly reported the frequency of *35 as the frequency of *34.

[0378] Analysis of CYP2D6 / CYP2D7 differentiating bases

[0379] A total of 208 single - nucleotide differences between CYP2D6 / 7 were extracted from the reference genome. Among the 1kGP samples with a total CN of CYP2D6 + CYP2D7 of 4, i.e., un - called structural variants, the percentage of samples in which the CN of the CYP2D6 base was called as 2 at 208 loci was queried ( Figure 26)。Many loci showed a low percentage of samples with two copies of the CYP2D6 base, indicating that the CYP2D6 / CYP2D7 base difference is not fixed in the population and thus the base difference cannot be used to distinguish the two genes. Read alignment relying on these loci will generate a large amount of noise when distinguishing the two genes. A total of 118 highly stable loci were selected, where >98% of the samples showed two copies of the CYP2D6 base for CYP2D6 / CYP2D7 distinction, which allowed obtaining the cleanest signal for calling SV.

[0380] Additional figures and tables

[0381] Figure 26 The CYP2D6 / CYP2D7 base difference loci are shown to have high variability in the population. The Y-axis shows the sample frequency in which the CN of the CYP2D6 base is called 2 among all samples with a total CYP2D6 + CYP2D7 CN of 4. The X-axis shows the genomic coordinates in hg38. The CYP2D6 exons are plotted as gray boxes above the figure. The black horizontal line represents the 98% cutoff.

[0382] Figure 27 The raw CYP2D6 CN across the CYP2D6 / 7 distinction loci in an example with SV is shown. The raw CYP2D6 CN is calculated as the total CYP2D6 + CYP2D7 CN multiplied by the ratio of CYP2D6 supporting reads to the total number of CYP2D6 and CYP2D7 supporting reads. The large diamonds represent the copy number of the CYP2D6-derived gene (which can be a full CYP2D6 or a fusion gene ending with CYP2D6) at the end of the gene, calculated as the total CN of CYP2D6 + CYP2D7 minus the CN of the CYP2D7 spacer region (see Figure 23 )。To detect SV, the CYP2D6 CN is called at each locus, and a change in the CYP2D6 CN within the gene indicates the presence of SV. For example, in HG01161, the CYP2D6 CN changed from 2 to 1 between exon 7 and exon 9, indicating a CYP2D7-CYP2D6 hybrid gene. In HG00553, the CYP2D6 CN changed from 2 to 3 between exon 1 and exon 2, indicating a CYP2D6-CYP2D7 hybrid gene.

[0383] Figure 28PacBio data shows confirmation of the *10D fusion in HG00421. In the comparison, a sample with *36 (HG00612) is shown. PacBio reads containing the fusion are those with shaded bases, which represent soft clips prepared by the aligner and are derived from the CYP2D7 portion of the fusion. The fusion breakpoints are close to each other, but the breakpoint of *36 is upstream of the base difference in exon 9 (those within the black box), while the breakpoint of *10D is downstream, thus keeping the CYP2D6 gene intact.

[0384] Figure 29 PacBio data shows a false *61 (CYP2D6 / CYP2D7 hybrid) call made by Aldy in HG02622. The expected genotype was *17 / *45, but Aldy called *61-like / *78 (both *61 and *78 are star alleles with SVs). PacBio data shows no structural variants in this region (each read aligns completely, and no soft clips indicate unaligned portions).

[0385] Figure 30A and Figure 30B shows a novel *10+*36+*36+*83 haplotype in HG00597. Figure 30A . The depth plot, as shown in Figure 27 , shows that HG00597 has three copies of the *36-like fusion, all of which have breakpoints in the homologous region between exon 7 and exon 9. Figure 30B . An IGV screenshot of the PacBio data, which shows all reads containing the fusion, i.e., those that align to soft clips. One copy of the fusion gene does not have g.42130692G>A, an SNP that is in *36 but not in *83, as shown in the region flanked by two black vertical lines. This copy is *83 and is different from the copy reported in PharmVar, which is a fusion gene with REP7 instead of REP6, otherwise the copy number in the region downstream of exon 9 would be 3 instead of Figure 30A 2 in

[0386] Figure 31A and Figure 31BCompare between 1kGP and pharmGKB frequencies. Each point represents a haplotype with frequency >= 0.5% in 1kGP or pharmGKB. SV-related haplotypes are marked, including the two haplotypes with the largest deviations (*10+*36 in East Asians and *4+*68 in Europeans). Other haplotypes with deviation values are annotated (*2, *41, *34, *39, *2, and *29). Draw a diagonal line for each subplot. The correlation coefficients for each population are listed (*10+*36 is excluded in East Asians and *4+*68 is excluded in Europeans for calculation). Figure 31B Values in the low value range (< 5%) are shown.

[0387] Figure 32 For non-limiting exemplary IGV snapshots, which show the de novo assembly of PacBio reads in HG00733 excluding the *68 fusion.

[0388] Table 24. Comparison of haplotype frequencies called by Cyrius and pharmGKB frequencies

[0389]

[0390]

[0391]

[0392]

[0393] Discussion

[0394] These embodiments describe Cyrius, a method capable of accurately diploidizing difficult CYP2D6 regions. The unique feature of this embodiment is the use of long-read data to validate both haplotypes and SVs. Long reads provide a unique opportunity to confirm the breakpoint regions of common SVs (CYP2D6 deletions and duplications, and CYP2D6 / 7 fusion genes) and to phase the CYP2D6 gene. Using 188 samples (including 50 samples with long-read validation data) as an orthogonal validation dataset, it is shown that Cyrius outperforms other CYP2D6 genotypers, achieving 97.9% accuracy, while Aldy achieved 88.8% and Stargazer achieved 85.6%. Specifically, compared to these existing CYP2D6 callers, Cyrius allows for the possibility that reads may not align in regions of high similarity between CYP2D6 / 7. Ambiguous read alignments in these regions can lead to incorrect copy number estimates and errors in small variant calling. By considering potentially misaligned reads and selecting a set of reliable CYP2D6 / 7 differentiating sites, Cyrius is able to better identify star alleles with SVs, resulting in 97.2% accuracy, compared to Aldy's 88.9% accuracy and Stargazer's 77.8% accuracy.

[0395] Among these 188 validation samples, a total of 41 different star alleles were validated, representing 38.7% of all star alleles listed in PharmGKB, including 53.4% of star alleles with known functional status. Even though the validation set included only 38.7% of the total known star alleles, based on the analysis of 1kGP samples in this embodiment, it is estimated that these alleles account for 96.5% of star alleles in the genome-wide population. Generally speaking, the allele frequencies calculated for 2504 1kGP samples from five ethnic groups are consistent with previous studies on simple star alleles. In contrast, for some star alleles defined by the presence of SVs, quite different frequencies were identified, which may be because many star alleles affected by SVs are difficult to resolve with conventional assays. This highlights the inherent errors in the combined results from studies using multiple different CYP2D6 assays, some of which may be designed to call only a subset of star alleles. For example, among the 5 assays used to generate the GeT-RM consensus genotype, the individual accuracies ranged from 47.1% to 75.2% compared to the consensus sequence (Table 25). A single method capable of resolving all known star alleles from a single assay is a better choice for constructing a population-level database.

[0396] Table 25. Accuracy of individual GeT-RM assays 。

[0397]

[0398] In addition, Cyrius was used to analyze 2,504 1kGP samples from five ethnic groups to determine star allele frequencies. The calculated allele frequencies were consistent with previous studies of simple star alleles, and Cyrius greatly improved the estimation of allele frequencies for star alleles involving structural variants that may be difficult to detect by conventional methods.

[0399] Some existing methods rely on accurate read alignment to distinguish CYP2D6 and CYP2D7, which can be error-prone due to several highly sequence-similar regions between the two genes, especially the intron 1-exon 2 and exon 7-exon 9 regions. Fuzzy alignment can lead to noise in the depth map, resulting in incorrect CNV calls. Additionally, incorrect read alignment can lead to false positive or negative variant calls. In contrast, Cyrius first calls the total CN of CYP2D6+CYP2D7 by counting all reads aligned to either gene, and a total CN not equal to 4 clearly indicates the presence of an SV. To determine the exact location of the SV, not all differences based on the reference genome are used. Many CYP2D6 / CYP2D7 base differences are not fixed, so not all of these positions can be used to reliably distinguish CYP2D6 from CYP2D7 ( Figure 26 ). Cyrius uses the selected 118 CYP2D6 / CYP2D7 discriminatory positions to determine the exact location of the SV. By first calling the total CN and then differentiating them using a subset of well-discriminating bases, Cyrius is able to achieve more accurate SV calls. For small variant calls, Cyrius overcomes the dependence on unambiguous alignment by looking for variant reads at both CYP2D6 positions and corresponding positions in CYP2D7, resulting in the most accurate small variant calls.

[0400] In this example, long-read data was used to validate both haplotypes and SV calls. The PacBio data in this example provided a clear picture of the CYP2D6-CYP2D7 region with high-quality long reads (10 kb to 20 kb). Specifically, the PacBio data helped to resolve the breakpoint regions of common structural variants (CYP2D6 deletions and duplications, and CYP2D6-CYP2D7 fusion genes). Even for PacBio reads, genotyping CYP2D6 may not be straightforward and may require targeted analysis, especially for structural variants involving duplications (CYP2D6 duplications and CYP2D6-CYP2D7 duplicated fusions), where the duplicated region >10 kb. For example, de novo assembly methods could not capture the *68 fusion in sample HG00733 ( Figure 31A and Figure 31B)。In addition, the length of PacBio reads is not sufficient to cover more than one copy of the repetitive sequence, and PacBio reads are too long for read count-based CN calling (for short reads), making it difficult to estimate the copy number. Whole-genome sequencing with short reads provides the most accurate solution for genotyping CYP2D6.

[0401] In the analysis of 1kGP samples, Cyrius was able to call a definitive genotype in greater than 97.6% of the samples. In some embodiments, Cyrius can resolve the remaining 2.4% of the samples. For example, in samples where multiple haplotype configurations are possible, it may be useful to employ a probabilistic method to derive the most likely genotype given the observed variants. Additionally, continued sequencing and testing of more samples will help confirm the ability to genotype rare star alleles and will also identify new variants that can be used to distinguish ambiguous diplotypes. This process is demonstrated in this example, where improvements were made to better call three star alleles that were initially mis-called in 188 validation samples. These improvements are beneficial for population-level genotyping because these three star alleles are present in nearly 1% (23 out of 2504) of the 1kGP samples.

[0402] As new star alleles are identified, these new star alleles can be added to the Cyrius database. One consideration in adding new star alleles defined by new variants is that these variants are unlikely to have been considered in previous star allele definitions. Thus, there may be novel combinations of new and existing variants that do not match any known combinations, resulting in no calls. For example, Cyrius includes the option to genotype 25 new star alleles added in PharmVar v4 (not included in GeT-RM, Aldy, or Stargazer). However, five of the 25 new star alleles (*119, *122, *135, *136, *139) have new variants that, when included, result in no calls for samples that were previously callable, indicating the existence of common novel star alleles with variant combinations not captured in PharmVar. Thus, these five star alleles were removed along with two other alleles (*127, which has a gene conversion variant in the homologous region, and *131, which has a variant in a noisy site), leaving the remaining 18 alleles. When new variants / star alleles are identified, novel star alleles are possible. Public WGS datasets such as the 2504 1kGP samples analyzed herein can be an important part of integrating new variants into star allele definitions because this data will allow for rapid assessment of variants in many samples with different genotypes.

[0403] Whole genome sequencing (WGS) provides a valuable opportunity to map all genetic variations across the entire genome, but many clinically important regions / variants are beyond the capabilities of most secondary analysis pipelines. CYP2D6 is one of the difficult regions in the genome that is both clinically important and requires targeted bioinformatics solutions beyond normal WGS pipelines. Such targeted approaches have been successfully applied to some difficult regions, such as the SMN1 gene responsible for spinal muscular atrophy, as shown in Example 1. More targeted methods such as Cyrius can accelerate pharmacogenetics, enabling personalized medicine.

[0404] Additional considerations

[0405] In at least some of the foregoing embodiments, one or more elements used in one embodiment may be interchangeably used in another embodiment, unless such substitution is technically infeasible. Those skilled in the art will understand that various other omissions, additions, and modifications may be made to the above methods and structures without departing from the scope of the claimed subject matter. All such modifications and changes are intended to fall within the scope of the subject matter defined by the appended claims.

[0406] Those skilled in the art will understand that for such processes and methods as well as other processes and methods disclosed herein, the functions performed in these processes and methods may be implemented in a different order. In addition, the steps and operations outlined are provided only as examples, and some of these steps and operations may be optional, combined into fewer steps and operations, or extended into additional steps and operations without detracting from the essence of the disclosed embodiments.

[0407] Regarding the use of substantially any plural and / or singular terms herein, those skilled in the art may appropriately convert from plural to singular and / or from singular to plural depending on the context and / or application. For clarity, various singular / plural permutations may be explicitly shown herein. As used in this specification and 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 is configured to" are intended to include one or more of the said devices. Such one or more of the said devices may also be jointly configured to perform the stated expression. For example, "a processor configured to perform expressions A, B, and C" may include a first processor configured to perform expression A and work in cooperation with a second processor configured to perform expressions B and C. Unless otherwise indicated, any reference to "or" herein is intended to include "and / or".

[0408] Those skilled in the art should understand that, generally speaking, the terms used herein, especially the terms in the appended claims (e.g., the subject matter of the appended claims), are generally intended to be "open" terms (e.g., the term "comprising" should be interpreted as "comprising but not limited to", the term "having" should be interpreted as "having at least", the term "including" should be interpreted as "including but not limited to", etc.). Those skilled in the art should also understand that if the specific quantity of the introduced claim expression is intended, such intention will be clearly expressed in the claim, and in the absence of such expression, there is no such intention. For example, for the sake of understanding, the following appended claims may include the use of the introductory phrases "at least one" and "one or more" to introduce claim expressions. However, even when the same claim includes the introductory phrases "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 as meaning "at least one" or "one or more"), the use of such phrases should not be understood as implying that introducing a claim expression by the indefinite article "a" or "an" will limit any particular claim containing such introduced claim expression to an embodiment containing only one such expression; the same applies to the use of a definite article to introduce a claim expression. In addition, even if the specific quantity of the introduced claim expression is clearly stated, those skilled in the art will also recognize that such expression should be interpreted as meaning at least the stated quantity (e.g., in the absence of other modifiers, a direct statement of "two expressions" means at least two expressions, or two or more expressions). In addition, in those cases where a convention similar to "at least one of A, B, and C, etc." is used, generally speaking, such convention is intended to be used in the sense that those skilled in the art will understand the meaning of the convention (e.g., "a system having at least one of A, B, and C" will include but not be limited to a system having A alone, B alone, C alone, A and B together, A and C together, B and C together, and / or A, B, and C together, etc.). In those cases where a convention similar to "at least one of A, B, or C, etc." is used, generally speaking, such convention is intended to be used in the sense that those skilled in the art will understand the meaning of the convention (e.g., "a system having at least one of A, B, or C" will include but not be limited to a system having A alone, B alone, C alone, A and B together, A and C together, B and C together, and / or A, B, and C together, etc.). Those skilled in the art should also understand that, in fact, no matter in the specification, claims or drawings, any disjunctive words and / or phrases presenting two or more alternative terms should be understood as considering the possibility of including one of the terms, any one of the terms, or both of these terms. For example, the phrase "A or B" will be understood as including the possibility of "A" or "B" or "A and B".

[0409] In addition, where the features or aspects of the present disclosure are described in terms of Markush groups, those skilled in the art will recognize that the disclosure is also thereby described in terms of any single member or subgroup of members of the Markush group.

[0410] As will be understood by those skilled in the art, for any and all purposes, such as in providing a written description, all ranges disclosed herein also include any and all possible subranges and combinations of subranges thereof. Any listed range can be readily identified as being sufficiently described and enabling the same range to be broken down into at least equal halves, thirds, quarters, fifths, tenths, etc. By way of non-limiting example, each range discussed herein can be readily broken down into a lower third, middle third, and upper third, etc. As will also be understood by those skilled in the art, all language such as "at most", "at least", "greater than", "less than", etc. includes the recited numbers and refers to ranges that can then be broken down into the subranges as described above. Finally, as will be understood by those skilled in the art, a range includes each individual member. Thus, for example, a group having 1 - 3 elements refers to a group having 1, 2, or 3 elements. Similarly, a group having 1 - 5 elements refers to a group having 1, 2, 3, 4, or 5 elements, etc.

[0411] It should be understood that, for purposes of illustration, various embodiments of the present disclosure have been described herein, and various modifications can 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, where the true scope and spirit are indicated by the following claims.

[0412] It should be understood that not 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 certain embodiments can be configured to operate in a manner that achieves or optimizes one advantage or a group of advantages as presented herein without necessarily achieving other objectives or advantages as may be presented or suggested herein.

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

[0414] Many other variations besides those described herein will be apparent from this disclosure. For example, according to an embodiment, certain actions, events, or functions of any of the algorithms described herein may be performed in a different order, may be added, combined, or entirely omitted (e.g., not all of the described actions or events are necessary for the practice of the algorithm). Additionally, in certain embodiments, actions or events may be performed concurrently rather than sequentially, for example, by multithreading, interrupt processing, or multiple processors or processor cores, or on other parallel architectures. Further, different tasks or processes may be performed by different machines and / or computing systems that may run together.

[0415] The various illustrative logical blocks and modules described in connection with the embodiments disclosed herein may be implemented or performed by a machine designed to perform the functions described herein, such as a processing unit or processor, a digital signal processor (DSP), an application specific integrated circuit (ASIC), a field programmable gate array (FPGA) or other programmable logic device, discrete gate or transistor logic, discrete hardware components, or any combination thereof. A processor may be a microprocessor, but in the alternative, the processor may be a controller, a microcontroller, or a state machine, combinations thereof, etc. The processor may include 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 conjunction with a DSP core, or any other such configuration. Although described primarily in connection with digital technology 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 in hybrid analog and digital circuitry. By way of example, a computing environment may include any type of computer system, including but not limited to a microprocessor-based computer system, mainframe computer, digital signal processor, portable computing device, device controller, or computing engine within a device.

[0416] Any process descriptions, elements, or boxes in the flowcharts described and / or shown herein should be understood as potentially representing code modules, segments, or portions that include one or more executable instructions for implementing specific logical functions or elements in that process. As will be understood by those skilled in the art, alternative specific implementations are included within the scope of the embodiments described herein, where elements or functions may be deleted, performed in the order shown or discussed (including substantially concurrently or in reverse order), depending on the functions involved.

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

Claims

1. A system for determining the copy number of the survival motor neuron 1 (SMN1) gene, comprising: A non-transitory memory configured to store executable instructions and sequence data, the sequence data including a plurality of sequence reads obtained from a sample of a subject and aligned with the survival motor neuron 1 (SMN1) gene or the survival motor neuron 2 (SMN2) gene; and A hardware processor in communication with the non-transitory memory, the hardware processor programmed by the executable instructions to perform: Receiving sequence data including a plurality of sequence reads obtained from a sample of a subject and aligned with the survival motor neuron 1 (SMN1) gene or the survival motor neuron 2 (SMN2) gene; Determining (i) a first number of sequence reads of the plurality of sequence reads aligned with a first SMN1 or SMN2 region respectively including at least one of exons 1 to 6 of the SMN1 gene or the SMN2 gene and (ii) a second number of sequence reads of the plurality of sequence reads aligned with a second SMN1 or SMN2 region respectively including at least one of exons 7 and 8 of the SMN1 gene or the SMN2 gene; Using (i) the length of the first SMN1 or SMN2 region and (ii) the length of the second SMN1 or SMN2 region respectively to determine (i) a first normalized number of sequence reads aligned with the first SMN1 or SMN2 region and (ii) a second normalized number of sequence reads aligned with the second SMN1 or SMN2 region; Using a Gaussian mixture model including a plurality of Gaussian functions each representing a different integer copy number to determine (i) the copy number of the total survival motor neuron SMN genes each being a complete SMN1 gene, a complete SMN2 gene, a truncated SMN1 gene or a truncated SMN2 gene and (ii) the copy number of any complete SMN genes each being the complete SMN1 gene or the complete SMN2 gene, respectively considering (i) the first normalized number of sequence reads aligned with the first SMN1 or SMN2 region and (ii) the second normalized number of sequence reads aligned with the second SMN1 or SMN2 region; For a base among a plurality of SMN1 gene-specific bases associated with the complete SMN1 gene, determining the most likely combination among a plurality of possible combinations of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene each including a total number of copies of the determined any complete SMN genes, considering (a) the number of sequence reads of the plurality of sequence reads having a base supporting the SMN1 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads having a base supporting the SMN2 gene-specific base corresponding to the SMN1 gene-specific base; And Determine the copy number of the SMN1 gene using the most likely combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene determined for the SMN1 gene-specific bases.

2. The system according to claim 1, wherein the sequence data comprises whole genome sequencing WGS data or short read WGS data.

3. The system according to any one of claims 1 to 2, wherein the subject is a neonatal subject, a pediatric subject, an adolescent subject or an adult subject.

4. The system according to any one of claims 1 to 2, wherein the sample comprises cells or cell-free DNA.

5. The system according to any one of claims 1 to 2, wherein the sequence reads of the plurality of sequence reads are aligned with the first SMN1 or SMN2 region or the second SMN1 or SMN2 region, wherein the alignment quality score is approximately zero.

6. The system according to any one of claims 1 to 2, wherein the first SMN1 or SMN2 region respectively comprises exons 1 to 6 of the SMN1 gene or the SMN2 gene and has a length of approximately 22.2 kb, and wherein the second SMN1 or SMN2 region respectively comprises exons 7 and 8 of the SMN1 gene or the SMN2 gene and has a length of approximately 6 kb.

7. The system according to any one of claims 1 to 2, wherein determining (i) a first normalized quantity of the sequence reads aligned to the first SMN1 or SMN2 region and (ii) a second normalized quantity of the sequence reads aligned to the second SMN1 or SMN2 region comprises: Use (i) the length of the first SMN1 or SMN2 region and (ii) the length of the second SMN1 or SMN2 region respectively to determine (i) the first normalized number of the sequence reads aligned with the first SMN1 or SMN2 region and (ii) the second normalized number of the sequence reads aligned with the second SMN1 or SMN2 region, and to determine (iii) the depth of the sequence reads of the region of the genome of the subject other than the locus containing the SMN1 gene and the SMN2 gene in the sequence data.

8. The system according to claim 7, wherein determining (i) the first normalized number of the sequence reads aligned with the first SMN1 or SMN2 region and (ii) the second normalized number of the sequence reads aligned with the second SMN1 or SMN2 region comprises: Use (i) the length of the first SMN1 or SMN2 region and (ii) the length of the second SMN1 or SMN2 region respectively to determine (i) the first SMN1 or SMN2 region length-normalized number of the sequence reads aligned with the first SMN1 or SMN2 region and (ii) the second SMN1 or SMN2 region length-normalized number of the sequence reads aligned with the second SMN1 or SMN2 region; And Using the depth of sequence reads of regions of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene, the first normalized depth of the sequence reads aligned to the first SMN1 or SMN2 region and the second normalized depth of the sequence reads aligned to the second SMN1 or SMN2 region are determined respectively according to (i) the first SMN1 or SMN2 region length-normalized quantity and (ii) the second SMN1 or SMN2 region length-normalized quantity, and the first normalized quantity of the sequence reads aligned to the first SMN1 or SMN2 region and the second normalized quantity of the sequence reads aligned to the second SMN1 or SMN2 region are the first normalized depth and the second normalized depth respectively.

9. The system according to any one of claims 1 to 2, wherein determining (i) a first normalized quantity of the sequence reads aligned to the first SMN1 or SMN2 region and (ii) a second normalized quantity of the sequence reads aligned to the second SMN1 or SMN2 region comprises: Using (i) the GC content of the first SMN1 or SMN2 region and (ii) the GC content of the second SMN1 or SMN2 region respectively to determine (i) the first normalized quantity of the sequence reads aligned to the first SMN1 or SMN2 region and (ii) the second normalized quantity of the sequence reads aligned to the second SMN1 or SMN2 region, and to determine (iii) the depth of sequence reads of regions of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene in the sequence data, and to determine (iv) the GC content of the regions of the genome.

10. The system according to claim 7, wherein the depth of the region comprises the average depth or the median depth of sequence reads of regions of the subject's genome other than the locus containing the SMN1 gene and the SMN2 gene in the sequence data.

11. The system according to claim 10, wherein the region comprises about 3000 preselected regions each having a length of about 2 kb and each spanning the subject's genome.

12. The system according to any one of claims 1 to 2, wherein (i) the first normalized quantity of the sequence reads aligned to the first SMN1 or SMN2 region and / or (ii) the second normalized quantity of the sequence reads aligned to the second SMN1 or SMN2 region is about 30 to about 40.

13. The system according to any one of claims 1 to 2, wherein the Gaussian mixture model comprises a one-dimensional Gaussian mixture model.

14. The system according to any one of claims 1 to 2, wherein the plurality of Gaussian functions of the Gaussian mixture model represent integer copy numbers from 0 to 10.

15. The system according to any one of claims 1 to 2, wherein the mean of each Gaussian function in the plurality of Gaussian functions is the integer copy number represented by the Gaussian function.

16. The system according to any one of claims 1 to 2, wherein determining (i) the copy number of the total SMN gene and (ii) the copy number of any full-length SMN gene comprises, respectively, taking into account (i) the first normalized quantity of the sequence reads aligned to the first SMN1 or SMN2 region and (ii) the second normalized quantity of the sequence reads aligned to the second SMN1 or SMN2 region, and using the Gaussian mixture model and a first predetermined posterior probability threshold to determine (i) the copy number of the total SMN gene and (ii) the copy number of any full-length SMN gene.

17. The system according to claim 16, wherein the first predetermined posterior probability threshold is 0.

95.

18. The system according to any one of claims 1 to 2, wherein the hardware processor is programmed by the executable instructions to perform: determining the copy number of the truncated SMN gene using (i) the determined copy number of the total SMN gene and (ii) the determined copy number of the full-length SMN gene.

19. The system according to claim 18, wherein the copy number of the truncated SMN gene is the difference between (i) the determined copy number of the total SMN gene and (ii) the determined copy number of the full-length SMN gene.

20. The system according to any one of claims 1 to 2, wherein the SMN1 gene-specific base is a splicing enhancer.

21. The system according to any one of claims 1 to 2, wherein the SMN1 gene-specific base is the base at c.840 of the SMN1 gene.

22. The system according to any one of claims 1 to 2, wherein, taking into account (a) the number of sequence reads of the plurality of sequence reads having bases supporting the SMN1 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads having bases supporting the corresponding SMN2 gene-specific base, the most likely combination of the possible copy number of the SMN1 gene and the possible copy number of the SMN2 gene is associated with the highest posterior probability relative to other combinations in the plurality of combinations.

23. The system according to any one of claims 1 to 2, wherein determining the most likely combination of the possible copy number of the SMN1 gene and the possible combination of the SMN2 gene comprises: Determining the most likely combination among the plurality of possible combinations of the possible copy number of the SMN1 gene and the possible copy number of the SMN2 gene, each including a total equal to the determined copy number of any full-length SMN gene, taking into account the ratio of (a) the number of sequence reads of the plurality of sequence reads having bases supporting the SMN1 gene-specific base to (b) the number of sequence reads of the plurality of sequence reads having bases supporting the SMN2 gene-specific base corresponding to the SMN1 gene-specific base.

24. The system according to any one of claims 1 to 2, wherein determining the most likely combination of the possible copy number of the SMN1 gene and the possible combination of the SMN2 gene comprises: Determine (a) the number of sequence reads of the plurality of sequence reads having bases that support the SMN1 gene-specific bases and (b) the number of sequence reads of the plurality of sequence reads having bases that support the SMN2 gene-specific bases corresponding to the SMN1 gene-specific bases of the SMN2 gene; Determine the ratio of (a) the number of sequence reads of the plurality of sequence reads having bases that support the SMN1 gene-specific bases to (b) the number of sequence reads of the plurality of sequence reads having bases that support the SMN2 gene-specific bases corresponding to the SMN1 gene-specific bases of the SMN2 gene; and Based on the ratio of (a) the number of sequence reads of the plurality of sequence reads having bases that support the SMN1 gene-specific bases to (b) the number of sequence reads of the plurality of sequence reads having bases that support the SMN2 gene-specific bases corresponding to the SMN1 gene-specific bases of the SMN2 gene, determine the most likely combination among the plurality of possible combinations of the possible copy number of the SMN1 gene and the possible copy number of the SMN2 gene, each including a total copy number of any complete SMN gene determined to be; 25. The system according to any one of claims 1 to 2, The most likely combinations for determining the possible copy number of the SMN1 gene and the possible combinations of the SMN2 gene include: For each of the plurality of SMN1 gene-specific bases, taking into account (a) the number of sequence reads of the plurality of sequence reads having bases that support the SMN1 gene-specific bases and (b) the number of sequence reads of the plurality of sequence reads having bases that support the SMN2 gene-specific bases corresponding to the SMN1 gene-specific bases of the SMN2 gene, determine the most likely combination associated with the highest posterior probability among the plurality of possible combinations of the possible copy number of the SMN1 gene and the possible copy number of the SMN2 gene, each including a total copy number of any complete SMN gene determined to be, and wherein determining the copy number of the SMN1 gene includes: determining the copy number of the SMN1 gene based on the possible copy number of the SMN1 gene in the most likely combination of the possible copy number of the SMN1 gene and the possible copy number of the SMN2 gene determined for each of the plurality of SMN1 gene-specific bases.

26. The system according to claim 25, wherein the SMN1 gene-specific base has identity with each of the plurality of SMN1 gene-specific bases other than the SMN1 gene-specific bases exceeding a pre-determined identity threshold.

27. The system according to claim 26, wherein the identity threshold is 97%.

28. The system according to claim 25, wherein the plurality of SMN1 gene-specific bases includes 8 SMN1 gene-specific bases.

29. The system according to claim 25, wherein each of the plurality of SMN1 gene-specific bases can be located on intron 6, exon 7, intron 7, or exon 8 of the SMN1 gene.

30. The system according to claim 25, wherein if the subject is of a first race, the plurality of SMN1 gene-specific bases are different, if the subject is of a second race, the plurality of SMN1 gene-specific bases are different, and if the subject is of an unknown race, the plurality of SMN1 gene-specific bases are different.

31. The system according to claim 25, wherein the race of the subject is unknown, and wherein the plurality of SMN1 gene-specific bases are not race-specific.

32. The system according to claim 25, wherein the race of the subject is known, and wherein the plurality of SMN1 gene-specific bases are specific to the race of the subject.

33. The system according to claim 25, wherein the hardware processor is programmed by the executable instructions to perform: receiving race information of the subject; and selecting the plurality of SMN1 gene-specific bases from the plurality of SMN1 gene-specific bases based on the received race information.

34. The system according to any one of claims 1 to 2, wherein determining the copy number of the SMN1 gene comprises: Determining the copy number of the SMN1 gene and the copy number of the SMN2 gene using the most likely combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene determined for each of the plurality of SMN1 gene-specific bases.

35. The system according to any one of claims 1 to 2, wherein determining the copy number of the SMN1 gene comprises: Determining the copy number of the SMN1 gene using the most likely combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene determined for the SMN1 gene-specific base and a second predetermined posterior probability threshold of the combination of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene.

36. The system according to claim 35, wherein the second predetermined posterior probability threshold is 0.6 or 0.

8.

37. The system according to claim 25, wherein the possible copy numbers of most of the determined SMN1 genes are consistent, and wherein the copy number of the SMN1 gene determined is the consistent possible copy number of the SMN1 gene.

38. The system according to claim 37, wherein the hardware processor is programmed by the executable instructions to perform: Considering (a) the number of sequence reads of the plurality of sequence reads having bases supporting any one of the plurality of SMN1 gene-specific bases and (b) the number of sequence reads of the plurality of sequence reads having bases supporting any one of the plurality of corresponding SMN2 gene-specific bases, determining the possible combinations of the possible copy numbers of the SMN1 gene and the possible copy numbers of the SMN2 gene including the copy numbers of any complete SMN genes determined in total; and Determining the possible copy number of the possible combination as the consistent possible copy number of the SMN1 gene.

39. The system according to any one of claims 1 to 2, wherein determining the copy number of the SMN1 gene comprises determining that the copy number of the SMN1 gene is zero, one, or more than one.

40. The system according to any one of claims 1 to 2, wherein the hardware processor is programmed by the executable instructions to perform: determining the spinal muscular atrophy SMA status of the subject based on the copy number of the SMN1 gene.

41. The system according to claim 40, wherein the SMA status of the subject comprises SMA, an SMA carrier but not SMA, and not an SMA carrier.

42. The system according to any one of claims 1 to 2, wherein the hardware processor is programmed by the executable instructions to perform: using the number of sequence reads of the plurality of sequence reads aligned with g.27134 of the SMN1 gene and the bases of the sequence reads aligned with g.27134 of the SMN1 gene to determine that the subject is a silent SMA carrier.

43. The system according to any one of claims 1 to 2, wherein the hardware processor is programmed by the executable instructions to perform: determining a treatment recommendation for the subject based on the determined copy number of the SMN1 gene.

44. The system according to claim 43, wherein the treatment recommendation comprises administering Nusinersen and / or Zolgensma to the subject.

45. A system for genotyping cytochrome P450 family 2 subfamily D member 6 CYP2D6 gene, comprising: A non-transitory memory configured to store executable instructions and sequence data, the sequence data comprising a plurality of sequence reads obtained from a sample of a subject and aligned with cytochrome P450 family 2 subfamily D member 6 CYP2D6 gene or cytochrome P450 family 2 subfamily D member 7 CYP2D7 gene; and A hardware processor communicatively coupled to the non-transitory memory, the hardware processor being programmed by the executable instructions to perform: Receiving sequence data comprising a plurality of sequence reads obtained from a sample of a subject and aligned with cytochrome P450 family 2 subfamily D member 6 CYP2D6 gene or cytochrome P450 family 2 subfamily D member 7 CYP2D7 gene; Determining (i) a first number of sequence reads of the plurality of sequence reads aligned with the CYP2D6 gene or the CYP2D7 gene; Respectively using (i) the length of the CYP2D6 gene or the CYP2D7 gene to determine (i) a first normalized number of the sequence reads aligned with the CYP2D6 gene or the CYP2D7 gene; Considering (i) the first normalized number of the sequence reads aligned with the CYP2D6 gene or the CYP2D7 gene, using a Gaussian mixture model comprising a plurality of Gaussian functions each representing a different integer copy number to determine (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene; For one of the multiple CYP2D6 gene - specific bases, determine the most likely combination among multiple possible combinations of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene, each including a total copy number that is the sum of the determined CYP2D6 gene and CYP2D7 gene copy numbers, taking into account (a) the number of sequence reads of the multiple sequence reads having bases that support the CYP2D6 gene - specific base and (b) the number of sequence reads of the multiple sequence reads having bases that support the CYP2D7 gene - specific base corresponding to the CYP2D6 gene - specific base; and Use the most likely combination of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene determined for the CYP2D6 gene - specific base to determine the alleles of the CYP2D6 gene that the subject has.

46. The system according to claim 45, wherein the sequence data comprises whole - genome sequencing (WGS) data or short - read WGS data.

47. The system according to any one of claims 45 to 46, wherein the subject is a neonatal subject, a pediatric subject, an adolescent subject, or an adult subject.

48. The system according to any one of claims 45 to 46, wherein the sample comprises cells or cell - free DNA.

49. The system according to any one of claims 45 to 46, wherein the sequence reads of the multiple sequence reads are aligned with the CYP2D6 gene or the CYP2D7 gene, and the alignment quality score is approximately zero.

50. The system according to any one of claims 45 to 46, wherein determining (i) a first number of sequence reads of the plurality of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene comprises: Determine (i) a first number of the sequence reads of the multiple sequence reads that align with at least one exon or intron of the CYP2D6 gene or at least one exon or intron of the CYP2D7 gene.

51. The system according to any one of claims 45 to 46, wherein determining (i) a first normalized quantity of the sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene comprises: Respectively use (i) the length of the CYP2D6 gene or the CYP2D7 gene to determine (i) a first normalized number of the sequence reads that align with the CYP2D6 gene or the CYP2D7 gene, and determine (iii) the depth of the sequence reads in the region of the subject's genome other than the locus containing the CYP2D6 gene and the CYP2D7 gene in the sequence data.

52. The system according to claim 51, wherein determining (i) the first normalized number of the sequence reads that align with the CYP2D6 gene or the CYP2D7 gene comprises: Respectively use (i) the length of the CYP2D6 gene or the CYP2D7 gene to determine (i) a first CYP2D6 - gene - or CYP2D7 - gene - length - normalized number of the sequence reads that align with the CYP2D6 gene or the CYP2D7 gene; and The depth of sequence reads using regions of the genome of the subject other than those containing the CYP2D6 gene and the locus of the CYP2D7 gene is determined according to the number of (i) CYP2D6 gene or CYP2D7 gene length normalizations to determine (i) the first normalized depth of the sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene. The first normalized depth of the sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene is the first normalized number of the sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, respectively.

53. The system according to any one of claims 45 to 46, wherein determining (i) the first normalized quantity of the sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene comprises: Use (i) the GC content of the CYP2D6 gene or the CYP2D7 gene to determine (i) the first normalized number of the sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, and to determine (iii) the depth of sequence reads in regions of the genome of the subject other than those containing the CYP2D6 gene and the locus of the CYP2D7 gene in the sequence data, and (iv) to determine the GC content of the regions of the genome.

54. The system according to claim 51, wherein the depth of the region includes the average depth or median depth of sequence reads in regions of the genome of the subject other than those containing the CYP2D6 gene and the locus of the CYP2D7 gene in the sequence data.

55. The system according to claim 54, wherein the region contains approximately 3000 preselected regions of the genome of the subject, each with a length of approximately 2 kb and each spanning the genome of the subject.

56. The system according to any one of claims 45 to 46, wherein (i) the first normalized number of the sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene is from approximately 30 to approximately 40.

57. The system according to any one of claims 45 to 46, wherein the Gaussian mixture model includes a one-dimensional Gaussian mixture model.

58. The system according to any one of claims 45 to 46, wherein the plurality of Gaussian functions of the Gaussian mixture model represent integer copy numbers from 0 to 10.

59. The system according to any one of claims 45 to 46, wherein the mean of each Gaussian function in the plurality of Gaussian functions is the integer copy number represented by the Gaussian function.

60. The system according to any one of claims 45 to 46, wherein determining (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene comprises: Taking into account (i) the first normalized number of the sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, use the Gaussian mixture model and a first predetermined posterior probability threshold to determine (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene.

61. The system according to claim 60, wherein the first predetermined posterior probability threshold is 0.

95.

62. The system according to any one of claims 45 to 46, wherein the most likely combination of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene is associated with the highest posterior probability, taking into account (a) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D6 gene-specific bases and (b) the number of sequence reads of the plurality of sequence reads having bases supporting the corresponding CYP2D7 gene-specific bases.

63. The system according to any one of claims 45 to 46, wherein the most likely combination for determining the likely copy number of the CYP2D6 gene and the likely copy number of the CYP2D7 gene comprises: Determine the most likely combination among the plurality of possible combinations of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene, each including a total copy number of the determined CYP2D6 gene and CYP2D7 gene, based on the ratio of (a) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D6 gene-specific bases to (b) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D7 gene-specific bases corresponding to the CYP2D6 gene-specific bases.

64. The system according to any one of claims 45 to 46, wherein determining the most likely combination of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene comprises: Determining (a) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D6 gene-specific bases and (b) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D7 gene-specific bases corresponding to the CYP2D6 gene-specific bases; Determining the ratio of (a) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D6 gene-specific bases to (b) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D7 gene-specific bases corresponding to the CYP2D6 gene-specific bases; And Determining the most likely combination among the plurality of possible combinations of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene, each including a total copy number of the determined CYP2D6 gene and CYP2D7 gene, taking into account the ratio of (a) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D6 gene-specific bases to (b) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D7 gene-specific bases corresponding to the CYP2D6 gene-specific bases.

65. The system according to any one of claims 45 to 46, wherein determining the alleles of the CYP2D6 gene that the subject has comprises: Use the most likely combination of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene determined for the CYP2D6 gene-specific bases to determine one or more structural variants of the CYP2D6 gene possessed by the subject.

66. The system according to claim 65, The most likely combinations for determining the possible copy number of the CYP2D6 gene and the possible copy number of the CYP2D7 gene include: For each of the plurality of CYP2D6 gene-specific bases, determining the most likely combination associated with the highest posterior probability among a plurality of possible combinations of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene, each including a total number of copies of the determined CYP2D6 gene and CYP2D7 gene, taking into account (a) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D6 gene-specific base and (b) the number of sequence reads of the plurality of sequence reads having bases supporting the CYP2D7 gene-specific base corresponding to the CYP2D6 gene-specific base of the CYP2D7 gene, and where determining the one or more structural variants of the CYP2D6 gene that the subject has includes: using the most likely combination of the possible copy numbers of the CYP2D6 gene and the possible copy numbers of the CYP2D7 gene determined for each of the plurality of CYP2D6 gene-specific bases to determine the one or more structural variants of the CYP2D6 gene that the subject has.

67. The system according to claim 66, wherein determining the one or more structural variants of the CYP2D6 gene that the subject has comprises: Determining the one or more structural variants of the CYP2D6 gene that the subject has based on the copy numbers of the CYP2D6 gene of the most likely combinations determined for two or more different ones of the plurality of CYP2D6 gene-specific bases and the positions of the two or more CYP2D6 gene-specific bases.

68. The system according to claim 66, wherein the CYP2D6 gene-specific base is consistent with each of the plurality of CYP2D6 gene-specific bases other than the CYP2D6 gene-specific base exceeding a predetermined consistency threshold.

69. The system according to claim 68, wherein the consistency threshold is 97%.

70. The system according to claim 66, wherein the plurality of CYP2D6 gene-specific bases includes 118 CYP2D6 gene-specific bases.

71. The system according to claim 66, wherein if the subject is of a first race, the plurality of CYP2D6 gene-specific bases are different, if the subject is of a second race, the plurality of CYP2D6 gene-specific bases are different, and if the subject is of an unknown race, the plurality of CYP2D6 gene-specific bases are different.

72. The system according to claim 66, wherein the race of the subject is unknown and wherein the plurality of CYP2D6 gene-specific bases are not race-specific.

73. The system according to claim 66, wherein the race of the subject is known and wherein the plurality of CYP2D6 gene-specific bases are specific to the race of the subject.

74. The system according to claim 66, wherein the hardware processor is programmed by the executable instructions to perform: Receiving race information of the subject; and Select the plurality of CYP2D6 gene-specific bases based on the received ethnicity information.

75. The system according to claim 65, wherein the hardware processor is programmed by the executable instructions to perform: Determine a second quantity of the sequence reads of the plurality of sequence reads aligned to the spacer region between (ii) the CYP2D7 gene and the repetitive element REP7 downstream of the CYP2D7 gene; Use the length of (ii) the spacer region to determine a second normalized quantity of the sequence reads aligned to the spacer region; and Using the Gaussian mixture model, determine the copy number of (ii) the spacer region in consideration of the second normalized quantity of the sequence reads aligned to the spacer region; Determining the structural variant of the CYP2D6 gene that the subject has includes: Determine the alleles of the CYP2D6 gene that the subject has using the most likely combination of the possible copy number of the CYP2D6 gene, the possible copy number of the CYP2D7 gene, and the copy number of the spacer region determined for the CYP2D6 gene-specific bases.

76. The system according to claim 75, wherein the one or more structural variants comprise a CYP2D6 / CYP2D7 fusion allele having the spacer region and a repetitive element REP7 downstream of the CYP2D6 / CYP2D7 fusion allele.

77. The system according to any one of claims 45 to 46, wherein the hardware processor is programmed by the executable instructions to perform: determining one or more minor variants of the CYP2D6 gene that the subject has using the received sequence data.

78. The system according to claim 77, wherein determining the one or more minor variants of the CYP2D6 gene that the subject has comprises: For a minor variant position of the CYP2D6 gene associated with a minor variant allele of the CYP2D6 gene, determine the most likely combination of the possible copy number of the minor variant allele of the CYP2D6 gene and the possible copy number of the reference allele of the CYP2D6 gene at the minor variant position that together count as the copy number of the CYP2D6 gene at the minor variant position, considering (a) the quantity of sequence reads having bases supporting the minor variant allele of the CYP2D6 gene at the minor variant position and (b) the quantity of sequence reads having bases supporting the reference allele of the CYP2D6 gene at the minor variant position, wherein the possible copy number of the minor variant allele of the CYP2D6 gene in the most likely combination indicates the one or more minor variants of the CYP2D6 gene.

79. The system according to claim 77, wherein determining the one or more minor variants of the CYP2D6 gene that the subject has comprises: For each of a plurality of variant positions of the CYP2D6 gene, where the variant position is associated with a variant allele of the CYP2D6 gene, determine the most likely combination of the possible copy numbers of the variant allele of the CYP2D6 gene at the variant position and the possible copy numbers of the reference allele of the CYP2D6 gene at the variant position that together amount to the copy number of the CYP2D6 gene at the variant position, taking into account (a) the number of sequence reads having a base supporting the variant allele of the CYP2D6 gene at the variant position and (b) the number of sequence reads having a base supporting the reference allele of the CYP2D6 gene at the variant position, wherein the possible copy numbers of the variant alleles of the CYP2D6 gene at the plurality of variant positions of the most likely combination indicate the one or more variants of the CYP2D6 gene.

80. The system according to any one of claims 45 to 46, wherein the hardware processor is programmed by the executable instructions to perform: For a variant position of the CYP2D6 gene that is associated with a variant allele of the CYP2D6 gene, determine the most likely combination of the possible copy numbers of the variant allele of the CYP2D6 gene at the variant position and the possible copy numbers of the reference allele of the CYP2D6 gene at the variant position that together amount to the copy number of the CYP2D6 gene at the variant position, taking into account (a) the number of sequence reads that align with the CYP2D6 gene and overlap the variant position and have a base supporting the variant allele of the CYP2D6 gene at the variant position and (b) the number of sequence reads that align with the CYP2D6 gene and overlap the variant position and have a base supporting the reference allele of the CYP2D6 gene at the variant position; and Use the determined possible copy numbers of the variant alleles of the CYP2D6 gene of the most likely combination to determine one or more variants of the CYP2D6 gene.

81. The system according to any one of claims 45 to 46, wherein the hardware processor is programmed by the executable instructions to perform: For each of the multiple minor variant positions of the CYP2D6 gene, where the minor variant position is associated with a minor variant allele of the CYP2D6 gene, determine the most likely combination of the possible copy numbers of the minor variant allele of the CYP2D6 gene at the minor variant position and the possible copy numbers of the reference allele of the CYP2D6 gene at the minor variant position that together amount to the copy number of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads that align with the CYP2D6 gene and overlap the minor variant position and have bases supporting the minor variant allele of the CYP2D6 gene at the minor variant position and (b) the number of sequence reads that align with the CYP2D6 gene and overlap the minor variant position and have bases supporting the reference allele of the CYP2D6 gene at the minor variant position; and Use the possible copy numbers of the minor variant alleles of the CYP2D6 gene at the multiple minor variant positions in the determined most likely combination to determine one or more minor variants of the CYP2D6 gene.

82. The system according to claim 78, wherein the minor variant position is in the CYP2D6 / CYP2D7 homology region, and wherein determining the most likely combination includes determining the most likely combination of the possible copy numbers of the minor variant allele of the CYP2D6 gene at the minor variant position and the possible copy numbers of the reference allele of the CYP2D6 gene at the minor variant position that together amount to the copy number of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads that align with the CYP2D6 gene or the CYP2D7 gene and have bases supporting the minor variant allele of the CYP2D6 gene at the minor variant position and / or (b) the number of sequence reads that align with the CYP2D6 gene or the CYP2D7 gene and have bases supporting the reference allele of the CYP2D6 gene at the minor variant position.

83. The system according to claim 78, wherein the minor variant position is not in the CYP2D6 / CYP2D7 homology region, and wherein determining the most likely combination includes determining the most likely combination of the possible copy numbers of the minor variant allele of the CYP2D6 gene at the minor variant position and the possible copy numbers of the reference allele of the CYP2D6 gene at the minor variant position that together amount to the copy number of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads that align with the CYP2D6 gene and not the CYP2D7 gene and have bases supporting the minor variant allele of the CYP2D6 gene at the minor variant position and / or (b) the number of sequence reads that align with the CYP2D6 gene and not the CYP2D7 gene and have bases supporting the reference allele of the CYP2D6 gene at the minor variant position.

84. The system according to claim 78, comprising determining the copy number of the CYP2D6 gene at the minor variant position.

85. The system according to claim 78, wherein the copy number of the CYP2D6 gene at the minor variant position comprises the copy number of the CYP2D6 gene.

86. The system according to claim 78, wherein the copy number of the CYP2D6 gene at the minor variant position comprises the copy number of the CYP2D6 gene of the most likely combination of the possible copy numbers of the CYP2D6 gene determined.

87. The system according to claim 78, wherein the copy number of the CYP2D6 gene at the minor variant position comprises the copy number of the CYP2D6 gene of the most likely combination and closest to the minor variant position of the possible copy numbers of the CYP2D6 gene.

88. The system according to claim 78, wherein the copy number of the CYP2D6 gene at the minor variant position comprises the copy number of the CYP2D6 gene at the 5' position or 3' position of the minor variant position.

89. The system according to claim 78, wherein the hardware processor is programmed by the executable instructions to perform: (a) determining the number of sequence reads of bases supporting the minor variant allele of the CYP2D6 gene; and (b) determining the number of sequence reads of bases supporting the reference allele of the CYP2D6 gene.

90. The system according to any one of claims 45 to 46, wherein determining the alleles of the CYP2D6 gene that the subject has comprises: Determine the allele of the CYP2D6 gene that the subject has.

91. The system according to any one of claims 45 to 46, wherein determining the alleles of the CYP2D6 gene that the subject has comprises: Use the one or more structural variants and / or the one or more minor variants of the determined CYP2D6 gene to determine the star allele and / or haplotype of the CYP2D6 gene that the subject has, optionally wherein the star allele is associated with a known function.

92. The system according to any one of claims 45 to 46, wherein the hardware processor is programmed by the executable instructions to perform: using the determined allele of the CYP2D6 gene to determine the level of CYP2D6 enzyme activity of the subject.

93. The system according to claim 92, wherein the enzyme activity is poor, moderate, normal or ultra-strong.

94. The system according to any one of claims 45 to 46, wherein the hardware processor is programmed by the executable instructions to perform: determining a treatment dose recommendation and / or treatment recommendation for the subject based on the allele of the CYP2D6 gene that the subject has.

95. A system for paralog genotyping, comprising: A non-transitory memory configured to store executable instructions and sequence data, the sequence data comprising a plurality of sequence reads obtained from a sample of a subject and aligned with a first paralog or a second paralog; and A hardware processor in communication with the non-transitory memory, the hardware processor being programmed by the executable instructions to perform: Considering (i) the first quantity of sequence reads aligned to the first region, a Gaussian mixture model using a plurality of Gaussian functions each representing a different integer copy number is used to determine the copy number of the first type of paralog. For one of the plurality of first paralog-specific bases, considering (a) the quantity of sequence reads of the plurality of sequence reads having bases supporting the first paralog-specific base and (b) the quantity of sequence reads of the plurality of sequence reads having bases supporting the second paralog-specific base corresponding to the first paralog-specific base of the second paralog, the most likely combination among a plurality of possible combinations respectively including the possible copy numbers of the first type of first paralog and the possible copy numbers of the first type of second paralog totaling the copy number of the first type of paralog determined is determined; and The most likely combination of the possible copy numbers of the first paralog and the possible copy numbers of the second paralog determined for the first paralog-specific base is used to determine the copy number or allele of the first paralog.

96. The system according to claim 95, wherein the hardware processor is programmed by the executable instructions to perform: determining (i) the first quantity of sequence reads of the plurality of sequence reads obtained from a sample of a subject in the sequence data and aligned to the first region.

97. The system according to any one of claims 95 to 96, wherein the hardware processor is programmed by the executable instructions to perform: using (i) the length of the first region to determine (i) a first normalized number of the sequence reads aligned to the first region, wherein determining the copy number of the first type of paralog comprises: Considering (i) the first normalized quantity of the sequence reads aligned to the first region, the Gaussian mixture model is used to determine the copy number of the first type of paralog.

98. The system according to any one of claims 95 to 96, wherein the hardware processor is programmed by the executable instructions to perform: receiving the sequence data including the plurality of sequence reads aligned to the first region.

99. The system according to any one of claims 95 to 96, wherein the hardware processor is programmed by the executable instructions to perform: considering (ii) the second quantity of sequence reads aligned to a second region, the Gaussian mixture is used to determine the copy number of one or more paralogs of the second type.

100. The system according to claim 99, wherein determining the copy number or allele of the first paralog comprises: The most likely combination of the possible copy numbers of the first paralog and the possible copy numbers of the second paralog determined for the first paralog-specific base and the copy number of the one or more paralogs of the second type is used to determine the copy number or allele of the first paralog.

101. The system according to claim 99, wherein the hardware processor is programmed by the executable instructions to perform: determining the copy number of a third type of paralog from the copy number of the paralogs of the first type and the copy number of the paralogs of the second type, and wherein determining the copy number or allele of the first paralog comprises: The most likely combination of the possible copy numbers of the first paralog and the possible copy numbers of the second paralog determined for the first paralog-specific base is used to determine the copy number or allele of the first paralog.

102. The system according to claim 99, wherein the first paralog is the survival motor neuron 1 SMN1 gene, wherein the second paralog is the survival motor neuron 2 SMN2 gene, wherein the first region comprises at least exon 1 to exon 6 of the SMN1 gene and at least exon 1 to exon 6 of the SMN2 gene, wherein the second region comprises at least one of exon 7 and exon 8 of the SMN1 gene and at least one of exon 7 and exon 8 of the SMN2 gene, wherein the first type of the paralog comprises the complete SMN1 gene and the complete SMN2 gene, wherein the second type of the one or more paralogs comprises the complete SMN1 gene, the complete SMN2 gene, a truncated SMN1 gene or a truncated SMN2 gene, and wherein the copy number of the first paralog comprises the copy number of the SMN1 gene.

103. The system according to claim 99, wherein the first paralog is the cytochrome P450 family 2 subfamily D member 6 CYP2D6 gene, wherein the second paralog is the cytochrome P450 family 2 subfamily D member 7 CYP2D7 gene, wherein the first region comprises the CYP2D6 gene and the CYP2D7 gene, wherein the second region comprises the spacer region between the CYP2D7 gene and the repeat element REP7 downstream of the CYP2D7 gene, wherein the first type of the paralog comprises the CYP2D6 gene and the CYP2D7 gene, wherein the second type of the one or more paralogs comprises the CYP2D6 / CYP2D7 fusion allele having the spacer region and the repeat element REP7 downstream of the CYP2D6 / CYP2D7 fusion allele, and wherein the allele of the first paralog comprises the allele of the CYP2D6 gene, which is a minor variant or a structural variant of the CYP2D6 gene.

104. The system according to any one of claims 95 to 96, wherein the first paralog and the second paralog have at least 90% sequence identity.

Images