RU2807604C2 - Methods and systems for diagnostics according to whole genome sequencing data - Google Patents

Abstract

FIELD: biologic information technology.

SUBSTANCE: method for determining the copy number of the survival motor neuron 1 (SMN1) gene is described. Under the control of a hardware processor, sequencing data containing multiple sequence reads obtained from a sample from a subject aligned to the survival motor neuron gene 1 (SMN1) or survival motor neuron gene 2 (SMN2) is obtained. A first number of sequence reads is determined from a plurality of sequence reads aligned to the first region SMN1 or SMN2 containing at least 1 of the first 6 exons of the SMN1 gene or the SMN2 gene, respectively, and a second number of sequence reads from a plurality of sequence reads aligned to the second region SMN1 or SMN2 containing at least exon 7 or exon 8 of the SMN1 gene or the SMN2 gene. The first normalized number of sequence reads and the second normalized number of sequence reads are determined. The copy number of full-length survival motor neuron (SMN) genes and the copy number of any intact SMN genes are determined. For one set of SMN1 gene-specific bases associated with an intact SMN1 gene, the most likely combination is determined from a plurality of possible combinations, each containing a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene, which add up to the number of copies of any identified intact SMN genes. The number of copies of the SMN1 gene is determined using the most probable combination of the possible number of copies of the SMN1 gene and the possible number of copies of the SMN2 gene, determined for a specific base for the SMN1 gene.

EFFECT: invention makes it possible to determine the number of copies of the SMN1 gene.

45 cl, 32 dwg, 25 tbl, 2 ex

Description

CROSS-REFERENCES TO RELATED APPLICATIONS

[0001] This application claims the benefit of U.S. Provisional Patent Application No. 62/896548, filed September 5, 2019, U.S. Provisional Patent Application No. 62/908555, filed September 30, 2019, and U.S. Provisional Patent Application No. 63/006651, filed April 7, 2020. The contents of each of the related applications are incorporated herein by reference in their entirety.

State of the art

Field of technology

[0002] The present disclosure relates generally to the field of paralog genotyping and, in particular, to paralog genotyping using sequencing data.

State of the art

[0003] Genotyping is a complex task. For example, spinal muscular atrophy is caused by loss of functionality of the survival of motor neuron 1 gene ( SMN1 ), but retention of the paralogous gene SMN2 . Due to the nearly identical sequences of SMN1 and its paralog SMN2, analysis of this region is challenging. As another example, CYP2D6 is involved in the metabolism of 25% of all drugs. Genotyping of CYP2D6 is challenging due to its high polymorphism, the presence of common structural variants (SVs), and high sequence similarity to the pseudogene paralogue CYP2D7 .

disclosure of the invention

[0004] 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 controlled by a processor (such as a hardware processor or virtual processor) and includes: obtaining sequencing data comprising multiple sequence reads from a sample from a subject 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 from a plurality of sequence reads aligned to a first region of the SMN1 or SMN2 gene containing at least one of the first 6 exons of the SMN1 gene or the SMN2 gene, respectively, and (ii) a second number of sequence reads from a plurality of sequence reads aligned to a second region of SMN1 or SMN2 containing at least exon 7 or exon 8 of the SMN1 gene or the SMN2 gene, respectively. The method may include: determining (i) a first normalized read count of a sequence aligned to a first region of the SMN1 or SMN2 gene, and (ii) a second normalized read count of a sequence aligned to a second SMN1 or SMN2 region, using (i) the length of the first SMN1 region or SMN2 and (ii) the length of the second region SMN1 or SMN2 , respectively. The method may include: determining (i) the copy number of full-length survival motor neuron (SMN) genes, each of which is an intact SMN1 gene, an intact SMN2 gene, a truncated SMN1 gene, or a truncated SMN2 gene; and (ii) the copy number of any intact SMN genes, each of which is an intact SMN1 gene or an intact SMN2 gene, using a mixture normal distribution model containing a plurality of normal distributions, each representing a different integer copy number, subject to (i) the first a normalized read count of the sequence aligned to the first SMN1 or SMN2 region, and (ii) a second normalized read count of the sequence aligned to the second SMN1 or SMN2 region, respectively. The method may include: for one set of SMN1 gene-specific bases associated with an intact SMN1 gene, determining the most likely combination from a plurality of possible combinations, each of which contains a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene, summed with the number of copies of any determined intact SMN genes, considering (a) the number of sequence reads from a set of sequence reads with bases that contain a base specific to the SMN1 gene, and (b) the number of sequence reads from a set of sequence reads with bases that contain a base specific to the SMN2 gene , corresponding to a base specific to the SMN1 gene. The method may include: determining the copy number of the SMN1 gene using the most likely combination of the possible SMN1 gene copy number and the possible SMN2 gene copy number determined for the SMN1 gene-specific base.

[0005] In some embodiments, the sequencing data comprises whole genome sequencing data (whole genome sequencing, WGS) or short-read WGS data. In some embodiments, the subject is a fetal subject, a neonatal subject, a pediatric subject, an adolescent subject, or an adult subject. The sample may contain cells or extracellular DNA. The sample may contain fetal cells or extracellular fetal DNA.

[0006] In some embodiments, a sequence read from a plurality of sequence reads is aligned to a first region SMN1 or SMN2 or a second region SMN1 or SMN2 with an alignment quality score of approximately zero. The first region of SMN1 or SMN2 may contain exon 1 to 6 of the SMN1 gene or the SMN2 gene, respectively, and is approximately 22.2 kb in length. The second region of SMN1 or SMN2 may contain exon 7 and exon 8 of the SMN1 gene or the SMN2 gene, respectively, and be approximately 6 kb in length.

[0007] In some embodiments, determining (i) a first normalized read count of a sequence aligned to a first SMN1 or SMN2 region and (ii) a second normalized read count of a sequence aligned to a second region comprises: determining (i) a first normalized read count sequence aligned to the first SMN1 or SMN2 region, and (ii) a second normalized read count of the sequence aligned to the 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, and (iii) the sequence read depth of the subject's genome region, excluding genetic loci containing the SMN1 gene and the SMN2 gene, as determined by sequencing. Determining (i) a first normalized read count of a sequence aligned to a first SMN1 or SMN2 region and (ii) a second normalized read count of a sequence aligned to a second SMN1 or SMN2 region may include: determining (i) a normalized read count over the length of the SMN1 region or SMN2 sequence aligned to the first SMN1 or SMN2 region, and (ii) the normalized read count to the length of the SMN1 or SMN2 region of the sequence aligned to the second SMN1 or SMN2 region, using (i) the length of the first SMN1 or SMN2 region and (ii) the length of the second region SMN1 or SMN2 , respectively. Determining (i) a first normalized read count of a sequence aligned to a first SMN1 or SMN2 region, and (ii) a second normalized read count of a sequence aligned to a second SMN1 or SMN2 region may include: determining (i) a first normalized read depth of a sequence aligned with a first SMN1 or SMN2 region, and (ii) a second normalized read depth of a sequence aligned to a second SMN1 or SMN2 region, based on (i) a first normalized read count along the length of the SMN1 or SMN2 region and (ii) a second normalized read count along the length region of SMN1 or SMN2 , respectively, using the sequence read count of a region of the subject's genome other than the genetic loci containing the SMN1 gene and the SMN2 gene, the first normalized read count of the sequence aligned to the first SMN1 or SMN2 region, and the second normalized read count of the sequence aligned with a second region SMN1 or SMN2 , which represent the first normalized intensity and the second normalized intensity, respectively.

[0008] In some embodiments, determining (i) a first normalized read count of a sequence aligned to a first SMN1 or SMN2 region and (ii) a second normalized read count of a sequence aligned to a second region comprises: determining (i) a first normalized read count sequence aligned to the first SMN1 or SMN2 region, and (ii) a second normalized read count of the sequence aligned to the second SMN1 or SMN2 region, using (i) the GC content of the first SMN1 or SMN2 region and (ii) the GC content of the second region SMN1 or SMN2 , respectively, and (iii) the sequence read depth of a region of the subject's genome that is distinct from the genetic loci containing the SMN1 gene and the SMN2 gene as determined by sequencing, and (iv) the GC content of the genomic region.

[0009] In some embodiments, the region read depth is the average depth or median sequence read depth of a region of the subject's genome other than the genetic loci containing the SMN1 gene and the SMN2 gene, as determined by sequencing. The region may contain approximately 3000 preselected regions of approximately 2 kb in length. each in the subject's genome. In some embodiments, (i) the first normalized read count of the sequence aligned to the first SMN1 or SMN2 region and/or (ii) the second normalized read count of the sequence aligned to the second SMN1 or SMN2 region is from about 30 to about 40.

[0010] In some embodiments, the mixture normal distribution model is a one-dimensional mixture normal distribution model. The plurality of distributions of a mixture of normal distributions may represent integer copy numbers from 0 to 10. The mean for each of the plurality of normal distributions may represent the integer copy numbers represented by the normal distributions.

[0011] In some embodiments, determining (i) the copy number of all SMN genes and (ii) the copy number of any intact SMN genes includes determining (i) the copy number of all SMN genes and (ii) the copy number of any intact SMN genes using a mixture model normal distributions, and a first predetermined posterior probability threshold given (i) a first normalized read count of the sequence aligned to the first SMN1 or SMN2 region, and (ii) a second normalized read count of the sequence aligned to the second SMN1 or SMN2 region, respectively. The first predefined posterior probability threshold may be 0.95.

[0012] In some embodiments, the method includes: determining the copy number of truncated SMN genes using (i) a certain copy number of the total number of SMN genes and (ii) a certain copy number of intact SMN genes. The copy number of truncated SMN genes may represent the difference of (i) the total copy number of certain SMN genes and (ii) the specific copy number of intact SMN genes.

[0013] In some embodiments, the SMN1 gene-specific base is a splicing enhancer. The SMN1 gene-specific base may be the base at c.840 of the SMN1 gene. In some embodiments, the most likely combination of a possible SMN1 gene copy number and a possible SMN2 gene copy number is associated with the highest posterior probability, compared to other combinations of the plurality of combinations given (a) the number of sequence reads from the plurality of sequence reads with bases that contain a specific for the SMN1 gene base, and (b) the number of sequence reads from the plurality of sequence reads with bases that contain the corresponding SMN2 gene-specific base.

[0014] In some embodiments, determining the most likely combination of a possible SMN1 gene copy number and a possible SMN2 gene combination includes: determining the most likely combination, from a plurality of possible combinations, each of which contains a possible SMN1 gene copy number and a possible SMN2 gene copy number, summed with the copy number of any defined intact SMN genes, taking into account the ratio of (a) the number of sequence reads from the set of sequence reads with bases that contain the SMN1 gene-specific base, and (b) the number of sequence reads from the set of sequence reads with bases that contain SMN2 a gene-specific base of the SMN2 gene corresponding to a specific base for the SMN1 gene. Determining the most likely combination of a possible SMN1 gene copy number and a possible SMN2 gene combination may involve: determining (a) the number of sequence reads from a plurality of sequence reads with bases that contain a SMN1 gene-specific base, and (b) the number of sequence reads from a plurality of sequence reads with bases that contain the SMN2 gene-specific base of the SMN2 gene corresponding to the SMN1 gene-specific base; determining the ratio of (a) the number of sequence reads from a plurality of sequence reads with bases that contain an SMN1 gene-specific base, and (b) the number of sequence reads from a plurality of sequence reads with bases that contain an SMN2 gene-specific base corresponding to the SMN2 gene-specific SMN1 gene base; and determining the most likely combination, from a plurality of possible combinations, each of which contains a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene, summed with the number of copies of any intact SMN gene determined based on the ratio of (a) the number of sequence reads out of the plurality of sequence reads with bases that contain the SMN1 gene-specific base, and (b) the number of sequence reads from the plurality of sequence reads with bases that contain the SMN2 gene-specific base of the SMN2 gene corresponding to the SMN1 gene-specific base.

[0015] In some embodiments, determining the most likely combination of a possible number of copies of the SMN1 gene and a possible combination of the SMN2 gene includes: for each of a plurality of gene-specific SMN1 bases, determining the most likely combination from a plurality of possible combinations, each of which contains a possible number of copies of the gene SMN1 and the possible copy number of the SMN2 gene, summed with the copy number of any identified intact SMN genes that is associated with the largest posterior probability given by (a) the number of sequence reads from the set of sequence reads with bases that contain the SMN1 gene-specific base, and (b ) the number of sequence reads from a set of sequence reads with bases that contain the SMN2 gene-specific base of the SMN2 gene corresponding to the SMN1 gene-specific base. 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, 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 a plurality of SMN1 gene-specific bases.

[0016] In some embodiments, the SMN1 gene-specific base has a match with each of a plurality of SMN1 gene-specific bases other than the SMN1 gene-specific base above a predetermined match threshold. The compliance threshold may be 97%. The SMN1 gene-specific base set may contain 8 SMN1 gene-specific bases. Each of a variety of bases specific to the SMN1 gene may be located on intron 6, exon 7, intron 7, or exon 8 of the SMN1 gene. The base set specific to the SMN1 gene if the subject is a first race, the base set specific to the SMN1 gene if the subject is a second race, and the base set specific to the SMN1 gene if the subject is an unknown race may vary. The subject's race may be unknown, and many bases specific to the SMN1 gene may be non-race specific. The race of the subject may be known, and a plurality of bases specific to the SMN1 gene may be specific to the race of the subject. In some embodiments, the method includes obtaining information about the race of the subject. The method may include: selecting a plurality of bases specific for the SMN1 gene from the plurality of bases specific for the SMN1 gene based on the obtained race information.

[0017] In some embodiments, determining the SMN1 gene copy number includes: determining the SMN1 gene copy number and the SMN2 gene copy number using the most likely combination of the possible SMN1 gene copy number and the possible SMN2 gene copy number determined for each of a plurality of bases specific for SMN1 gene. The copy number determination may include: determining the copy number of the SMN1 gene using the most likely combination of the possible SMN1 gene copy number and the possible SMN2 gene copy number determined for an SMN1 gene-specific base and a second prespecified posterior probability threshold for the combination of the possible SMN1 gene copy number and possible number of copies of the SMN2 gene. The second predefined posterior probability threshold can be 0.6 or 0.8.

[0018] In some embodiments, the majority of the determined possible copy numbers of the SMN1 gene are the same. The copy number of a particular SMN1 gene may be consistent with the possible copy number of the SMN1 gene. The method may include: determining a candidate combination containing a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene, summed with the number of copies of any intact SMN gene, determined taking into account (a) the number of sequence reads from a plurality of sequence reads with bases that contain any of a plurality of bases specific to the SMN1 gene; and (b) a number of sequence reads from a plurality of sequence reads with bases that contain any of a plurality of corresponding bases specific to the SMN2 gene. The method may include: determining a possible copy number of a candidate combination that is a consensus possible copy number of the SMN1 gene.

[0019] In some embodiments, determining the copy number of the SMN1 gene includes: determining the copy number of the SMN1 gene to be zero, one, or greater than one. In some embodiments, the method includes: determining spinal muscular atrophy (SMA) status in a subject based on the copy number of the SMN1 gene. The SMA status for a subject may include SMA, SMA carrier/no SMA carrier, and non-SMA carrier. In some embodiments, the method includes identifying a subject as a silent carrier of SMA using a set of sequence reads from a plurality of sequence reads aligned to g.27134 of the SMN1 gene, and based on sequence reads aligned to g.27134 of the SMN1 gene.

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

[0021] Described herein are methods for genotyping the cytochrome P450 family 2 subfamily D member 6 gene ( CYP2D6 ). In some embodiments, a method for genotyping a CYP2D6 gene is controlled by a processor (such as a hardware processor or virtual processor) and includes: obtaining sequencing data comprising multiple sequence reads from a sample from a subject aligned to a CYP2D6 gene or a cytochrome P450 family 2 subfamily D member 7 gene ( CYP2D7 ). The method may include: determining (i) a first number of sequence reads from a plurality of sequence reads aligned to a CYP2D6 gene or a CYP2D7 gene. The method may include: determining (i) the first normalized read count of the sequence aligned to the CYP2D6 gene or the CYP2D7 gene, using (i) the length of the CYP2D6 gene or the CYP2D7 gene, respectively. The method may include: determining (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene using a mixture of normal distributions model containing a plurality of normal distributions, each representing a different integer copy number, based on (i) the first normalized read count of the sequence aligned to CYP2D6 genome or CYP2D7 genome. The method may include: for one of a plurality of bases specific to the CYP2D6 gene, determining the most likely combination of a plurality of possible combinations, each of which contains a possible number of copies of the CYP2D6 gene and a possible number of copies of the CYP2D7 gene, summed to the total number of copies of the CYP2D6 gene and the particular CYP2D7 gene, considering (a) multiple sequence reads with bases that contain a CYP2D6 gene-specific base, and (b) multiple sequence reads with bases that contain a CYP2D7 gene-specific base corresponding to a CYP2D6 gene-specific base. The method may include: determining the allele of the CYP2D6 gene present in the subject using the most likely combination of the possible number of copies of the CYP2D6 gene and the possible number of copies of the CYP2D7 gene determined for the base specific for the CYP2D6 gene.

[0022] In some embodiments, the sequencing data comprises whole genome sequencing (WGS) data or short-read WGS data. The subject may be a fetal subject, a neonatal subject, a pediatric subject, an adolescent subject, or an adult subject. The sample may contain cells or extracellular DNA. The sample may contain cells or extracellular DNA.

[0023] In some embodiments, a sequence read from a plurality of sequence reads aligned to a CYP2D6 gene or a CYP2D7 gene with an alignment quality score of about zero. In some embodiments, determining (i) a first number of sequence reads from a plurality of sequence reads aligned to a CYP2D6 gene or a CYP2D7 gene comprises: determining (i) a first number of sequence reads from a plurality of sequence reads aligned to at least one exon or intron of the gene CYP2D6 or at least one of the exons or introns of the CYP2D7 gene.

[0024] In some embodiments, determining (i) a first normalized count of sequence reads aligned to a CYP2D6 gene or a CYP2D7 gene comprises: determining (i) a first normalized count of sequence reads aligned to a CYP2D6 gene or a CYP2D7 gene using (i) the length of the CYP2D6 gene or the CYP2D7 gene, respectively, and (iii) the sequence read depth of a region of the subject's genome other than the genetic loci containing the CYP2D6 gene and the CYP2D7 gene as determined by sequencing. Determining (i) a first normalized count of sequence reads aligned to a CYP2D6 gene or a CYP2D7 gene, and (ii) a second normalized count of sequence reads aligned to a second region may include: determining (i) a first CYP2D6 gene or a length-normalized count of a CYP2D7 gene sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, using (i) the length of the CYP2D6 gene or the CYP2D7 gene, respectively. Determining (i) the first normalized depth of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, and (ii) the second normalized depth of sequence reads aligned to the second region may include: determining (i) the first normalized depth of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, from (i) the normalized count along the length of the CYP2D6 gene or the CYP2D7 gene, using the sequence read depth of a region of the subject's genome other than the genetic loci containing the CYP2D6 gene and the CYP2D7 gene, the first normalized read depth of the sequence aligned to the CYP2D6 gene or gene CYP2D7 , represents the first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, respectively.

[0025] In some embodiments, determining (i) a first normalized count of sequence reads aligned to a CYP2D6 gene or a CYP2D7 gene comprises: determining (i) a first normalized count of sequence reads aligned to a CYP2D6 gene or a CYP2D7 gene using (i) the GC content of the CYP2D6 gene or the CYP2D7 gene and (iii) the sequence read depth of the subject's genomic region other than the genetic loci containing the CYP2D6 gene and the CYP2D7 gene as determined by sequencing data and (iv) the GC content of the genomic region. The region read depth may include the average depth or median sequence read depth of a region of the subject's genome other than the genetic loci containing the CYP2D6 gene and the CYP2D7 gene as determined by sequencing. The region may contain approximately 3000 preselected regions of approximately 2 kb in length. each in the subject's genome. In some embodiments, (i) the first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene and/or (ii) the second normalized number of sequence reads aligned to the second region is from about 30 to about 40.

[0026] In some embodiments, the mixture normal distribution model is a one-dimensional mixture normal distribution model. The plurality of distributions of a mixture of normal distributions may represent integer copy numbers from 0 to 10. The mean for each of the plurality of normal distributions may represent the integer copy numbers represented by the normal distributions.

[0027] In some embodiments, determining (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene includes determining (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene using a mixture of normal distributions model and a first prespecified posterior probability threshold, provided that (i) the first normalized sequence read count matches the CYP2D6 gene or the CYP2D7 gene. The first predefined posterior probability threshold may be 0.95.

[0028] In some embodiments, the most probable combination of a possible CYP2D6 gene copy number and a possible CYP2D7 gene copy number is associated with the highest posterior probability, compared to other combinations of the plurality of combinations given (a) the number of sequence reads from the plurality of sequence reads with bases, that contain a CYP2D6 gene-specific base, and (b) the number of sequence reads from a plurality of sequence reads with bases that contain the corresponding CYP2D7 gene-specific base.

[0029] In some embodiments, determining the most likely combination containing a possible number of copies of the CYP2D6 gene and a possible number of copies of the CYP2D7 gene comprises: determining the most likely combination, from a plurality of possible combinations, each of which contains a possible number of copies of the CYP2D6 gene and a possible number of copies of the CYP2D7 gene, summed with the total copy number of the CYP2D6 gene and the specific CYP2D7 gene, taking into account the ratio of (a) the number of sequence reads from multiple sequence reads with bases that contain a CYP2D6 gene-specific base, and (b) the number of sequence reads from multiple sequence reads with bases that contain a CYP2D7 gene-specific base corresponding to a CYP2D6 gene-specific base. Determining the most likely combination containing the possible copy number of the CYP2D6 gene and the possible copy number may involve: determining (a) the number of sequence reads from a set of sequence reads with bases that contain a CYP2D6 gene-specific base, and (b) the number of sequence reads from a set sequence reads with bases that contain a CYP2D7 gene-specific base corresponding to a CYP2D6 gene-specific base; determining the ratio of (a) the number of sequence reads from a plurality of sequence reads with bases that contain a CYP2D6 gene-specific base, and (b) the number of sequence reads from a plurality of sequence reads with bases that contain a CYP2D7 gene-specific base corresponding to a CYP2D6 gene-specific base basis; and determining the most likely combination of a plurality of possible combinations, each containing a possible number of copies of the CYP2D6 gene and a possible number of copies of the CYP2D7 gene, summed with the total number of copies of the CYP2D6 gene and the particular CYP2D7 gene, taking into account the ratio of (a) the number of sequence reads from the plurality of sequence reads with bases that contain a CYP2D6 gene-specific base, and (b) the number of sequence reads from a plurality of sequence reads with bases that contain a CYP2D7 gene-specific base corresponding to a CYP2D6 gene-specific base.

[0030] In some embodiments, determining the allele of the CYP2D6 gene present in a subject includes: determining one or more structural variants of the CYP2D6 gene in the subject 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 specific gene CYP2D6 bases. In some embodiments, determining the most likely combination of a possible number of copies of a CYP2D6 gene and a possible number of copies of a CYP2D7 gene includes determining the most likely combination for each of a plurality of CYP2D6 gene-specific bases, a plurality of possible combinations, each of which contains a possible number of copies of a CYP2D6 gene and a possible number of copies of the CYP2D7 gene, summed with the total copy number of a particular CYP2D6 gene and the CYP2D7 gene, which is associated with the largest posterior probability given by (a) the number of sequence reads from the set of sequence reads with bases that contain the CYP2D6 gene-specific base, and (b) the number sequence reads from a plurality of sequence reads with bases that contain a CYP2D7 gene-specific base of the CYP2D7 gene corresponding to a CYP2D6 gene-specific base. Determining one or more structural variants of the CYP2D6 gene present in a subject may involve determining one or more structural variants 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 each of a plurality of CYP2D6 gene-specific bases. In some embodiments, determining one or more structural variants of the CYP2D6 gene present in a subject includes: determining one or more structural variants of the CYP2D6 gene in the subject, based on the copy number of the CYP2D6 gene, the most likely combinations determined for two or more of a plurality of different bases, specific for the CYP2D6 gene, and positions of two or more bases specific for the CYP2D6 gene.

[0031] In some embodiments, the CYP2D6 gene-specific base has a match with each of a plurality of CYP2D6 gene-specific bases, except for the CYP2D6 gene-specific base, above a predetermined matching threshold. The compliance threshold may be 97%. The CYP2D6 gene-specific base set may contain 118 CYP2D6 gene-specific bases. The base set specific to the CYP2D6 gene if the subject is a first race, the base set specific to the CYP2D6 gene if the subject is a second race, and the base set specific to the CYP2D6 gene if the subject is an unknown race may vary. The subject's race may be unknown, and many bases specific to the CYP2D6 gene may not be race specific. The race of the subject may be known, and many bases specific to the CYP2D6 gene may be specific to the race of the subject. In some embodiments, the method includes obtaining information about the race of the subject. The method may include: selecting a plurality of CYP2D6 gene-specific bases from the plurality of CYP2D6 gene-specific bases based on the obtained race information.

[0032] In some embodiments, the method includes: determining (ii) a second number of sequence reads from a plurality of sequence reads aligned to the spacer region between the CYP2D7 gene and the REP7 repeat element downstream of the CYP2D7 gene. The method may include: determining (ii) a second normalized read count of the sequence aligned to the spacer region using (ii) the length of the spacer region. The method may include: determining (ii) the copy number of the spacer region using a mixture normal distribution model defined by (ii) a second normalized number of sequence reads aligned to the spacer region. Determining one or more structural variants of the CYP2D6 gene present in a subject may include: determining one or more structural variants of the CYP2D6 gene in the subject using the most likely combination of possible CYP2D6 gene copy number and possible CYP2D7 gene copy number determined for a CYP2D6 gene-specific base , and the number of copies of the spacer region. One or more structural variants may contain a CYP2D6 / CYP2D7 fusion allele with a spacer region and a REP7 repeat element downstream of the CYP2D6 / CYP2D7 fusion allele.

[0033] In some embodiments, the method includes: identifying one or more small variants of the CYP2D6 gene in a subject using the obtained sequencing data. In some embodiments, determining one or more minor variants of the CYP2D6 gene present in a subject includes: for the minor variant allele position of the CYP2D6 gene associated with the minor variant allele of the CYP2D6 gene, determining the most likely combination of the possible number of copies of the minor variant allele of the CYP2D6 gene at the position of the minor variant and the possible number of copies of the reference allele of the CYP2D6 gene, summed with the number of copies of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads with bases that contain the small variant allele of the CYP2D6 gene at the minor variant position, and (b) the number of sequence reads with bases that contain the reference allele of the CYP2D6 gene at the position of the minor variant, the possible number of copies of the minor variant allele of the CYP2D6 gene, most likely, combinations in the minor variant indicate one or more minor variants of the CYP2D6 gene. In some embodiments, determining one or more minor variants of the CYP2D6 gene present in a subject includes: for each of a plurality of minor variant positions of the CYP2D6 gene, the minor variant position is associated with a minor variant allele of the CYP2D6 gene, determining the most likely combination of the possible number of copies of the minor allele the CYP2D6 gene at the minor variant position and the possible number of copies of the reference allele of the CYP2D6 gene at the minor variant position, summed with the copy number of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads with bases that contain the small variant allele of the CYP2D6 gene at the minor position variant, and (b) the number of sequence reads with bases that contain the reference allele of the CYP2D6 gene at the minor variant position, the possible copy number of small variant alleles of the CYP2D6 gene, most likely combinations at multiple minor variant positions indicating one or more minor variants of the CYP2D6 gene .

[0034] In some embodiments, the method includes: for a position of a minor variant of the CYP2D6 gene associated with a minor allele of a variant of the CYP2D6 gene, determining the most likely combination of the possible copy number of the small allele of the CYP2D6 gene at the position of the minor variant and the possible number of copies of the reference allele of the CYP2D6 gene at the position minor variant summed with the copy number of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads aligned to the CYP2D6 gene overlapping with the minor variant position and the base that contains the small variant allele of the CYP2D6 gene at the minor variant position and (b ) number of sequence reads aligned to the CYP2D6 gene overlapping the minor variant position and the base that contains the reference allele of the CYP2D6 gene at the minor variant position; and identifying one or more small variants of the CYP2D6 gene using the possible number of copies of the small variant allele of the CYP2D6 gene of the most likely determined combination. In some embodiments, the method includes: for each of a plurality of minor variant positions of the CYP2D6 gene, the minor 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 allele of the CYP2D6 gene at the position of the minor variant and the possible copy number of the reference allele of the CYP2D6 gene at the minor variant position, summed with the copy number of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads aligned to the CYP2D6 gene overlapping with the minor variant position, and the base that contains the CYP2D6 gene minor variant allele at the minor variant position, and (b) the number of sequence reads aligned to the CYP2D6 gene at the minor variant position, and the base that contains the reference allele of the CYP2D6 gene at the minor variant position; and identifying one or more small variants of the CYP2D6 gene using the possible copy number of small variant alleles of the CYP2D6 gene of the most likely combinations at a plurality of defined positions of the small variants.

[0035] In some embodiments, the position of the minor variant is in the homologous region of CYP2D6 / CYP2D7 , determining the most likely combination includes determining the most likely combination of the possible copy number of the small variant allele of the CYP2D6 gene at the position of the minor variant and the possible number of copies of the reference allele of the CYP2D6 gene at the position of the small variant summed with the copy number of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene with a base that contains the minor variant allele of the CYP2D6 gene at the minor variant position, and/or (b) the number sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, with a base that contains the CYP2D6 reference allele at the minor variant position. In some embodiments, the minor variant position is not in the CYP2D6 / CYP2D7 homologous region, determining the most likely combination includes determining 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, summed with the number of copies of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads aligned to the CYP2D6 gene, rather than the CYP2D7 gene, to the base that contains the small variant allele of the CYP2D6 gene at the minor variant position, and/or (b ) number of sequence reads aligned to the CYP2D6 gene rather than the CYP2D7 gene, with a base that contains the CYP2D6 reference allele at the minor variant position.

[0036] 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 CYP2D6 gene copy number at the minor variant position may include the CYP2D6 gene copy number of the possible CYP2D6 gene copy numbers of the most likely specific combinations. The CYP2D6 gene copy number at the minor variant position may include the CYP2D6 gene copy number of the possible CYP2D6 gene copy numbers of the most likely specific combinations and located closest to the minor variant position. The number of copies of the CYP2D6 gene at the minor variant position may include the number of copies of the CYP2D6 gene at the 5' position or at the 3' position of the small variant position. In some embodiments of the present invention, the method includes: (a) determining the number of sequence reads with bases that contain a minor allele variant of the CYP2D6 gene; and (b) determining the number of sequence reads with bases that contain the reference allele of the CYP2D6 gene.

[0037] In some embodiments, determining an allele of a CYP2D6 gene in a subject includes: determining alleles (eg, 2, 3, 4, 5 or more alleles) of the CYP2D6 gene present in the subject. In some embodiments, determining the allele of the CYP2D6 gene present in a subject includes: determining the star allele and/or haplotype of the CYP2D6 gene in the subject using one or more specific structural variants of the CYP2D6 gene, and/or one or more minor variants of the CYP2D6 gene, wherein the star allele is not necessarily associated with a known function.

[0038] In some embodiments, the method includes: determining the level of CYP2D6 enzymatic activity in a subject using a specific allele of the CYP2D6 gene. Enzyme activity can be low, medium, normal or ultra-fast. In some embodiments, the method includes determining a treatment dosage recommendation and/or treatment recommendation for a subject based on the CYP2D6 gene allele present in the subject.

[0039] Paralog genotyping is described herein. In some embodiments, the paralog genotyping system comprises: a non-volatile memory configured to store executable commands and sequencing data containing a plurality of sequence reads obtained from a sample from a subject aligned to a first or second paralog. The system may comprise: a processor (such as a hardware processor or virtual processor) communicating with non-volatile memory, the processor being programmed with executable instructions to perform: determining the number of copies of paralogs of the first type using a mixture of normal distributions model containing a plurality of distributions, each of which represents another integer copy number expressed by (i) the first number of sequence reads aligned to the first region. A hardware processor programmed with executable instructions to: determine the most likely combination for one of a plurality of paralog-specific bases from a plurality of possible combinations, each containing a possible number of copies of a first paralog of the first type and a possible number of copies of a second paralog of the first type, summed with the number of copies of certain paralogs of the first type, given (a) the number of sequence reads from the set of sequence reads with bases that contain a base specific to the first paralog, and (b) the number of sequence reads from the set of sequence reads with bases that contain a base specific to the second paralog corresponding to a base specific to the first paralogue. The hardware processor is programmed through executable instructions to: determine the copy number or allele of the first paralog using the most likely combination of the possible copy number of the first paralog and the possible copy number of the second paralog determined for the base specific to the first paralog. In some embodiments, the first paralog and the second paralog have at least 90% sequence identity.

[0040] In some embodiments, a hardware processor, through executable instructions, to: determine (i) a first number of sequence reads from a plurality of sequence reads in sequencing data obtained from a sample from a subject aligned to the first region. The method may include: determining (i) a first normalized number of sequence reads aligned to the first region using (i) the length of the first region, wherein determining the copy number of the first type of paralogs includes: determining the copy number of the first type of paralog using a mixture normal distribution model, given by (i) the first normalized number of sequence reads aligned to the first region. The hardware processor may be programmed using executable instructions to perform: may include: receiving sequencing data containing a plurality of sequence reads aligned to the first region.

[0041] In some embodiments, the hardware processor is programmed, through executable instructions, to: determine the copy number of one or more second type paralogs using a mixture of normal distributions defined by (ii) a second number of sequence reads aligned to the second region. Determining the copy number or allele of the first paralog may include: determining the copy number or allele of the first paralog using the most likely combination of the possible copy number of the first paralog and the possible copy number of the second paralog, determined for the base specific to the first paralog, and the copy number of one or more paralogs second type. The method may include: determining the number of copies of the third type of paralogs from the number of copies of the first type of paralogs and the number of copies of the second type of paralogs. Determining the copy number or allele of the first paralog includes: determining the copy number or allele of the first paralog using the most likely combination of the possible copy number of the first paralog and the possible copy number of the second paralog, determined for the base specific to the first paralog,

[0042] In some embodiments, the first paralogue is the survival motor neuron 1 ( SMN1 ) gene. The second paralog may be the survival motor neuron 2 ( SMN2 ) gene. The first region may comprise at least one exon 1 to 6 of the SMN1 gene and at least one exon 1 to 6 of the SMN2 gene. The second region may comprise at least exon 7 or 8 of the SMN1 gene and at least one of exon 7 or 8 of the SMN2 gene. The first type of paralogues may include an intact SMN1 gene and an intact SMN2 gene. One or more type two paralogs may include an intact SMN1 gene, an intact 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.

[0043] In some embodiments, the first paralog is a cytochrome P450 family 2 subfamily D member 6 gene ( CYP2D6 ). The second paralog may be the cytochrome P450 family 2 subfamily D member 7 gene ( CYP2D7 ). The first region may contain the CYP2D6 gene and the CYP2D7 gene. The second region may contain a spacer region between the CYP2D7 gene and the REP7 repeat element downstream of the CYP2D7 gene. Paralogues of the first type may contain the CYP2D6 gene and the CYP2D7 gene. One or more type two paralogues may contain a CYP2D6 / CYP2D7 fusion allele with a spacer region and a REP7 repeat element downstream of the CYP2D6/CYP2D7 fusion allele. The copy number of the first paralog may contain an allele of the CYP2D6 gene present in the subject that is a minor variant or structural variant of the CYP2D6 gene.

[0044] Described herein are embodiments of a system (eg, a computer system) comprising non-volatile memory configured to store executable instructions; and a processor (eg, a hardware processor or a virtual processor) in connection with the non-volatile memory, the hardware processor being programmed by executable instructions to perform any method described herein. Described herein are embodiments of a device (eg, an electronic device) comprising non-volatile memory configured to store executable instructions; and a processor (eg, a hardware processor or a virtual processor) in connection with the non-volatile memory, the hardware processor being programmed by executable instructions to perform any method described herein. Embodiments described herein are computer-readable media comprising executable instructions that, when executed by a processor (e.g., a hardware processor or virtual processor) of a system or device, cause the hardware processor to execute any method described herein.

[0045] A detailed description of one or more embodiments is provided in the accompanying drawings and description below. Other features, aspects and advantages will become apparent from the description, drawings and claims. Neither this summary nor the following detailed description is intended to define or limit the scope of the invention.

Brief description of graphic materials

[0046] In FIG. 1A-1E provide explanations for determining the copy number of SMA and SMN in accordance with one embodiment of the method described herein.

[0047] In FIG. 2A-2C show SMN1/2 copy number population distributions determined using one embodiment of the method described herein.

[0048] In FIG. Figure 3 shows the SMA identified in two triplets in the Next Generation Children project and validated using MLPA.

[0049] In FIG. 4 shows population frequencies determined using one embodiment of the method described herein, consistent with previous studies.

[0050] In FIG. 5 shows a non-limiting example of an IGV snapshot showing that CYP2D6 is highly polymorphic and is located downstream of CYP2D7 , a pseudogenic paralog of CYP2D6 .

[0051] In FIG. 6 shows a non-limiting example of a schematic illustration of deletions, duplications and fusions of the CYP2D6 /7 genes.

[0052] In FIG. 7 is a non-limiting example of a plot showing that allelic frequencies determined by the method are consistent with the PharmVar Database Pharmacogene Variation (PharmVar) Consortium.

[0053] In FIG. 8 is a flow diagram showing an example of a method for determining the copy number of the survival motor neuron 1 ( SMN1 ) gene using sequencing data.

[0054] In FIG. 9 is a flow diagram showing an example of a method for genotyping the cytochrome P450 family 2 subfamily D member 6 ( CYP2D6 ) gene using sequencing data.

[0055] In FIG. 10 is a flow diagram showing an example of a method for genotyping a paralog using sequencing data.

[0056] In FIG. 11 is a block diagram of an exemplary computing system configured to implement paralog genotyping using sequencing data.

[0057] In FIG. 12A and 12B show non-limiting example graphs illustrating common ICCs affecting the SMN1 / SMN2 loci. In FIG. 12A shows depth profiles in the SMN1/SMN2 regions. Samples with total SMN1 + SMN2 copy numbers of 2, 3, 4, and 5 are shown as dots, respectively. For each copy number category, the depth of 50 samples is summed. Each point represents normalized depth values in a 100-bp window. Read counts were calculated in each 100-bp window, the reads for SMN1 and SMN2 were summed, and normalized to the depth of the wild-type samples (CN=4). SMN exons are represented as magenta boxes. The two x-axes show the coordinates in SMN1 (bottom) and SMN2 (top). In FIG. 12B shows depth profiles pooled from 50 samples carrying deletions of exons 7 and 8, shown as dots. Read depth values were calculated in the same manner as shown in FIG. 12A.

[0058] In FIG. 13 shows a non-limiting example of a scatterplot of total SMN copy number ( SMN1 + SMN2 ) (x-axis, denotes read depth of exons 1-6) and intact SMN copy number (y-axis, denotes read depth of exons 7-8).

[0059] In FIG. 14A-14D show the copy number distributions of SMN1 / SMN2 / SMN * in the population. In FIG. 14A is a non-limiting example of an illustrative graph illustrating the percentage of samples showing copy number assignment agreement with c.840C>T across 16 SMN1-SMN2 base difference sites in African and non-African populations. Site 13* is a splice variant site c.840C>T. The black horizontal line indicates 85% agreement. In FIG. 14B shows non-limiting examples of histograms of copy number distributions of SMN1 , SMN2 and SMN* across five populations in the 1kGP and NIHR BioResource cohort (numbers given in Table 15). In FIG. 14C shows a non-limiting example of a plot of SMN1 copy number versus total SMN2 copy number (intact SMN2 + SMN* ). In FIG. Figure 14D shows two trios with an SMA proband identified by a specialist and orthogonally confirmed in the NIHR BioResource cohort. The copy numbers per allele of SMN1 , SMN2 , and SMN* are phased and labeled for each member of the trio.

[0060] In FIG. Figure 15 shows non-limiting example plots, each illustrating the posterior probability distribution for modeling SMN1 copy number using a single site at different read depths and SMN1:SMN2 copy number combinations.

[0061] In FIG. 16 shows a non-limiting example of an IGV snapshot of the SMN2 region in a sample with a deletion in exon 7-8. Horizontal lines pair two reads on the central alignment track. BLAT results for two split reads spanning the breakpoint are shown in the bottom track, showing two segments of the same alignment for the read on either side of the deletion breakpoint.

[0062] In FIG. 17 shows non-limiting examples of graphs illustrating the correlation between raw SMN1 copy numbers at 15 base differences near c840.C>T and raw SMN1 copy numbers at the c840.C>T site. The raw copy number of SMN1 at each site was calculated as the copy number of intact SMN multiplied by the fraction of SMN1 supporting the number of SMN1 + SMN2 reads from the number of reads. Correlation coefficients are given in the header of each graph.

[0063] In FIG. 18A and 18B show non-limiting examples of SMN1 / SMN2 haplotype plots in samples with SMN1 :2 SMN2 :0 and SMN1 :2 SMN2 :1 in 1kGP. The y-axis shows the raw copy numbers of SMN1 as shown in FIG. 16. The X-axis shows the 16 sites, the indexes of which are listed and explained in Table 8. Index #13 represents the site c840.C>T. Samples with SMN1: 2 SMN2 :0 are shown together in the top left graph. Samples with SMN1 :2 SMN2 :1 are shown as 5 clusters. FIG. 18A Non-African Population FIG. 18B African population

[0064] In FIG. 19 shows a non-limiting example of an IGV snapshot showing a 1.9 kb deletion. in SMN1 in MB509.

[0065] In FIG. 20 is a non-limiting illustrative graph illustrating SMN1/SMN2 /SMN* copy numbers in the 1kGP and NIHR cohorts.

[0066] In FIG. 21A and 21B show discrepancies and missing determinations in the test samples.

[0067] In FIG. 22 shows copy number determinations obtained from BWA and Isaac BAM.

[0068] In FIG. 23 is a non-limiting example plot showing the quality of WGS data in the CYP2D6 /7 region. The average mapping quality for the 1kGP samples is plotted for each position in the CYP2D6 /7 region. In a window 200 bp long. a median filter is used. REP6, REP7, and CYP2D6 /7 exons 9 are shown as boxes on the left ( CYP2D6 ) and right ( CYP2D7 ). Two 2.8 kb repeat regions. downstream, CYP2D6 (REP6) and CYP2D7 (REP7) are identical and essentially incompatible. The dotted frame indicates the spacer region between CYP2D7 and REP7. The two major homologous regions in the genes are shaded.

[0069] In FIG. Figure 24 shows structural variants confirmed by PacBio CCS reads. PacBio reads confirm deletion (*5), duplication, and fusion (*36, *68, and *13). Graphs were generated using sv-viz2 (zotero.org/google-docs/?xAunA6). For deletions and duplications, due to significant homology in the REP regions, the exact location of the breakpoints in the REP is not available. The break points in A and B are for illustration purposes only.

[0070] In FIG. 25 shows a non-limiting example plot showing CYP2D6 allele frequencies for five ethnic populations for the ten most common CYP2D6 altered function haplotypes. One haplotype (*2x2) is characterized by increased function, two haplotypes (*4 and *4 + *68) by lack of function, and the remaining haplotypes by decreased function.

[0071] In FIG. 26 shows that sites of CYP2D6 / CYP2D7 base differences are highly variable in the population. The y-axis shows the frequency of samples in which the CN for the CYP2D6 base is identified in 2 of all samples having a total CYP2D6 + CYP2D7 CN of 4. The x-axis shows the genome coordinates in hg38. CYP2D6 exons are shown as gray boxes above the graph. The black horizontal line indicates the 98% cutoff.

[0072] In FIG. Figure 27 shows unprocessed CYP2D6 CNs at CYP2D6/7 differentiation sites in SV examples. Raw CYP2D6 CNs were calculated as the total CYP2D6+CYP2D7 CNs multiplied by the ratio of CYP2D6 supporting reads from CYP2D6 to CYP2D7 supporting reads. The large diamond indicates the number of copies of genes derived from CYP2D6 at the end of the gene (can be a complete CYP2D6 gene or a fusion gene ending with CYP2D6 ), calculated as the total CN of CYP2D6+CYP2 D7 minus the CN of the CYP2D7 spacer region (see FIG. 23). To detect SV, the CN of CYP2D6 was determined at each site, and the change in CN of CYP2D6 in the gene indicated the presence of SV. For example, in HG01161, CN CYP2D6 changed from 2 to 1 between exon 7 and exon 9, indicating a CYP2D7-CYP2D6 hybrid gene. In HG00553, CN CYP2D6 changed from 2 to 3 between exon 1 and exon 2, indicating a CYP2D6-CYP2D7 hybrid gene.

[0073] In FIG. 28 shows that PacBio data support the *10D fusion in HG00421. For comparison, a sample with *36 (HG00612) is shown. PacBio reads containing fusions are reads with shaded bases, bases that represent a software cross-link made by the aligner and were derived from part of the CYP2D7 fusion. The fusion breakpoints are close to each other, but the breakpoint for *36 is located upstream of the base differences in exon 9 (located within the black block), and the breakpoint for *10D is located downstream, leaving the CYP2D6 gene intact.

[0074] In FIG. 29 shows that the PacBio data had a false *61 ( CYP2D6 / CYP2D7 hybrid) produced 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). The PacBio data showed no structural variant in this region (each read was fully aligned, with no soft stitches indicating unaligned parts).

[0075] In FIG. 30A and 30B show the new haplotype *10+*36+*36+*83 in HG00597. FIG. 30A The depth chart shown in FIG. 27 shows that HG00597 had three copies of *36-like fusions, all of which had a breakpoint in the homologous region between exon 7 and exon 9. FIG. 30B screenshot of IGV data from PacBio showing all reads containing fusions, i.e. those aligned by soft stitching. One copy of the fusion gene lacked g.42130692G>A, a SNP that was located at *36 but not at *83, as shown in the region flanked by the two black vertical lines. This copy was *83, and unlike what was reported in PharmVar, it was a fusion gene with REP7 rather than REP6, otherwise the copy number of the region downstream of exon 9 would have been 3 instead of 2 in FIG. 30A.

[0076] In FIG. 31A and 31B compared the frequencies of 1kGP and pharmGKB. Each dot represents a haplotype with a frequency that is greater than or equal to 0.5% in either 1kGP or pharmGKB. SV-associated haplotypes were noted, including the two haplotypes with the greatest deviation (*10+*36 in East Asians and *4+*68 in Europeans). Other haplotypes with deviant values are labeled (*2, *41, *34, *39, *2, and *29). A diagonal line is drawn for each panel. Correlation coefficients are given for each population (*10+*36 excluded for East Asians and *4+*68 excluded for Europeans for calculation). In FIG. 31B shows values in the lower range of values (<5%).

[0077] In FIG. 32 shows a non-limiting example of an IGV snapshot showing de novo assembly of PacBio reads into HG00733 not including the *68 fusion.

[0078] Throughout the drawings, reference numbers may be reused to indicate correspondence between reference elements. The drawings are presented to illustrate the embodiments described herein and are not intended to limit the scope of the description.

Detailed description of the invention

[0079] The following detailed description contains references to related drawings that form a part of this document. In graphic materials, similar symbols generally indicate similar components unless the context otherwise requires. The illustrative embodiments described in the detailed description, drawings and claims are not intended to be limiting. Other embodiments and other changes may be used without departing from the spirit or scope of the claimed subject matter presented herein. It should be understood that aspects of this specification generally set forth herein and illustrated in the figures may be rearranged, substituted, combined, separated, and constructed in a wide variety of different configurations, all of which are expressly contemplated herein and are part of the specification in this document.

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

[0081] Described herein are methods for determining the copy number of the survival motor neuron 1 ( SMN1 ) gene and/or 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 controlled by a processor (such as a hardware processor or virtual processor) and includes: obtaining sequencing data comprising a plurality of sequence reads from a sample from a subject aligned to the SMN1 gene or the SMN2 gene. The method may include: determining (i) a first number of sequence reads from a plurality of sequence reads aligned to a first region of the SMN1 or SMN2 gene containing at least one of the first 6 exons of the SMN1 gene or the SMN2 gene, respectively, and (ii) a second number of sequence reads from a plurality of sequence reads aligned to a second region of SMN1 or SMN2 containing at least exon 7 or exon 8 of the SMN1 gene or the SMN2 gene, respectively. The method may include: determining (i) a first normalized read count of a sequence aligned to a first region of the SMN1 or SMN2 gene, and (ii) a second normalized read count of a sequence aligned to a second SMN1 or SMN2 region, using (i) the length of the first SMN1 region or SMN2 and (ii) the length of the second region SMN1 or SMN2 , respectively. The method may include: determining (i) the copy number of full-length survival motor neuron (SMN) genes, each of which is an intact SMN1 gene, an intact SMN2 gene, a truncated SMN1 gene, or a truncated SMN2 gene; and (ii) the copy number of any intact SMN genes, each of which is an intact SMN1 gene or an intact SMN2 gene, using a mixture normal distribution model containing a plurality of normal distributions, each representing a different integer copy number, subject to (i) the first a normalized read count of the sequence aligned to the first SMN1 or SMN2 region, and (ii) a second normalized read count of the sequence aligned to the second SMN1 or SMN2 region, respectively. The method may include: for one set of SMN1 gene-specific bases associated with an intact SMN1 gene, determining the most likely combination from a plurality of possible combinations, each of which contains a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene, summed with the number of copies of any determined intact SMN genes, considering (a) the number of sequence reads from a set of sequence reads with bases that contain a base specific to the SMN1 gene, and (b) the number of sequence reads from a set of sequence reads with bases that contain a base specific to the SMN2 gene , corresponding to a base specific to the SMN1 gene. The method may include: determining the copy number of the SMN1 gene and/or the SMN2 gene using the most likely combination of the possible SMN1 gene copy number and the possible SMN2 gene copy number determined for the SMN1 gene-specific base.

[0082] Described herein are methods for genotyping the cytochrome P450 family 2 subfamily D member 6 gene ( CYP2D6 ). In some embodiments, a method for genotyping a CYP2D6 gene is controlled by a processor (such as a hardware processor or virtual processor) and includes: obtaining sequencing data comprising multiple sequence reads from a sample from a subject aligned to a CYP2D6 gene or a cytochrome P450 family 2 subfamily D member 7 gene ( CYP2D7 ). The method may include: determining (i) a first number of sequence reads from a plurality of sequence reads aligned to a CYP2D6 gene or a CYP2D7 gene. The method may include: determining (i) the first normalized read count of the sequence aligned to the CYP2D6 gene or the CYP2D7 gene, using (i) the length of the CYP2D6 gene or the CYP2D7 gene, respectively. The method may include: determining (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene using a mixture of normal distributions model containing a plurality of normal distributions, each representing a different integer copy number, based on (i) the first normalized read count of the sequence aligned to CYP2D6 genome or CYP2D7 genome. The method may include: for one of a plurality of bases specific to the CYP2D6 gene, determining the most likely combination of a plurality of possible combinations, each of which contains a possible number of copies of the CYP2D6 gene and a possible number of copies of the CYP2D7 gene, summed to the total number of copies of the CYP2D6 gene and the particular CYP2D7 gene, considering (a) multiple sequence reads with bases that contain a CYP2D6 gene-specific base, and (b) multiple sequence reads with bases that contain a CYP2D7 gene-specific base corresponding to a CYP2D6 gene-specific base. The method may include: determining the allele of the CYP2D6 gene present in the subject using the most likely combination of the possible number of copies of the CYP2D6 gene and the possible number of copies of the CYP2D7 gene determined for the base specific for the CYP2D6 gene.

[0083] Methods for genotyping paralogs are described herein. In some embodiments, a method for genotyping paralogs is controlled by a processor (such as a hardware processor or virtual processor) and includes: obtaining sequencing data comprising a plurality of sequence reads from a sample from a subject aligned to a first paralog or a second paralog. The method may include: determining the copy number of the first type of paralogs using a mixture normal distribution model comprising a plurality of distributions, each of which represents a different integer copy number given by (i) a first number of sequence reads aligned to the first region. The method may include: determining the most likely combination for one of a plurality of bases specific to the first paralog from a plurality of possible combinations, each of which contains a possible number of copies of the first paralog of the first type and a possible number of copies of the second paralog of the first type, summed with the number of copies of the specified paralogs of the first type, given (a) the number of sequence reads from the set of sequence reads with bases that contain a base specific to the first paralog, and (b) the number of sequence reads from the set of sequence reads with bases that contain a base specific to the second paralog corresponding to the base specific to the first paralog. The method may include: determining the copy number or allele of the first paralog using the most likely combination of the possible copy number of the first paralog and the possible copy number of the second paralog determined for a base specific to the first paralog.

[0084] Described herein are embodiments of a system (eg, a computer system) comprising non-volatile memory configured to store executable instructions; and a processor (eg, a hardware processor or a virtual processor) in connection with the non-volatile memory, the hardware processor being programmed by executable instructions to perform any method described herein. Described herein are embodiments of a device (eg, an electronic device) comprising non-volatile memory configured to store executable instructions; and a processor (eg, a hardware processor or a virtual processor) in connection with the non-volatile memory, the hardware processor being programmed by executable instructions to perform any method described herein. Embodiments described herein are computer-readable media comprising executable instructions that, when executed by a processor (e.g., a hardware processor or virtual processor) of a system or device, cause the hardware processor to execute any method described herein.

Diagnosis of spinal muscular atrophy and carrier testing based on whole genome sequencing data

[0085] Spinal muscular atrophy (SMA) is characterized by weakening of skeletal muscles and is a leading genetic cause of death in early childhood with an incidence of 1 in 6,000-10,000 live births and a carrier rate of 1:40-8,012. SMA is caused by mutations in the SMN1 (survival motor neuron 1) gene (FIG. 1A). The duplicated SMN2 gene differs from SMN1 by only a few base pairs, one of which, the c.840C>T splice variant in exon 7, has functional consequences. By interrupting the splicing enhancer, the c.840C>T mutation results in increased exon 7 skipping and decreased full-length transcripts in SMN2 3 (FIG. 1B-1D). The genomic region is subject to unequal crossing over and gene conversion, resulting in variable copy numbers of SMN1 and SMN2 (FIG. 1B). Due to the high incidence and severity of the disease, extensive screening for SMA is recommended, and a key factor for this screening is determination of SMN1 copy number for SMA diagnosis and carrier testing. In addition, SMN2 copy number determines the severity of SMA and is important for clinical classification and prognosis.

[0086] Standard SMA carriage tests use PCR-based methods such as multiplex ligation probe amplification (MLPA), quantitative PCR (qPCR), and digital PCR. These methods mainly target the c.840C>T site. Incorporating SMA screening into high-throughput NGS-based assays that can simultaneously profile large numbers of genes or even the entire genome may be beneficial. The near-perfect sequence identity between SMN1 and SMN2 makes variant detection difficult by standard NGS-based methods.

[0087] Described herein is an SMN copy number detector based on a bioinformatics method that determines the copy number of SMN1 and SMN2 with whole genome sequencing (WGS) data (FIG. 1E). The method may include determining the copy number of SMN1 + SMN2 in two regions, exons 1-6 and exons 7-8, by summing the reads in SMN1 and SMN2 . The method may include differentiating SMN1 from SMN2 using read counts at fixed base differences. In some embodiments, the method does not include realigning the aligned sequences to a modified reference pattern. The method represents the first SMN copy number tool that can identify both SMA patients and carriers from WGS data. Some embodiments of the method are not limited to exons 7 and 8 and are not primarily focused on c.840C>T. This method uses a gene-wide approach and provides the most comprehensive set of determinations, including copy number of full-length SMN1 and SMN2 , as well as the truncated form of SMN with deletion of exons 7 and 8. This method can be easily applied to any WGS data and will be a useful tool for SMA diagnosis and screening carriers for inclusion in high-throughput population-wide WGS screening.

[0088] In FIG. 1A-1E show explanations of the SMA copy number determination of SMN in accordance with one embodiment of the bioinformatics method described herein. Table 1 shows the differentiation of SMN1 from SMN2 based on a fixed single nucleotide polymorphism (SNP) in accordance with an embodiment of the method. SMN1 copy number detection occurs at 16 sites near c.840C>T. Nine sites with high percentage matches to c.840C>T selected to jointly determine SMN1 copy number. In FIG. 2A-2C and Table 2 show the population distribution of specific SMN1/2 copy numbers. More copies of SMN1 were observed with fewer copies of SMN2 , suggesting gene conversion as a mechanism for SMN1 and SMN2 copy number variation. In table 3 shows the verification of the copy number determined using the bioinformatics method in comparison with the copy number determined by the digital PCR method. Validation against digital PCR showed 100% copy number agreement in SMN1 and 98% in SMN2 . In FIG. Figure 3 shows the SMA identified in two triplets in the Next Generation Children project and validated using MLPA. In FIG. 4 and Table 4 show the frequencies in the populations determined using the bioinformatics method, consistent with previous studies.

Table 1. Differentiation of SMN1 from SMN2 based on fixed single nucleotide polymorphism (SNP)

Site no. Location Selected SMN1 SMN2 Percentage of compliance with c.840 C>T Position, hg19 Grounds Position, hg19 Grounds 1 Intron 6 70244142 A 69368717 G 85.2 2 Intron 6 70245876 T 69370451 C 85.5 3 Intron 6 70246016 G 69370591 A 94.8 4 Intron 6 70246019 T 69370594 C 94.1 5 Intron 6 70246156 G 69370731 A 93.9 6 Intron 6 70246167 T 69370742 C 59.8 7 Intron 6 70246320 G 69370895 A 96.3 8 Intron 6 Yes 70246793 G 69371368 A 99 9 Intron 6 Yes 70246919 A 69371499 C 98.2 10 Intron 6 Yes 70247219 G 69371799 A 98.8 eleven Intron 6 Yes 70247290 T 69371870 C 99 12 Intron 6 Yes 70247724 G 69372304 A 99.6 13 Exon 7
(c.840 C>T) Yes 70247773 C 69372353 T 100 14 Intron 7 Yes 70247921 A 69372501 G 99.5 15 Intron 7 Yes 70248036 A 69372616 G 99.6 16 Exon 8 Yes 70248501 G 69373081 A 97.9

Table 2. Distribution of SMN1/2 copy numbers in populations

Number of subjects Full size SMN1 Full size SMN2 N < 2235 9 1 1 24 1 2 10 1 3 2 1 4 162 2 0 827 2 1 1005 2 2 45 2 3 1 2 4 15 3 0 75 3 1 41 3 2 7 3 3 1 3 4 5 4 0 3 4 1 3 4 2

Table 3. Verification of the copy number determined using the bioinformatics method in comparison with the copy number determined by the digital PCR method

Sample ID SMA status SMN Copy Number Determinant Digital PCR Agreed Full Length Full Length Number of SMN copies* SMN1 copy number SMN2 copy number SMN1 copy number SMN2 copy number NA03813 Affected 0 3 0 0 3 Yes NA09677 Affected 0 3 0 0 3 Yes NA23689 Affected 0 3 0 0 3 Yes NA00232 Affected 0 2 0 0 2 Yes NA10684 Affected 0 2 0 0 2 Yes NA23687 Carrier 1 2 0 1 2 Yes NA23688 Carrier 1 2 0 1 2 Yes NA03815 Carrier 1 1 0 1 1 Yes

Table 4. Population frequencies determined using bioinformatics methods, consistent with previous studies

Carrier SMN1 SMN1 SMN1 Frequency (%) CN=1 CN=2 CN>=3 Representatives of the Caucasian race This study 2.2 92.2 5.6 Hendrickson et al. ^a 2.7 91 6.3 Sugarman et al. ^b 2.02 90.9 7.05 Representatives of the Negroid race This study 0.44 44.8 54.8 Hendrickson et al. ^a * 1.1 52.1 46.8 Sugarman et al. ^b * 0.98 51.9 47.1

aHendrickson et al. Differences in SMN1 allele frequencies among ethnic groups within North America. J Med Genet. 2009;46(9):641-644. doi:10.1136/jmg.2009.066969.

bSugarman et al. Pan-ethnic carrier screening and prenatal diagnosis for spinal muscular atrophy: clinical laboratory analysis of >72,400 specimens. Eur J Hum Genet. 2012; 20(1): 27-32.doi: 10.1038/Ejhg.2011.134.

*African Americans

Characterization of medically applicable variants from 2500 publicly available high-depth genomes of various origins

[0089] Population-based whole genome sequencing (WGS) data is becoming increasingly available. For example, publicly available sequencing data is available, such as high depth (>30X) WGS data for >2500 samples from the 1000 Genomes Project (1kGP). This has significantly improved the clinical interpretation of simple single nucleotide variants (SNVs) and insertions/deletions (indices). However, many medically important regions and variants, such as triplet repeats and homologs, are not included in WGS-based databases because annotation of these regions and variants requires specialized bioinformatics techniques. Characterization of known clinical variants at the population level is therefore necessary to maximize the impact of population sequencing experiments. In some embodiments, the methods described herein address three shortcomings of standard secondary analysis processes: 1) spinal muscular atrophy (SMA) detection and carrier screening, 2) genotyping of CYP2D6 for pharmacogenomic applications, and 3) detection of repeat triplet expansions. The methods can be used to determine SMN1 /2 copy number, CYP2D6 star alleles, repeat expansions in a 1kGP population, and quantify differences between subpopulations. This paper describes subpopulation frequency distributions and perpendicular validation of these methods using supporting data generated from high-quality long reads.

CYP2D6

[0090] CYP2D6 is an important drug metabolizing enzyme that is highly polymorphic (FIG. 5). CYP2D6 has high sequence similarity to its pseudogenic paralogue ( CYP2D7 ). Genotyping CYP2D6 by WGS is challenging due to common gene conversions between CYP2D6 and CYP2D7 (hereafter referred to as CYP2D6 /7), common SVs ( CYP2D6 /7 gene deletions, duplications, and fusion genes; See FIG. 6 for illustration), and also the sequence similarity of CYP2D/7 , which leads to ambiguous read alignment for any of the genes (FIG. 5). This document describes a CYP2D6 detector based on a bioinformatics method that is capable of identifying (eg, pinpointing) diplotypes targeting star alleles (eg, all star alleles) with known functions. In some embodiments, the method includes the following steps:

[0091] 1. Determination of the total copy number of CYP2D6 + CYP2D7 .

[0092] 2. Determination of CNVs/hybrids based on copy number determinations at CYP2D6 / CYP2D7 differentiation sites.

[0093] 3. Determine 56 SNPs/indexes from BAM (or other file containing sequence reads).

- Use of information about the number of copies.

- reading results in both positions of CYP2D6 and CYP2D7 in homologous regions.

[0094] 4. Determination of star alleles and diplotypes based on all identified variants.

[0095] Table 5 shows the results of testing the CYP2D6 star allele determinations made by this method. The CYP2D6 star allele definitions obtained by the method for 92 of 96 samples are consistent with the GeT-RM consensus definitions from multiple platforms. This method outperformed identifiers such as Aldy ( CYP2D6 star allele determination for 89 of 96 samples agreed with the GeT-RM consensus) and Stargazer ( CYP2D6 star allele determination for 83 of 96 samples agreed with the GeT-RM consensus).

Table 5. Validation of the CYP2D6 determinant.

Sample Determination of CYP2D6 GeT-RM Consensus Aldy Stargazer NA24008 *1/*4+*68 *1/*4 *1/*4+*68 *1/*4+*68 NA21781 *2x2/*4+*68 *2x2/*68+*2 *2x2/*4+*68 *2x2/*4+*68 NA23874 *4/*4+*68 *4/*4 Undefined *4/*4+*68 NA18565 *10/*10+*36 *10/*36x2 *10/*10+*36 *10/*10+*36

[0096] In FIG. 7 shows that the allelic frequencies determined by the method are consistent with the PharmVar Database Pharmacogene Variation (PharmVar) Consortium.

Determining motor neuron survival 1 gene copy number using sequencing data

[0097] In FIG. 8 is a flow diagram showing an example of a method 800 for determining the survival copy number of a motor neuron 1 gene using sequencing data, such as whole genome sequencing data. The method 800 may be implemented as a set of executable program instructions stored on a computer-readable medium, such as one or more computer system disk drives. For example, the computing system 1100 shown in FIG. 11 and described in more detail below, may execute a set of executable program instructions to implement method 800. When method 800 is initiated, the executable program instructions may be loaded into a storage device, such as RAM, and executed by one or more processors of the computing system 1100. Although method 800 is described with respect to the computing system 1100 shown in FIG. 11, the description is provided for illustrative purposes only and is not intended to be limiting. In some embodiments, method 800 or portions thereof may be performed sequentially or in parallel by multiple computing systems.

[0098] After method 800 begins at step 804, method 800 proceeds to step 808, where a computing system (such as computing system 1100 described with reference to FIG. 11) determines (i) the first number of reads of the sequence from the plurality of reads sequences aligned to the first region of the survival motor neuron 1 ( SMN1 ) or survival motor neuron 2 ( SMN2 ) gene containing at least one of exons 1-6 of the SMN1 gene or the SMN2 gene, respectively, and (ii) a second number of sequence reads from a plurality of sequence reads aligned to a second region of the SMN1 or SMN2 gene containing at least exon 7 or exon 8 of the SMN1 gene or the SMN2 gene, respectively. The first number of sequence reads aligned to the first SMN1 or SMN2 region (or the second number of sequence reads aligned to the second SMN1 or SMN2 region) may be or approximately equal to, 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 more.

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

[0100] In some embodiments, the computing system receives sequencing data containing multiple sequence reads obtained from a sample from a subject aligned to the SMN1 gene or the SMN2 gene. Sequencing data may include whole genome sequencing (WGS) data or short-read WGS data. In some embodiments, the subject is a fetal subject, a neonatal subject, a pediatric subject, an adolescent subject, or an adult subject. The sample may contain cells or extracellular DNA. The sample may contain fetal cells or extracellular fetal DNA.

[0101] In some embodiments, a sequence read from a plurality of sequence reads is aligned to a first region SMN1 or SMN2 or a second region SMN1 or SMN2 with an alignment quality score of approximately zero. The alignment quality may be, for example, 0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10 or more (on a scale of 0 to 1 from the alignment score).

[0102] Method 800 proceeds from step 808 to step 812, where the computing system determines (i) a first normalized read count of the sequence aligned to the first SMN1 or SMN2 region, and (ii) a second normalized read count of the sequence aligned to the second SMN1 or SMN2 region. SMN2 , using (i) the length of the first region SMN1 or SMN2 and (ii) the length of the second region SMN1 or SMN2 , respectively. 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) may be or approximately equal to, for example, 1, 2, 3, 4, 5, 6, 7, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or more. The length of the first region of SMN1 or SMN2 may be or is approximately equal to, for example, 3 kb, 6 kb, 9 kb, 12 kb, 15 kb n., 18 kb., 21 kb., 22.2 kb., 24 kb. or more. The length of the second region of SMN1 or SMN2 may be or approximately equal to, for example, 3 kb, 6 kb. or more.

[0103] 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, the computing system may determine (i) a first normalized the 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 (i) the length of the first SMN1 or SMN2 region and (ii) the length of the second SMN1 region or SMN2 , respectively, and (iii) the depth of sequence reads of the subject's genomic region, excluding genetic loci containing the SMN1 gene and the SMN2 gene, as determined by sequencing data. The sequence read depth of the subject's genome region, excluding genetic loci containing the SMN1 gene and the SMN2 gene, as determined by sequencing data may be or approximately equal to, for example, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or more.

[0104] To determine (i) a first normalized number of sequence reads aligned to the first region SMN1 or SMN2 , and (ii) a second normalized number of sequence reads aligned to the second region SMN1 or SMN2 , the computing system determines (i) the normalized number of sequence reads along the length of the first SMN1 or SMN2 region aligned to the first SMN1 or SMN2 region, and (ii) the normalized sequence read count along the length of the second SMN1 or SMN2 region aligned to the second SMN1 or SMN2 region, using (i) the length of the first SMN1 region or SMN2 and (ii) the length of the second region SMN1 or SMN2 , respectively. The computing system may determine (i) a first normalized depth of sequence reads aligned to a first SMN1 or SMN2 region, and (ii) a second normalized depth of sequence reads aligned to a second SMN1 or SMN2 region, based on (i) a normalized number of sequence reads by length the first SMN1 or SMN2 region and (ii) the normalized number of sequence reads along the length of the second SMN1 or SMN2 region, respectively, using the sequence read depth of the subject's genomic region excluding genetic loci containing the SMN1 gene and the SMN2 gene. 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 may represent the first normalized depth and the second normalized depth, respectively.

[0105] 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, the computing system may determine (i) a first normalized the 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 (i) the GC content of the first SMN1 or SMN2 region and (ii) the GC content of the second SMN1 or SMN2 region, respectively, and (iii) the sequence read depth of the subject's genomic region, excluding genetic loci containing the SMN1 gene and the SMN2 gene, as determined by sequencing, and (iv) the GC content of the genomic region. The GC content of the first region of SMN1 or SMN2 (or the GC content of the second region of SMN1 or SMN2 ) may be at or about equal to, 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 sequence read depth of a region of a subject's genome, excluding genetic loci containing the SMN1 gene and the SMN2 gene, as determined by sequencing may be or approximately equal to, for example, 3, 4, 5, 10, 20, 30, 40, 50, 100, or more. The GC content of a region of a subject's genome, excluding genetic loci containing the SMN1 gene and the SMN2 gene, as determined by sequencing may be or approximately equal to, for example, 40%, 41%, 42%, 43%, 44%, 45%, 46% , 47%, 48%, 49%, 50%, 51%, 52%, 53%, 54%, 55%, 56%, 57%, 58%, 59% or 60%.

[0106] In some embodiments, the region depth includes the average sequence read depth of a region of the subject's genome, excluding genetic loci containing the SMN1 gene and the SMN2 gene, as determined by sequencing data. The region depth may include the median sequence read depth of a region of the subject's genome, excluding genetic loci containing the SMN1 gene and the SMN2 gene, as determined by sequencing. The depth of the region may be or approximately equal to, for example, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or more. The region length may comprise approximately 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000 or more preselected regions of approximately 0.5 kb, 1 kb, 1.5 k in length. bp, 2 kb, 2.5 kb or 3 kb. in each throughout the subject's genome. For example, a region may contain approximately 3000 preselected regions of approximately 2 kb in length. in each throughout the subject's genome.

[0107] 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 more. For example, (i) the first normalized read count of the sequence aligned to the first SMN1 or SMN2 region and/or (ii) the second normalized read count of the sequence aligned to the second SMN1 or SMN2 region is from about 30 to about 40.

[0108] Method 800 proceeds from step 812 to step 816, where the computing system determines (i) the copy number of common motor neuron survival (SMN) genes and (ii) the copy number of any intact gene(s) ) SMN using a mixture normal distribution model containing a plurality of distributions, each representing a different integer copy number given by (i) a first normalized number of sequence reads aligned to the first SMN1 or SMN2 region, and (ii) a second normalized number of sequence reads, aligned with the second region of SMN1 or SMN2 , respectively. Common motor neuron survival genes may comprise an intact SMN1 gene, an intact SMN2 gene, a truncated SMN1 gene, and/or a truncated SMN2 gene. Any intact SMN gene(s) may contain an intact SMN1 gene and/or an intact SMN2 gene. The copy number of the total SMN gene(s) (or any gene(s) herein) may be or be approximately, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more. The copy number of any intact SMN gene(s) (or any gene(s) herein) may be or be approximately, for example, 1, 2, 3, 4, 5 , 6, 7, 8, 9, 10 or more.

[0109] In some embodiments, the mixture normal distribution model is a one-dimensional mixture normal distribution model. The set of normal distribution model distributions can represent integer copy numbers, for example, 0 to 5, 0 to 6, 0 to 7, 0 to 8, 0 to 9, 0 to 10, 0 to 11, 0 to 12, 0 to 13, 0 to 14, or 0 to 15. For example, a set of mixture normal distribution model distributions can represent integer copy numbers from 0 to 10. The mean (for example, 1, 2, 3, 4, 5 , 6, 7, 8, 9, 10, or more) of each of the plurality of distributions may be an integer number of copies represented by the distribution (e.g., copy number 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more). The standard deviation of the distributions may be or approximately be, for example, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 or more .

[0110] In some embodiments, to determine (i) the copy number of common SMN gene(s) and (ii) the copy number of any intact SMN gene(s), the computing system may determine (i) the copy number of the common SMN gene(s) and (ii) the copy number of any intact SMN gene(s) using a mixture normal distribution model, and the first predetermined a posterior probability threshold considering (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 (or any given posterior probability threshold of a given description) may be or approximately be, 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 more. For example, the first predetermined posterior probability threshold may be 0.95.

[0111] Method 800 proceeds from step 816 to step 820, where the computing system determines, for one of a plurality of SMN1 gene-specific bases (also referred to herein as SMN differentiating bases) associated with the intact SMN1 gene, the most likely combination of the plurality possible combinations, each containing a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene, summed with the number of copies of any intact SMN gene(s), taking into account (a) the number of sequence reads (e.g. unnormalized or normalized number of sequence reads) from a plurality of sequence reads with bases that contain a SMN1 gene-specific base, and (b) the number of sequence reads (e.g., unnormalized or normalized number of sequence reads) from a plurality of sequence reads with bases that contain a SMN1 gene-specific SMN2 gene base of the SMN2 gene corresponding to a base specific to the SMN1 gene. The possible number of copies of the SMN1 gene may be or approximately be, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more. The possible number of copies of the SMN2 gene may be or approximately be, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more.

[0112] In some embodiments, the most likely combination of a possible SMN1 gene copy number and a possible SMN2 gene copy number is associated with the highest posterior probability, compared to other combinations of the plurality of combinations given (a) the number of sequence reads from the plurality of sequence reads with bases, that contain the SMN1 gene-specific base, and (b) the number of sequence reads from the plurality of sequence reads with bases that contain the corresponding SMN2 gene-specific base. The highest posterior probability (or any probability of a given description) may be or approximately be, 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 more. The difference in posterior probability (or any probability of a given description) may be or approximately be, 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 more.

[0113] In some embodiments, to determine the most likely combination of a possible number of copies of the SMN1 gene and a possible combination of the SMN2 gene, the computer system may determine the most likely combination of a plurality of possible combinations, each of which contains a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene. summed with the copy number of any defined intact SMN genes, taking into account the ratio of (a) the number of sequence reads from the set of sequence reads with bases that contain the SMN1 gene-specific base, and (b) the number of sequence reads from the set of sequence reads with bases that contain SMN2 a gene-specific base of the SMN2 gene corresponding to a specific base for the SMN1 gene. To determine the most likely combination of the possible number of copies of the SMN1 gene and the possible combination of the SMN2 gene; the computing system may determine (a) the number of sequence reads from a plurality of sequence reads with bases that contain an SMN1 gene-specific base, and (b) the number of sequence reads from a plurality of sequence reads with bases that contain an SMN2 gene-specific base of the SMN2 gene corresponding to a specific for the SMN1 gene base. The computing system may determine the ratio of (a) the number of sequence reads from the plurality of sequence reads with bases that contain the SMN1 gene-specific base, and (b) the number of sequence reads from the plurality of sequence reads with bases that contain the SMN2 gene-specific base of the SMN2 gene corresponding to a base specific to the SMN1 gene. The computing system may determine the most likely combination, from a plurality of possible combinations, each of which contains a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene, summed with the number of copies of any intact SMN gene, determined based on the ratio of (a) the number of sequence reads from the set sequence reads with bases that contain the SMN1 gene-specific base, and (b) the number of sequence reads from the plurality of sequence reads with bases that contain the SMN2 gene-specific base of the SMN2 gene corresponding to the SMN1 gene-specific base.

[0114] In some embodiments, to determine the most likely combination of a possible number of copies of the SMN1 gene and a possible combination of the SMN2 gene, the computing system determines, for each of the plurality of SMN1 gene-specific bases: the most likely combination of the plurality of possible combinations, each of which contains the possible number of copies the SMN1 gene and the possible copy number of the SMN2 gene, summed with the copy number of any particular intact SMN gene that is associated with the largest posterior probability given by (a) the number of sequence reads from the set of sequence reads with bases that contain the SMN1 gene-specific base, and ( b) the number of sequence reads from a plurality of sequence reads with bases that contain the SMN2 gene-specific base of the SMN2 gene corresponding to the SMN1 gene-specific base. The number of sequence reads aligned to an SMN1 gene-specific base (or an SMN2 gene-specific base) may be or approximately, for example, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or more. To determine the copy number of the SMN1 gene, the computing system may determine the copy number of the SMN1 gene based on the possible SMN1 gene copy number of the most likely combination of the possible SMN1 gene copy number and the possible SMN2 gene copy number determined for each of a plurality of bases specific to the SMN1 gene.

[0115] In some embodiments, the SMN1 gene-specific base is a splicing enhancer. The SMN1 gene-specific base may be the base at c.840 of the SMN1 gene. In some embodiments, the SMN1 gene-specific base has a match with each of a plurality of SMN1 gene-specific bases other than the SMN1 gene-specific base above a predetermined matching threshold. The predetermined compliance threshold (or any threshold value herein) may be or approximately be, for example, 80%, 81%, 82%, 83%, 84%, 85%, 86%, 87%, 88%, 89%, 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, 99% or more. For example, the matching threshold could be 97%. The plurality of bases specific to the SMN1 gene may contain or approximately contain 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21 or more specific for the SMN1 gene bases. For example, the set of bases specific to the SMN1 gene may contain 8 bases specific to the SMN1 gene. Each of a variety of bases specific to the SMN1 gene may be located on intron 6, exon 7, intron 7, or exon 8 of the SMN1 gene.

[0116] A plurality of bases specific to the SMN1 gene if the subject is a first race (or ethnicity), a plurality of bases specific to the SMN1 gene if the subject is a second race (or ethnicity), and a plurality of bases specific to the SMN1 gene , if the subject is of an unknown race, may vary. Race may 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 ethnicity) of the subject may be unknown, and the set of bases specific to the SMN1 gene may not be race (or ethnicity) specific. The race (or ethnicity) of the subject may be known, and a plurality of bases specific to the SMN1 gene may be specific to the race (or ethnicity) of the subject. In some embodiments, the computing system may receive information about the race (or ethnicity) of a subject. Based on the received race (or ethnicity) information, the computing system may select sets of bases specific for the SMN1 gene from the plurality of bases specific for the SMN1 gene.

[0117] Method 800 proceeds from step 820 to step 824, where the computing system determines the copy number of the SMN1 gene using the most likely combination of the possible SMN1 gene copy number and the possible SMN2 gene copy number determined for the SMN1 gene-specific base. In an alternative or additional embodiment, the computing system determines the copy number of the SMN2 gene using the most likely combination of the possible SMN1 gene copy number and the possible SMN2 gene copy number determined for the SMN1 gene-specific base.

[0118] In some embodiments, to determine the copy number of the SMN1 gene, the computer system may determine the copy number of the SMN1 gene and the copy number of the SMN2 gene 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 each of the plurality of bases, specific for the SMN1 gene. To determine the copy number, the computing system may determine the copy number of the SMN1 gene using the most likely combination of the possible copy number of the SMN1 gene and the possible copy number of a particular SMN2 gene base specific for the SMN1 gene and a second predetermined posterior probability threshold for the combination of the possible copy number of the SMN1 gene and possible number of copies of the SMN2 gene. The second predetermined posterior probability threshold (or any predetermined posterior probability threshold of this description) may be or approximately be, 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 more. For example, the second predefined posterior probability threshold may be 0.6 or 0.8.

[0119] In some embodiments, the majority of the determined possible copy numbers of the SMN1 gene are the same. The copy number of a particular SMN1 gene may be consistent with the possible copy number of the SMN1 gene. The computing system may determine a possible combination containing a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene, summed with the number of copies of any intact SMN gene, determined by taking into account (a) the number of sequence reads from a set of sequence reads with bases that contain any of the set bases specific to the SMN1 gene, and (b) the number of sequence reads from a plurality of sequence reads with bases that contain any of a plurality of corresponding bases specific to the SMN2 gene. The computing system can determine the possible number of copies of a possible combination - the agreed upon possible number of copies of the SMN1 gene.

[0120] In some embodiments, to determine the copy number of the SMN1 gene, the computing system may determine the copy number of the SMN1 gene to be zero, one, or greater than one. In some embodiments, the computing system may determine a subject's spinal muscular atrophy (SMA) status based on the copy number of the SMN1 gene. The SMA status for a subject may include SMA, SMA carrier/no SMA carrier, and non-SMA carrier. In some embodiments, a computing system may determine a subject to be a silent carrier of SMA using a set of sequence reads from a plurality of sequence reads aligned to g.27134 of the SMN1 gene, and based on sequence reads aligned to g.27134 of the SMN1 gene.

[0121] For example, the computing system at step 820 may determine the ratio of SMN1 to SMN2 reads of overlapping locations where genes have different sequence bases. For positions at which SMN1 differs from SMN2 , the computing system may produce reads that overlap based on either SMN1 or SMN2 . Based on these reads, the computing system can count the number of SMN1- specific bases and the number of SMN2- specific bases. The computing system can determine the proportion of SMN1 or SMN2 reads. The computing system can calculate the number of copies of SMN1 and SMN2 at positions where SMN1 differs from SMN2 . The computing system may combine the number of full-length copies with the ratio of SMN1 to SMN2 to determine the copy number of SMN1 and SMN2 . The computing system at step 824 may combine the copy number from the plurality of fixed differences between SMN1 and SMN2 to obtain the exact copy number of SMN1 and SMN2 .

[0122] To determine SMA/no SMA or carrier/no carrier . In some embodiments, the computing system may determine the copy number of the truncated SMN gene(s) using (i) the copy number of the total specific SMN gene(s) and (ii) the copy number intact SMN gene(s). The copy number of the truncated SMN gene(s) may be the difference of (i) the copy number of the common SMN gene(s) and (ii) the copy number of the intact SMN gene(s ) .

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

[0124] Method 800 ends at step 828.

Genotyping of the cytochrome P450 family 2 subfamily D member 6 gene using sequencing data

[0125] In FIG. 9 is a flow diagram showing an example method 900 of genotyping a cytochrome P450 family 2 subfamily D member 6 gene using sequencing data, such as whole genome sequencing data. Method 900 may be implemented as a set of executable program instructions stored on a computer-readable medium, such as one or more computer system disk drives. For example, the computing system 1100 shown in FIG. 11 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 storage device, such as RAM, and executed by one or more processors of the computing system 1100. Although method 900 is described with respect to the computing system 1100 shown in FIG. 11, the description is provided for illustrative purposes only and is not intended to be limiting. In some embodiments, method 900 or portions thereof may be performed sequentially or in parallel by multiple computing systems.

[0126] The number of sequence reads (eg, unnormalized or normalized number of sequence reads) mapped to a CYP2D6 gene or a CYP2D7 gene can be used to determine the total copy number (CN) of the CYP2D6 gene and the CYP2D7 gene, using a mixture normal distribution model. The total copy number of the CYP2D6 gene and the CYP2D7 gene can be used to determine the copy number of CYP2D6 at different CYP2D6 / CYP2D7 differentiating bases (also referred to herein as CYP2D6 gene-specific bases) by repeating all possible combinations of CYP2D6 and CYP2D7 copy numbers at different CYP2D6/CYP2D7 differentiating bases. CYP2D7 . CYP2D6 copy number at different CYP2D6 / CYP2D7 differentiating bases can be used to designate structural variants. For example, at each of the CYP2D6/CYP2D7 differentiating bases (also referred to herein as CYP2D6 gene-specific bases), the number of chromosomes carrying the CYP2D6 gene and the number of chromosomes carrying the CYP2D7 gene can be determined by combining the total copy number of the CYP2D6 gene and the CYP2D7 gene with the number of reads supporting a base specific to each gene. Based on the determined total copy number, all possible combinations of CYP2D6 and CYP2D7 copy numbers can be repeated to obtain the combination that gives the highest posterior probability for the observed CYP2D6 and CYP2D7 read counts. Structural variants can be identified by identifying the bases at which the copy number of the CYP2D6 gene changes.

[0127] One or more minor variants may be defined. Minor variants can be determined, for each minor variant minor variant position, by repeating all possible copy number combinations of the variant allele and the reference (non-variant) allele to determine the most likely allele copy number using the sequence read with the minor variant position in the CYP2D6 or CYP2D7 gene. For example, if there are only three copies of the CYP2D6 gene and 10 reads of the variant allele, and 20 reads of the reference allele, the copy number of the variant allele can be determined to be one, i.e., there is one copy of the CYP2D6 gene carrying the minor variant. For example, you can search for small variants that define star alleles in sequencing data (for example, a BAM file). The small variants of interest can be divided into variants that fall within the homologous regions of CYP2D 6/CYP2D7 and variants that do not. For the former, variant reads aligned to the CYP2D6 gene or the CYP2D7 gene, overlapping each position of the small variant of the CYP2D6 gene of interest or the corresponding position in the CYP2D7 gene, can be found. In the latter case, reads aligned to the CYP2D6 gene and overlapping the position of a small variant of the CYP2D6 gene of interest can be found. The number of copies defined in an area can also be taken into account when defining small variants. Specific structural variants and minor variants can be compared to star allele definitions to identify star alleles that can be further grouped into haplotypes.

[0128] After method 900 begins at step 904, method 900 proceeds to step 908, where a computing system (e.g., computing system 1100 described with reference to FIG. 11) determines (i) a first number of sequence reads from a plurality of sequence reads , aligned with the cytochrome P450 family 2 subfamily D member 6 gene ( CYP2D6 ) or the cytochrome P450 family 2 subfamily D member 7 gene ( CYP2D7 ). The first number of sequence reads aligned to the first CYP2D6 gene or CYP2D7 gene (or any gene herein) may be or approximately 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 more).

[0129] The computing system may receive sequencing data containing multiple sequence reads obtained from a sample from a subject aligned to the CYP2D6 gene or the CYP2D7 gene. In some embodiments, the sequencing data comprises whole genome sequencing (WGS) data or short-read WGS data. The subject may be a fetal subject, a neonatal subject, a pediatric subject, an adolescent subject, or an adult subject. The sample may contain cells or extracellular DNA. The sample may contain cells or extracellular DNA.

[0130] In some embodiments, a sequence read from a plurality of sequence reads aligned to a CYP2D6 gene or a CYP2D7 gene with an alignment quality score of about zero. The alignment quality may be, for example, 0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10 or more (on a scale of 0 to 1 from the alignment score).

[0131] In some embodiments, to determine (i) a first number of sequence reads from a plurality of sequence reads aligned to a CYP2D6 gene or a CYP2D7 gene, a computing system may determine (i) a first number of sequence reads from a plurality of sequence reads aligned to at least one exon or intron of the CYP2D6 gene (for example, one of exons 1-9 or one of introns 1-8 of the CYP2D6 gene) and/or at least one of the exons or introns of the CYP2D7 gene (for example, one of exons 1-9 or one from introns 1-8 of the CYP2D7 gene).

[0132] Method 900 proceeds from step 908 to step 912, where the computing system determines (i) the first normalized number of 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 number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene (or any gene herein) may be or approximately be, for example, 1, 2, 3, 4, 5, 6, 7, 9, 10, 20, 30 , 40, 50, 60, 70, 80, 90, 100 or more. The length of the CYP2D6 gene may be or approximately 4.4 kb, for example. The length of the CYP2D7 gene may be or approximately 4.9 kb, for example.

[0133] In some embodiments, to determine (i) the first normalized read count of the sequences aligned to the CYP2D6 gene or the CYP2D7 gene, the computing system may determine (i) the first normalized read count of the sequences aligned to the CYP2D6 gene or the CYP2D7 gene using ( i) the length of the CYP2D6 gene or the CYP2D7 gene, respectively, and (iii) the sequence read depth of a region of the subject's genome other than the genetic loci containing the CYP2D6 gene and the CYP2D7 gene as determined by sequencing. The sequence read depth of a region of a subject's genome, excluding genetic loci containing the CYP2D6 gene and the CYP2D7 gene (or any genes herein), in the sequencing data may be or approximately equal to, for example, 3, 4, 5, 10, 20, 30, 40 , 50, 60, 70, 80, 90, 100 or more.

[0134] To determine (i) a first normalized number of sequence reads aligned to the CYP2D6 gene or a CYP2D7 gene, and (ii) a second normalized number of sequence reads aligned to the second region, the computing system may determine (i) a first normalized number of reads by length CYP2D6 gene or CYP2D7 gene sequences aligned to the CYP2D6 gene or CYP2D7 gene, using (i) the length of the CYP2D6 gene or CYP2D7 gene, respectively. The computing system may determine (i) a first normalized read depth of sequences aligned to a CYP2D6 gene or a CYP2D7 gene from (i) a length-normalized read of a CYP2D6 or CYP2D7 gene, using the read depth of sequences of a region of the subject's genome other than the genetic loci containing the gene CYP2D6 and CYP2D7 . The first normalized read depth of sequences aligned to a CYP2D6 gene or a CYP2D7 gene may be the first normalized read depth of sequences aligned to a CYP2D6 gene or a CYP2D7 gene, respectively.

[0135] In some embodiments, to determine (i) the first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, the computing system may determine (i) the first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene using ( i) the GC content of the CYP2D6 gene or the CYP2D7 gene and (iii) the depth of sequence reads in the subject's genomic region, excluding genetic loci containing the CYP2D6 gene and the CYP2D7 gene, as determined by sequencing, and (iv) the GC content of the genomic region. The GC content of the CYP2D6 gene or the CYP2D7 gene (or any gene herein) may be or approximately be, 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 sequence read depth of a region of a subject's genome, excluding genetic loci containing the CYP2D6 gene and the CYP2D7 gene, as determined by sequencing may be or approximately equal to, for example, 3, 4, 5, 10, 20, 30, 40, 50, 100, or more. The GC content of a region of a subject's genome, excluding genetic loci containing the CYP2D6 gene and the CYP2D7 gene (or any genes herein), as determined by sequencing may be or approximately be, for example, 40%, 41%, 42%, 43%, 44 %, 45%, 46%, 47%, 48%, 49%, 50%, 51%, 52%, 53%, 54%, 55%, 56%, 57%, 58%, 59% or 60%.

[0136] The region depth may be the average sequence read depth of a region of the subject's genome, excluding the genetic loci containing the CYP2D6 gene and the CYP2D7 gene, as determined by sequencing data. Region depth may include the median sequence read depth of a region of a subject's genome, excluding genetic loci containing the CYP2D6 gene and the CYP2D7 gene, as determined by sequencing data. The depth of the region may be or be approximately 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or more. The region length may comprise approximately 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000 or more preselected regions of approximately 0.5 kb, 1 kb, 1.5 k in length. bp, 2 kb, 2.5 kb or 3 kb. in each throughout the subject's genome. For example, a region may contain approximately 3000 preselected regions of approximately 2 kb in length. in each throughout the subject's genome.

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

[0138] Method 900 proceeds from step 912 to step 916, where the computing system determines (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene using a mixture normal distribution model containing a plurality of distributions, each representing a different integer copy number, given ( i) the first normalized number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene. The total copy number of the CYP2D6 gene and the CYP2D7 gene (or any genes herein) may be or be about 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more.

[0139] In some embodiments, the mixture normal distribution model is a one-dimensional mixture normal distribution model. The set of normal distribution model distributions can represent integer copy numbers, for example, 0 to 5, 0 to 6, 0 to 7, 0 to 8, 0 to 9, 0 to 10, 0 to 11, 0 to 12, 0 to 13, 0 to 14, or 0 to 15. For example, a set of mixture normal distribution model distributions can represent integer copy numbers from 0 to 10. The mean (for example, 1, 2, 3, 4, 5 , 6, 7, 8, 9, 10, or more) of each of the plurality of distributions may be an integer number of copies represented by the distribution (e.g., copy number 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more). The standard deviation of the distributions may be or approximately be, for example, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 or more .

[0140] In some embodiments, to determine (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene, the computing system may determine (i) the total copy number of the CYP2D6 gene and the CYP2D7 gene using a mixture of normal distributions model and a first predefined posterior probability threshold, given (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 given posterior probability threshold of a given description) may be or approximately be, 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 more. For example, the first predetermined posterior probability threshold may be 0.95.

[0141] Method 900 proceeds from step 916 to step 920, where the computing system determines, for one of a plurality of CYP2D6 gene-specific bases (also referred to herein as CYP2D6 / CYP2D7 differentiating bases), the most likely combination from a plurality of possible combinations, each of which contains a possible number of copies of a CYP2D6 gene and a possible number of copies of a CYP2D7 gene, summed with the total number of copies of a CYP2D6 gene and a particular CYP2D7 gene, taking into account (a) the number of sequence reads (e.g., unnormalized or normalized number of sequence reads) from the set of sequence reads with bases that contain a CYP2D6 gene-specific base, and (b) the number of sequence reads (eg, non-normalized or normalized number of sequence reads) from the set of sequence reads with bases that contain a CYP2D7 gene-specific base corresponding to a CYP2D6 gene-specific base. The possible number of copies of the CYP2D6 gene may be or be about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more. The possible number of copies of the CYP2D7 gene may be or be about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more.

[0142] In some embodiments, the most likely combination of a possible CYP2D6 gene copy number and a possible CYP2D7 gene copy number is associated with the highest posterior probability, compared to other combinations of the plurality of combinations given (a) the number of sequence reads from the plurality of sequence reads with bases, that contain a CYP2D6 gene-specific base, and (b) the number of sequence reads from a plurality of sequence reads with bases that contain the corresponding CYP2D7 gene-specific base. The highest posterior probability (or any probability of a given description) may be or approximately be, 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 more. The difference in posterior probability (or any probability of a given description) may be or approximately equal to, 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 more.

[0143] In some embodiments, to determine the most likely combination containing a possible number of copies of the CYP2D6 gene and a possible number of copies of the CYP2D6 gene, the computing system may determine the most likely combination from a plurality of possible combinations, each of which contains a possible number of copies of the CYP2D6 gene and a possible number of copies of the gene CYP2D7 , summed with the total copy number of the CYP2D6 gene and a specific CYP2D7 gene, taking into account the ratio of (a) the number of sequence reads in multiple sequence reads with bases that contain a CYP2D6 gene-specific base, and (b) the number of sequence reads in multiple sequence reads with bases , which contain a CYP2D7 gene-specific base corresponding to a CYP2D6 gene-specific base. To determine the most likely combination involving the possible copy number of the CYP2D6 gene and the possible copy number, the computing system can determine (a) the number of sequence reads from a set of sequence reads with bases that contain a CYP2D6 gene-specific base, and (b) from a set of sequence reads with bases that contain a CYP2D7 gene-specific base corresponding to a CYP2D6 gene-specific base. The computing system may determine the ratio of (a) the number of sequence reads from a plurality of sequence reads with bases that contain a CYP2D6 gene-specific base, and (b) the number of sequence reads from a plurality of sequence reads with bases that contain a CYP2D7 gene-specific base corresponding to a CYP2D7 gene-specific base. for the CYP2D6 gene base. The computing system may determine the most likely combination from a plurality of possible combinations, each containing a possible number of copies of the CYP2D6 gene and a possible number of copies of the CYP2D7 gene, summed with the total number of copies of the CYP2D6 gene and the particular CYP2D7 gene, taking into account the ratio of (a) the number of sequence reads from a plurality of sequence reads with bases that contain a CYP2D6 gene-specific base, and (b) a number of sequence reads from the plurality of sequence reads with bases that contain a CYP2D7 gene-specific base corresponding to a CYP2D6 gene-specific base.

[0144] In some embodiments, to determine the most likely combination of a possible number of copies of a CYP2D6 gene and a possible combination of a CYP2D7 gene, the computing system determines, for each of a plurality of bases specific to the CYP2D6 gene: the most likely combination of a plurality of possible combinations, each of which contains a possible number copies of the CYP2D6 gene and the possible number of copies of the CYP2D7 gene, summed with the total number of copies of a particular CYP2D6 gene and the CYP2D7 gene, which is associated with the largest posterior probability given by (a) the number of sequence reads from the set of sequence reads with bases that contain the CYP2D6 gene-specific base , and (b) the number of sequence reads from a plurality of sequence reads with bases that contain a CYP2D7 gene-specific base of the CYP2D7 gene corresponding to a CYP2D6 gene-specific base. The number of sequence reads aligned to an SMN1 gene-specific base (or an SMN2 gene-specific base) may be or approximately, for example, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 or more. To determine the allele of the CYP2D6 gene in a subject, the computing system may determine the allele of the CYP2D6 gene that the subject has a minor variant, a structural variant of the CYP2D6 gene, or neither, 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 each of the many bases specific to the CYP2D6 gene.

[0145] In some embodiments, the CYP2D6 gene-specific base has a match with each of a plurality of CYP2D6 gene-specific bases, except for the CYP2D6 gene-specific base, above a predetermined match threshold. The compliance threshold (or any threshold of this description) may be or be approximately, for example, 80%, 81%, 82%, 83%, 84%, 85%, 86%, 87%, 88%, 89%, 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, 99% or more. For example, a predefined compliance threshold could be 97%. The plurality of bases specific to the CYP2D6 gene may contain or contain approximately, for example, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 118, 120, 130, 140, 150, 160, 170 or more specific bases for the CYP2D6 gene. For example, a set of bases specific for the CYP2D6 gene may contain 118 bases specific for the CYP2D6 gene.

[0146] Method 900 proceeds from step 920 to step 924, in which the computing system determines one or more structural variants of the CYP2D6 gene present in the subject 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 base, specific for the CYP2D6 gene. For example, a computing system may determine the ratio of CYP2D6 to CYP2D7 reads overlapping locations where genes have different sequence bases. For positions at which CYP2D6 differs from CYP2D7 , the computing system may generate reads that overlap at either a CYP2D6 or CYP2D7 base. Based on these reads, the computing system can count the number of CYP2D6- specific bases and the number of CYP2D7- specific bases. The computing system can determine the proportion of CYP2D6 or CYP2D7 reads. A computing system can calculate the copy number of CYP2D6 and CYP2D7 at positions where CYP2D6 differs from CYP2D7 . A computing system may combine the total copy number of CYP2D6 and CYP2D7 with the ratio of CYP2D6 to CYP2D7 to determine the copy number of CYP2D6 and CYP2D7 . The computing system may identify minor variants using the copy numbers of CYP2D6 and CYP2D7 with one or more fixed differences between CYP2D6 and CYP2D7 . The computational system may detect structural variants by combining the copy numbers of CYP2D6 and CYP2D7 with multiple fixed differences between CYP2D6 and CYP2D7 to determine whether there is a transition between the copy numbers of CYP2D6 and CYP2D7 that determines the type of structural variant found in the sample.

[0147] REP-containing fusion genes . In some embodiments, the computing system may determine (ii) a second number of sequence reads from a plurality of sequence reads aligned to the spacer region between the CYP2D7 gene and the REP7 repeat element downstream of the CYP2D7 gene. The second number of sequence reads from the plurality of sequence reads aligned to the spacer region between the CYP2D7 gene and the REP7 repeat element downstream of the CYP2D7 gene may 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 more. The computing system may determine (ii) a second normalized read count of the sequence aligned to the spacer region using (ii) the length of the spacer region. The second normalized number of sequence reads aligned to the spacer region may be or approximately, for example, 1, 2, 3, 4, 5, 6, 7, 9, 10, 20, 30, 40, 50, 60, 70, 80 , 90, 100 or more. The length of the spacer region may be or be about, for example, 1.5 kb. The computing system may determine (ii) the copy number of the spacer region using a mixture normal distribution model defined by (ii) a second normalized number of sequence reads aligned to the spacer region. The number of copies of the spacer region may be or be about, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more. To determine the allele of the CYP2D6 gene present in a subject, the computing system may determine the allele of the CYP2D6 gene present in the subject, a minor variant or a structural variant of the CYP2D6 gene, or neither, using a combination of the possible copy number of the CYP2D6 gene and the possible copy number of the CYP2D7 gene. defined for the base specific for the CYP2D6 gene and the copy number of the spacer region. The structural variant may contain a CYP2D6 / CYP2D7 fusion allele with a spacer region and a REP7 repeat element downstream of the CYP2D6 / CYP2D7 fusion allele.

[0148] Method 900 proceeds from step 924 to step 928, in which the computing system may, for the position of the minor variant of the CYP2D6 gene associated with the minor allele of the variant of the CYP2D6 gene, determine the most likely combination of the possible number of copies of the minor allele of the CYP2D6 gene at the position of the minor variant and the possible copy number of the reference allele of the CYP2D6 gene at the minor variant position summed with the copy number of the CYP2D6 gene at the minor variant position, taking into account (a) the number of sequence reads (e.g., unnormalized or normalized number of sequence reads) aligned to the CYP2D6 gene overlapping with the minor variant position , and the base that contains the minor variant allele of the CYP2D6 gene at the position of the minor variant, and (b) the number of sequence reads (e.g., unnormalized or normalized number of sequence reads) aligned to the CYP2D6 gene overlapping with the position of the minor variant, and the base that contains the reference allele of the CYP2D6 gene in the position of the minor variant. The possible copy number of the small variant allele of the CYP2D6 gene, the most likely combination at the position of the small variant may indicate one or more small variants of the CYP2D6 gene.

[0149] The computing system may determine, for each of a plurality of minor variant positions of the CYP2D6 gene, the position of the minor variant associated with the minor variant allele of the CYP2D6 gene, determining the most likely combination of the possible copy number of the minor allele of the CYP2D6 gene at the position of the minor variant and the possible copy number of the reference allele of the CYP2D6 gene at the position of the minor variant, summed with the copy number of the CYP2D6 gene at the position of the minor variant, taking into account (a) the number of sequence reads (e.g., unnormalized or normalized number of sequence reads) aligned to the CYP2D6 gene overlapping with the position of the minor variant, and the base that contains the minor variant allele of the CYP2D6 gene at the minor variant position and (b) the number of sequence reads (e.g., non-normalized or normalized number of sequence reads) aligned to the CYP2D6 gene overlapping with the minor variant position and the base that contains the reference allele of the CYP2D6 gene at the minor position option. The possible copy number of small variant alleles of the CYP2D6 gene, the most likely combinations at multiple positions of the small variants may indicate one or more small variants of the CYP2D6 gene.

[0150] In some embodiments, the computing system may determine 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 CYP2D6 gene copy number at the minor variant position may include the CYP2D6 gene copy number of the possible CYP2D6 gene copy numbers of the most likely specific combinations. The CYP2D6 gene copy number at the minor variant position may include the CYP2D6 gene copy number of the possible CYP2D6 gene copy numbers of the most likely specific combinations and located closest to the minor variant position. The number of copies of the CYP2D6 gene at the minor variant position may include the number of copies of the CYP2D6 gene at the 5' position or at the 3' position of the small variant position.

[0151] In some embodiments, the computing system may (a) determine the number of sequence reads (eg, unnormalized or normalized number of sequence reads) with bases that contain a small variant allele of the CYP2D6 gene. The computing system may (b) determine the number of sequence reads (eg, unnormalized or normalized number of sequence reads) to bases that contain the reference allele of the CYP2D6 gene.

[0152] The method 900 proceeds from step 928 to step 932, where the computing system determines one or more small variants of the CYP2D6 gene using the possible number of copies of the small variant allele of the CYP2D6 gene of the most likely determined combination. The computing system may determine one or more small variants of the CYP2D6 gene using the possible copy number of small variant alleles of the CYP2D6 gene, the most likely combinations at a plurality of specific positions of the small variants.

[0153] In some embodiments, the position of the minor variant is in the homologous region of CYP2D6 / CYP2D7 . To determine the most likely combination, the computing system may determine the most likely combination of the possible number of copies of the minor variant allele of the CYP2D6 gene at the minor variant position and the possible number of copies of the reference allele of the CYP2D6 gene at the minor variant position, summed with the number of copies of the CYP2D6 gene at the minor variant position, taking into account ( a) the number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, with a base that contains the minor variant allele of the CYP2D6 gene at the position of the minor variant, and/or (b) the number of sequence reads aligned to the CYP2D6 gene or the CYP2D7 gene, with the base, which contains the CYP2D6 reference allele at the minor variant position. In some embodiments, the position of the minor variant is not in the CYP2D6 / CYP2D7 homologous region. To determine the most likely combination, the computing system may determine the most likely combination of the possible number of copies of the minor variant allele of the CYP2D6 gene at the minor variant position and the possible number of copies of the reference allele of the CYP2D6 gene at the minor variant position, summed with the number of copies of the CYP2D6 gene at the minor variant position, taking into account ( a) the number of sequence reads aligned to the CYP2D6 gene rather than the CYP2D7 gene with a base that contains the minor variant allele of the CYP2D6 gene at the position of the minor variant, and/or (b) the number of sequence reads aligned to the CYP2D6 gene rather than the CYP2D7 gene CYP2D7 , with a base that contains the CYP2D6 reference allele at the minor variant position.

[0154] For example, the computing system may first determine the nature and breakpoints of the SV (structural variant, e.g., deletion or duplication) based on the copy number of paralog-specific bases. Additionally or alternatively, the computing system can subsequently identify a predefined set of small variants (those are variants specific to a gene of interest, e.g., CYP2D6 , and are different from the differentiating bases of paralogs), based on read alignment, total copy number, and (sometimes ) the nature of the SV and the break points determined in the first stage. Because the alignment is not always precise, the computational system can extract bases of interest from reads that match any paralog.

[0155] Method 900 proceeds from step 932 to step 936, where the computing system can determine the star allele and/or haplotype of the CYP2D6 gene present in the subject using one or more defined structural variants of the CYP2D6 gene, and/or one or more minor gene variants CYP2D6 . The star allele may be associated with a known function. The star allele and/or haplotype of the CYP2D6 gene may 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.

[0156] Enzyme activity . In some embodiments, the computing system may determine the level of CYP2D6 enzymatic activity in a subject using a specific allele of the CYP2D6 gene. Enzyme activity can be low, medium, normal or ultra-fast. The computing system may determine a treatment dosage recommendation and/or a treatment recommendation for a subject based on one or more minor variants and/or one or more structural variants.

[0157] Method 900 ends at step 940.

Genotyping paralogs using sequencing data

[0158] In FIG. 10 is a flow diagram showing an example method 1000 for genotyping a paralog using sequencing data, such as whole genome sequencing data. Method 1000 may be implemented as a set of executable program instructions stored on a computer-readable medium, such as one or more computer system disk drives. For example, the computing system 1100 shown in FIG. 11 and described in more detail below, may execute a set of executable program instructions to implement method 1000. When method 1000 is initiated, the executable program instructions may be loaded into a storage device, such as RAM, and executed by one or more processors of the computer system 1100. While method 1000 is described with respect to the computing system 1100 shown in FIG. 11, the description is provided for illustrative purposes only and is not intended to be limiting. In some embodiments, the method 1000 or portions thereof may be performed sequentially or in parallel by multiple computing systems.

[0159] After method 1000 begins at step 1004, method 1000 proceeds to step 1008, where a computing system (e.g., computing system 1100 described with reference to FIG. 11) receives sequencing data containing a plurality of sequence reads obtained from the sample from the subject aligned with the first paralogue or the second paralogue. Techniques for generating sequence reads include sequencing by synthesis using, for example, the MINISEQ, MISEQ, NEXTSEQ, HISEQ, and NOVASEQ sequencing tools from Illumina, Inc. (San Diego, California).

[0160] Method 1000 proceeds from step 1008 to step 1012, where the computing system determines the copy number of the first type of paralogs using a mixture normal distribution model containing a plurality of distributions, each representing a different integer copy number given by (i) the first number of sequence reads , aligned with the first area. The number of copies of paralogues of the first type (or any type of this description) may be or be approximately, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more.

[0161] The set of distributions of a mixture of normal distributions model may 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 , from 0 to 12, from 0 to 13, from 0 to 14, or from 0 to 15. For example, a set of distributions of a mixture normal distribution model can represent integer copy numbers from 0 to 10. The mean value (for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or more) of each of the plurality of distributions may be an integer number of copies represented by the distribution (e.g., copy number 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more). The standard deviation of the distributions may be or approximately be, for example, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 or more .

[0162] In some embodiments, the computing system may determine (i) a first number of sequence reads from a plurality of sequence reads in sequencing data obtained from a sample from a subject aligned to the first region. The first number of sequence reads of the plurality of sequence reads in the sequencing data obtained from a sample from a subject aligned to the first region (or any region herein) may be or approximately be, 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 more. The computing system may determine (i) a first normalized number of sequence reads aligned to the first region using (i) the length of the first region. The first normalized number of sequence reads aligned to the first region (or any region herein) may 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 more. The length of the first region may be or be approximately, 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 t. bp, 19 kbp, 20 kbp, 21 kbp, 22 kbp, 23 kbp, 24 kbp. b., 25 kb., 26 kb., 27 kb., 28 kb., 29 kb., 30 kb. or more. To determine the copy number of the first type of paralogs, the computing system may determine the copy number of the first type of paralog using a mixture normal distribution model defined by (i) the first normalized number of sequence reads aligned to the first region.

[0163] In some embodiments, the computing system may determine the copy number of one or more second type paralogues using a mixture of normal distributions defined by (ii) a second number of sequence reads aligned to the second region. The number of copies of one or more paralogs of the second type (or any type of this description) may be or be approximately, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more. To determine the copy number or allele of the first paralog, the computing system may determine the copy number or allele of the first paralog using the most likely combination of the possible copy number of the first paralog and the possible copy number of the second paralog, determined for a base specific to the first paralog, and the copy number of one or more paralogues of the second type. The computing system can determine the number of copies of the third type of paralogs from the number of copies of the first type of paralogs and the number of copies of the second type of paralogs. The number of copies of the third type of paralogues (or any type of this description) may be or be approximately, for example, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 or more. To determine the copy number or allele of the first paralog, the computing system may determine the copy number or allele of the first paralog using the most likely combination of the possible copy number of the first paralog and the possible copy number of the second paralog, determined for a base specific to the first paralog.

[0164] Methods for aligning sequences to a reference genomic sequence may use transformers such as the Burrows-Wheeler transformer (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.

[0165] Method 1000 proceeds from step 1012 to step 1016, where the computing system determines, for one set of bases specific to the first paralog, from a plurality of possible combinations, each of which contains a possible number of copies of the first paralog of the first type and a possible number of copies of the second paralog of the first type, summed with the number of copies of certain paralogs of the first type, given (a) the number of sequence reads (e.g., the unnormalized or normalized number of sequence reads) from the set of sequence reads with bases that contain a base specific to the first paralog, and (b) the number of reads sequence reads (eg, unnormalized or normalized number of sequence reads) from a set of sequence reads with bases that contain a second paralog-specific base corresponding to a first paralog-specific base.

[0166] Method 1000 proceeds from step 1016 to step 1020, where the computing system determines the copy number or allele of the first paralog using the most likely combination of the possible copy number of the first paralog and the possible copy number of the second paralog determined for the base specific to the first paralog.

[0167] In some embodiments, the first paralogue is the survival motor neuron 1 ( SMN1 ) gene. The second paralogue may be the survival motor neuron 2 ( SMN2 ) gene. The first region may comprise at least one exon 1 to 6 of the SMN1 gene and at least one exon 1 to 6 of the SMN2 gene. The second region may comprise at least exon 7 or 8 of the SMN1 gene and at least one of exon 7 or 8 of the SMN2 gene. The first type of paralogues may include an intact SMN1 gene and an intact SMN2 gene. One or more second type paralogs may include an intact SMN1 gene, an intact 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. The computing system may determine the copy number of the SMN1 gene by implementing the method 800 (or a portion thereof) described with reference to FIG. 8.

[0168] In some embodiments, the first paralog is a cytochrome P450 family 2 subfamily D member 6 gene ( CYP2D6 ). The second paralog may be the cytochrome P450 family 2 subfamily D member 7 gene ( CYP2D7 ). The first region may contain the CYP2D6 gene and the CYP2D7 gene. The second region may contain a spacer region between the CYP2D7 gene and the REP7 repeat element downstream of the CYP2D7 gene. Paralogues of the first type may contain the CYP2D6 gene and the CYP2D7 gene. One or more type two paralogues may contain a CYP2D6 / CYP2D7 fusion allele with a spacer region and a REP7 repeat element downstream of the CYP2D6/CYP2D7 fusion allele. The first paralog allele may be an allele of the CYP2D6 gene present in the subject that is a minor variant or structural variant of the CYP2D6 gene. The computing system can determine the allele of the CYP2D6 gene by implementing the method 900 (or part thereof) described with reference to FIG. 9.

[0169] In various embodiments, the first and second paralogs may be different from each other. Examples of the first and second paralogs include, but are not limited to, the SMN1 gene and the SMN2 gene; CYP2D6 gene and CYP2D7 gene; double homeobox 4 gene (DUX4), DUX4c gene, DUX4-like protein 2 gene (DUX4L2), DUX4-like protein 3 gene (DUX4L3), DUX4-like protein 4 gene (DUX4L4), DUX4-like protein 5 gene (DUX4L5) , DUX4-like protein 6 (DUX4L6) gene, DUX4-like protein 7 (DUX4L7) gene, and double homeobox 2 (DUX2) gene; and the ribosomal protein S17 gene (RpS17) and the RpS17-like protein gene (RpS17L). In some embodiments, the computing system may determine the copy number or allele of the first paralog by implementing the method 800 (or a portion thereof) described with reference to FIG. 8, and/or the method 900 (or part thereof) described with reference to FIG. 9.

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

[0171] Method 1000 ends at step 1024.

Execution conditions

[0172] In FIG. 11 depicts the general architecture of an exemplary computing device 1100 configured to perform paralog genotyping. The general architecture of the computing device 1100 shown in FIG. 11 includes an arrangement of computer hardware and software components. Computing device 1100 may include many more (or fewer) elements than those shown in FIG. 11. However, it is not necessary that all of these generally conventional elements be shown to provide a description of the invention. As illustrated, 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 can communicate with each other via a communication bus. Network interface 1120 may provide connectivity to one or more networks or computing systems. Thus, the data processing unit 1110 can receive information and instructions from other computing systems or services via the network. The data processing unit 1110 may also communicate to and from memory 1170 and optionally provide output information to an optional display 1150 through an input/output device interface 1140. The input/output device interface 1140 may also receive input from an optional input device 1160, such as a keyboard. , mouse, digital pen, microphone, touch screen, gesture recognition system, voice recognition system, gamepad, accelerometer, gyroscope, or other input device.

[0173] Memory 1170 may contain computer program instructions (grouped as modules or components in some embodiments) that processing unit 1110 executes to implement one or more embodiments. Memory 1170 essentially includes RAM, ROM, and/or other read-only, non-transitory, or nonvolatile computer-readable media. Storage device 1170 may store an operating system 1172 that provides computer program instructions for use by processing unit 1110 in general administration and operation of computing device 1100. Storage device 1170 may further include computer program instructions and other information to implement aspects of this disclosure.

[0174] 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 with reference to FIG. 10. Alternatively or additionally, the paralog genotyping module 1174 may be or may include a module for determining SMN1 copy number using sequencing data, such as the method 800 described with reference to FIG. 8. Alternatively or additionally, paralog genotyping module 1174 may be or may include a module for genotyping the CYP2D6 gene using sequencing data, such as method 900 described with reference to FIG. 9. In addition, memory 1170 may include or communicate with data store 1190 and/or one or more other data stores that store sequencing data and/or genotyping results of one or more paralogs.

Examples

[0175] The following examples describe in more detail certain aspects described above, which are in no way intended to limit the scope of this disclosure.

Example 1

Diagnosis of spinal muscular atrophy and carrier screening based on whole genome sequencing data

[0176] Spinal muscular atrophy (SMA), caused by loss of the functional SMN1 gene but retention of the paralogous SMN2 gene, is a leading genetic cause of death in early childhood. Due to the nearly identical sequences of SMN1 and its paralogue SMN2 , analysis of this region using next-generation sequencing (NGS)-based assays is challenging. The American College of Medical Genetics recommends screening prospective parents for SMA before conception to determine the copy number (CN) of SMN1 .

[0177] This example describes a bioinformatics method that accurately identifies copy numbers of SMN1 and SMN2 using whole genome sequencing (WGS) data. The method calculates the copy number of SMN1 and SMN2 using read depth and eight informative reference genomic differences between SMN1 and SMN2 .

[0178] The SMN1/2 statuses of 12,747 short-read whole genomes sequenced to great depth (>30x) in five ethnic populations were characterized. From these samples, a total of 251 (1317) samples with complete loss of the SMN1 gene (conversion) and 6241 (374) samples with loss of the SMN2 gene (conversion) were identified. A pan-ethnic carriage rate of 2% was calculated, which is consistent with previous studies. In addition, specific copy numbers were confirmed, and all copy number determinations (48/48) of SMN1 and 98% (47/48) of SMN2 were consistent with those measured by digital PCR.

[0179] This WGS-based SMN copy number determination method can be used to identify the determination of both carrier and exposed SMA status, allowing SMA testing to be offered as a comprehensive test in neonatal diagnostics, as well as an accurate screening tool for carrier status in large-scale settings. WGS sequencing projects.

Introduction

[0180] With recent advances in next generation sequencing (NGS), it is now possible to profile large numbers of genes or even the entire genome with high throughput and in a clinically relevant time frame. Based on these advances, efforts are underway in many countries to sequence populations on a large scale, in which testing for rare genetic disorders, including carrier status, will be a major factor. Spinal muscular atrophy (SMA), an autosomal recessive neuromuscular disorder characterized by the loss of alpha motor neurons, causes severe muscle weakness and atrophy, which is atrophy that occurs during or shortly after birth. SMA is the leading genetic cause of childhood mortality after cystic fibrosis. The incidence of SMA is 1 in 6000–10000 live births, and the carrier rate is 1:40–80 among various ethnic groups. Four clinical types of SMA are classified according to age of onset and severity of the disease: very weak infants unable to sit without support (type I), weakly sitting but unable to stand (type II), outpatients with weaker legs than arms (type III) , and the presentation of SMA in adult patients is quite benign (type IV). Early detection of SMA may be critical to long-term quality of life due to the availability of two early treatments, Nusinersen and Zolgensma, which have received FDA approval for the relief of SMA.

[0181] The SMN region includes two paralogous genes: SMN1 and SMN2 . SMN2 is located at a distance of 875 kb. from SMN1 on chromosome 5q and is caused by an ancestral gene duplication unique to the human lineage. The genomic region around SMN1/2 undergoes uneven cross-linking and gene conversion, resulting in copy number variation (CV) of SMN1 and SMN2 . SMN2 shares more than 99.9% sequence identity with SMN1 , and one of the base differences, c.840C>T in exon 7, is of critical functional significance. By interrupting the splicing enhancer, c.840T promotes exon 7 skipping, resulting in the vast majority of SMN2- derived transcripts (70-85%, depending on the tissue) being unstable and not fully functional. Approximately 95% of SMA cases result from a biallelic lack of the functional nucleotide c.840C, caused by deletion of SMN1 or gene conversion to SMN2 (c.840T). In the remaining 5% of cases, patients with SMA have other pathogenic variants in SMN1 in trans configuration lacking the c.840C allele. 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.

[0182] Due to the high incidence and severity of the disease, extensive screening for SMA is recommended by the American College of Medical Genetics. The benefit of carrier screening in the population has been demonstrated in pilot studies. Screening for SMA includes: 1) determination of SMN1 copy number for SMA diagnosis and carrier testing and 2) determination of SMN2 copy number for clinical classification and prognosis. Traditionally, SMA testing and carrier testing are performed using polymerase chain reaction (PCR)-based assays such as quantitative PCR (qPCR), multiplex ligation probe amplification (MLPA), and digital PCR. These methods mainly determine the copy number of SMN1 based on the c.840C>T site, which differs between SMN1 and SMN2 . This example demonstrates that WGS can match or exceed the performance of these tests and indicates that both current and future precision medicine initiatives can use genomic data for population-level screening.

[0183] Reproducing the current SMA testing regime poses a challenge for high-throughput WGS due to the near-perfect sequence identity between SMN1 and SMN2 . In addition, frequent gene conversion between SMN1 and SMN2 is thought to result in the formation of hybrid genes. These problems require the application of a bioinformatics method to overcome the difficulties in this field. Two NGS-based tests for SMA carrier detection are presented. 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;16:100. doi:10.1186/s12881-015-0246-2) used a Bayesian hierarchical model to calculate the probability that the proportion of reads derived from SMN1 is equal to or less than 1/3 given three base differences between SMN1 and SMN2 . The method described in Larson's publication allows testing for SMA; although, since the method does not perform copy number determination, it is not an ideal solution for carrier screening. Conversely, 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; 19(8): 936-944 doi:10.1038/gim.2016.215) describes copy number determination for both SMN1 and SMN2 based on targeted sequencing data that closely mimics the current qPCR method. The Feng method is designed for targeted sequencing and therefore requires special normalization, which limits the method to a single analysis at a single site. The method calculates 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.840 C>T site. Using full coverage and the SMN1 : SMN2 ratio, the method determines the absolute copy number of SMN1 and SMN2 . Because this method is based on only one locus, it is unreliable for WGS data, where the depth variability for each locus can be very high.

[0184] Compared to targeted sequencing, WGS provides much more uniform genome coverage and provides a lower bias approach for detecting copy number variants (CNVs). In addition, WGS provides the opportunity to comprehensively profile the spectrum of population variation in the SMN region, the understanding of which at the sequence level leaves much to be desired. This example describes a new method that detects copy numbers of both SMN1 and SMN2 using WGS data. While most traditional assays only test for the absence of c.840C as a proxy for "exon 7 deletion", this example describes a method that can more fully characterize variability in a region including: 1) DNA deletions, including complete deletion/duplication gene and partial deletion of the region that includes exons 7 and 8; and 2) discovery of small variants, including the g.27134T>G SNP, that correlate with silent SMA carriers (two copies of SMN1 on the same haplotype). The accuracy of this method was demonstrated by comparing copy number determinations using digital PCR with the WGS-based determinations of the example. Concordance was shown to be 100% (48/48) for SMN1 and 98% (47/48) for SMN2 . In addition, the method was applied to 2504 unrelated samples from the 1000 Genomes project and 10243 unrelated samples from the NIHR BioResource project to report the copy number distribution of SMN1 and SMN2 . The carrier frequencies for SMA determined using the method described in the example are consistent with data obtained in previous PCR-based studies. In addition to demonstrating the accuracy of the way to quantify variants in the SMN region, this example highlights the importance of using ethnically diverse populations when developing new information methods to identify complex clinically relevant regions of the genome.

Materials and methods

Samples and data processing

[0185] Samples tested by digital PCR were obtained from the collection of the Motor Neurone Disease Research Laboratory (Nemours Alfred I. duPont Hospital for Children) and were obtained from cell lines as described above. This cohort contained 29 SMA samples (14 SMA type I, 1 SMA type I/II, 10 SMA type II, 3 SMA type III and 1 SMA of unknown clinical grade), six samples with non-SMA neuromuscular diseases ( including hereditary sensory and autonomic neuropathy 3, myotonic dystrophy type I, distal hereditary motor neuropathy type I and peripheral neuropathy type I and sensorimotor Charcot-Marie-Tooth neuropathy type IA), as well as 13 normal samples. WGS was performed using the TruSeq DNA PCR-free sample preparation kit using 150-bp paired-end reads sequenced on HiSeq X instruments from Illumina (San Diego, CA, USA). The GRCh37 genome assembly was used for read alignment.

[0186] The population-based studies included 13,343 subjects from the NIHR BioResource Rare Diseases project (EGAS00001001012) who performed WGS on people with rare diseases and their close relatives. Additional subjects (n = 840) from the Next Generation Children project (EGAD00001004357) were also studied who performed diagnostic trio WGS on patients and their parents from neonatal and pediatric intensive care units in the UK. WGS for these studies was performed using Illumina's TruSeq DNA PCR-Free Preparation kit with paired-end 100- or 125-bp reads sequenced on an Illumina HiSeq 2500 instrument or paired-end 150-bp reads sequenced on a HiSeq X instrument The GRCh37 genome assembly was used for read alignment. In the population analysis, related subjects and subjects of unknown origin were excluded, leaving 10,243 unrelated subjects.

[0187] For data from the 1000 Genomes Project (1kGP), WGS BAM was downloaded from ncbi.nlm.nih.gov/bioproject/PRJEB31736/. These BAM files were generated by sequencing 2 x 150 bp reads. on Illumina NovaSeq 6000 instruments from non-PCR libraries sequenced to an average depth of at least 30x and aligned to the human hs38 DH reference using BWA-MEM v0.7.15 (average genome coverage >30x).

Analysis of SMN copy number by orthogonal methods

[0188] For validation samples, copy numbers of SMN1 and SMN2 were measured using a QuantStudio 3 D digital PCR system (Life Technologies, Carlsbad, CA) using exon 7 allele-specific probes as described above. SMN1 and SMN2 copy numbers were normalized to RPPH1 ( RNase P ) copy numbers. The detected SMA patterns in the Next Generation Children project were confirmed using standard MLPA (SALSA MLPA P060 SMA Carrier Probemix, MRC-Holland).

Copy number determination for intact and truncated SMN

[0189] The SMN1 and SMN2 loci are affected by the copy numbers of two common variants, the copy numbers of genome-wide variants, and partial deletion of exons 7 and 8 (see results of this example). A truncated form of SMN with partial deletion of exons 7 and 8 was named SMN* . A method for determining the copy number of intact SMN1 + SMN2 genes (hereinafter referred to as SMN ) and truncated SMN genes ( SMN *) using the following steps.

[0190] Read detection and counting SMN1 And SMN2 : Read counts were calculated directly from the WGS aligned BAM file based on all reads mapped toSMN1 orSMN2, including reading with zero matching quality. Often the read data will align to these regions with zero mapping quality because the sequence of the two regions is identical. These two genes share sequence only with each other and not with other regions of the genome. The number of reads in the 22.2 kb region, including exons 1 to 6, was used to calculate the total copy numberSMN (SMN1,SMN2 AndSMN*), and the number of reads in the 6 kb region, including exon 7 and exon 8, was used to calculate the number of copies of the intactSMN (SMN1 AndSMN2).

[0191] Calculation of normalized depth of regions SMN : The read counts of the two regions described above were normalized by region length and further normalized by dividing by the average depth of 3000 preselected 2 kb regions. from the genome.

[0192] Convert normalized depth to copy number : Population normalized depth values were modeled using a univariate mixture of 11 distributions that centered around each integer copy number value representing copy number states ranging from 0 to 10. Copy number of total SMN and intact SMN was determined using a normal distribution mixture model (GMM) with a posterior probability threshold of 0.95.

[0193] Calculation of the number of copies of intact and shortened SMN : Number of copies of intactSMN was defined as the copy number of the 6.3 kb region spanning exons 7 and 8. The copy number of the truncatedSMN (SMN*) was obtained by subtracting the number of copies of the intactSMN from the total number of copiesSMN, calculated from a 22.2-kb region containing exons 1–6.

Genotyping the number of copies of alleles for individual bases

[0194] The number of chromosomes carrying the bases SMN1 and SMN2 was determined by combining the total copy number of SMN with the number of reads supporting each of the gene-specific bases. Based on determining the copy number of intact SMN at each position, the method iterated through all possible copy number combinations of SMN1 and SMN2 , outputting the combination that gave the highest posterior probability for the observed number of reads supporting SMN1 and SMN2 . In addition to determining the copy number of bases specific to SMN1 or SMN2 , this method can be used at variant positions to determine the copy number of SNPs that are known to be specific for one of the two genes, for example g27134T>G, as described below.

Number of copies of SMN1 and SMN2

[0195] For the 16 positions (localized from intron 6 to exon 8) that differ between SMN1 and SMN2 in the reference genome, it was examined whether these sites were indeed fixed in the population by comparing the specific copy numbers of the SMN1 alleles for these positions with the specific copy numbers of the splice -variant form of the base SMN1 c.840C. Eight positions, including c.840C>T, where SMN1 bases are fixed or fixed in the population, were identified based on matches to the splice variant form of the base (see Results section in this example, FIG. 14A). The remaining sites may be polymorphic in the population and may not be reliable for use in determining copy number.

[0196] To obtain a final copy number determination by this method, the following steps must be performed: 1) SMN1 copy number determinations are consistent at least 5 of 8 sites at a posterior probability cutoff of 0.8, or 2) at least 5 of 8 sites ( with a posterior probability >0.6) are consistent with the copy numbers obtained from overlapping reads from all 8 sites (with a posterior probability >0.9). Otherwise, copy number determinations were not performed for either SMN1 or SMN2 . SMA samples were identified as having zero copy of intact SMN1 , and carrier samples were identified as having one copy of intact SMN1 .

[0197] At higher copy number values, greater variability in read depth can be expected, resulting in less sensitive copy number determination (with a lower posterior probability) at individual sites and greater discrepancy between sites. As a result, the probability of missing detection was higher in samples with higher SMN1 / SMN2 copy numbers, i.e., both values greater than or equal to two (see FIG. 15). However, in such samples it is still possible to reliably determine whether the SMN1 copy number is 0 (SMA) or 1 (carrier), allowing for the determination of SMA/no-SMA or carrier/non-carrier. If SMN1 copy number was not determined, if at least seven of the SMN1 copy number determinations were significantly greater than zero, then the sample was designated as “no SMA.” Likewise, if at least seven SMN1 copy number determinations were significantly greater than one, the sample was designated as a “non-carrier.” In addition, if SMN1 copy number was not determined, the absence of the c.840C allele, which would indicate SMA, was directly tested. This was done by testing whether the number of reads supporting the SMN1 base (c.840C) was more likely in the absence of SMN1 or one copy of it.

results

Total copy number of variant affecting SMN1 / SMN2 loci

[0198] The SMN1 and SMN2 genes are located in the ~2 Mb region. in a reference genome with a large number of complex segmental and inverted segmental duplications. Although existing methods (eg, PCR-based methods) focus primarily on the c.840C>T site, this example illustrates a copy number approach based on sequencing data for complete genes. SMN1 copy number was defined as the number of SMN genes carrying the c.840C allele, and SMN2 copy number was defined as the number of SMN genes carrying the c.840T allele. Sequence analysis was performed using high depth (>30X) WGS data from 2,504 samples from the 1000 Genomes Project (1kGP), as well as 10,243 unrelated samples from the NIHR BioResource project (see methods in this example).

[0199] To formulate a copy number determination strategy, the copy numbers of two common variants that resulted in DNA deletions were first characterized. The primary copy number estimate for variants includes the entire SMN1 / SMN2 gene region. The read depth in the ~30 kb homologous region containing the SMN1 and SMN2 genes was examined. In FIG. Figure 12A shows normalized read depths in 100 bp sliding windows. in samples with different copy numbers of SMN1+SMN2 variants in a given region (representing both SMN1 and SMN2 ). The depth profile shows that the entire region has been removed or duplicated in these samples. The exact breakpoints of a given copy number of variants were expected to vary between samples due to extensive sequence homology within and outside a given region, and could only be distinguished at high resolution with long reads. To test for SMA, the analysis was limited (~30 kb) to regions including the SMN genes ( SMN1 or SMN2 ).

[0200] In addition to the copy number variants of the entire gene, a partial deletion of 6.3 kb was detected. gene spanning both exons 7 and 8 (FIG. 12B, FIG. 16). The breakpoint sequences are identical between SMN1 and SMN2 , so this deletion occurs at either chr5:70244114 - 70250420 in SMN1 or chr5:69368689 - 69375000 in SMN2 (FIG. 16, hg19). However, about 500 bp. below the breakpoint defining the end of this deletion, there are three base differences between the SMN1 and SMN2 loci (70250881A>69375425C, 70250981A>69375525G, 70250991A>69375535G). Among the samples containing this deletion, 245 read pairs were identified from 237 samples, where one read spanned the breakpoint and the other spanned at least two of the three SMN differentiating bases. Analysis of these read pairs showed that 100% corresponded to a deletion occurring in the background of the SMN2 sequence. This truncated form of SMN2 has been named SMN* , and since both exons 7 and 8 are deleted, SMN* likely has limited or no biological function. Thus, SMN* is an important variant to consider in any determination of SMN copy number.

[0201] In FIG. 12A and 12B show non-limiting example graphs illustrating common ICCs affecting the SMN1 / SMN2 loci. In FIG. 12A shows depth profiles in the SMN1/SMN2 regions. Samples with total SMN1 + SMN2 copy numbers of 2, 3, 4, and 5 are shown as dots, respectively. For each copy number category, the depth of 50 samples is summed. Each point represents normalized depth values in a 100-bp window. Read counts were calculated in each 100-bp window, the reads for SMN1 and SMN2 were summed, and normalized to the depth of the wild-type samples (CN=4). SMN exons are represented as magenta boxes. The two x-axes show the coordinates in SMN1 (bottom) and SMN2 (top). In FIG. 12B shows depth profiles pooled from 50 samples carrying deletions of exons 7 and 8, shown as dots. Read depth values were calculated in the same manner as shown in FIG. 12A.

[0202] After searching for anomalous read pairs, no other total copy number variants were found in the SMN region. Combining this information together, SMN gene copy numbers were designed to specifically identify the number of intact and truncated forms by dividing the genes into two regions: a 6.3-kb region containing exons 7-8, and a 22.2-kb region .bp containing exons 1-6. The copy number of these two regions was calculated from the read depth as described in the Methods section of this example. The copy number calculated based on the region of exons 7-8 provided the number of intact SMN genes. Samples with SMN* had higher levels of copy number detection from the exon 1-6 region compared to copy number detection from the exon 7-8 region, and their difference was represented by SMN* copy number. In FIG. Figure 13 shows the results of this calculation for 12,747 cohort samples, where 2,144 SMN* instances were identified, including 140 samples with two copies of SMN* and one sample with three copies of SMN* .

[0203] In FIG. 13 shows a non-limiting example of a scatterplot of total SMN copy number ( SMN1 + SMN2 ) (x-axis, denotes read depth of exons 1-6) and intact SMN copy number (y-axis, denotes read depth of exons 7-8).

Copy number differentiation of SMN1 from SMN2

[0204] After calculating the total copy number of the SMN genes, SMN1 and SMN2 were differentiated as described below. Since c.840C>T is the most important functional difference between SMN1 and SMN2 , the absolute copy number of these two genes can theoretically be obtained by the ratio between the number of reads supporting SMN1 and SMN2 at that site. However, read depth at a single diploid position is typically 30 - 40X for a WGS data set and sometimes does not provide sufficient power to clearly distinguish between different copy number states (see FIG. 15). Thus, when performing copy number determinations, additional base differences near c.840C>T were used so that information at these sites could be combined with c.840C>T. Because it is desirable to differentiate intact SMN1 from SMN2 , variants that occur within the 6.3-kb deletion were considered. SNPs in homopolymers and short tandem repeats (TRs), which may be more error-prone, were excluded, resulting in 16 base differences between SMN1 and SMN2 (Table 8).

[0205] For these 16 base differences, the copy number of the SMN1 and SMN2 alleles was determined independently (see Methods section of this example) and the copy number determinations for each position were compared with the copy number determinations at the splice variant site (FIG. 14A, FIG. 17). There was a noticeable difference between the agreement of definitions in African and non-African populations (FIG. 14A). For non-African samples, 13 sites with high (>85%) copy number consistent with the splice site were found. Conversely, for African samples, only seven splice site-matched high copy number regions were found, and the match values were lower than in non-African populations. This is consistent with intragenic variation at many of these positions and higher frequencies for these non-reference alleles in non-African populations. The splice variant and seven positions that were highly consistent with the splice variant in both African and non-African populations were selected to determine copy numbers of SMN1 and SMN2 . By limiting ourselves to two copy number states that allow easy identification of hybrid alleles ( SMN1 =CN2 and SMN2 =CN0 or SMN1 =CN2 and SMN2 =CN1), it became possible to estimate the allele frequencies of these sites in the SMN1 and SMN2 genes (Table 9, FIG. 18A and 18B). Based on this analysis across these eight positions, it was estimated that up to 0.5% of SMN1 genes contain the SMN2 allele. Conversely, it is estimated that up to 0.9% of SMN2 genes carry the SMN1 allele. These observations may be the result of gene conversion or the fact that many of these eight sites are polymorphic in the population. Most of these hybrid alleles originate from African populations (Table 9).

[0206] In FIG. 14A-14D show the copy number distributions of SMN1 / SMN2 / SMN * in the population. In FIG. 14A is a non-limiting example of an illustrative graph illustrating the percentage of samples showing copy number assignment agreement with c.840C>T across 16 SMN1-SMN2 base difference sites in African and non-African populations. Site 13* is a splice variant site c.840C>T. The black horizontal line indicates 85% agreement. In FIG. 14B shows non-limiting examples of histograms of copy number distributions of SMN1 , SMN2 and SMN* across five populations in the 1kGP and NIHR BioResource cohort (numbers given in Table 15). In FIG. 14C shows a non-limiting example of a plot of SMN1 copy number versus total SMN2 copy number (intact SMN2 + SMN* ). In FIG. Figure 14D shows two trios with an SMA proband identified by a specialist and orthogonally confirmed in the NIHR BioResource cohort. The copy numbers per allele of SMN1 , SMN2 , and SMN* are phased and labeled for each member of the trio.

[0207] Introducing more base differences increased the ability to differentiate SMN1 from SMN2 . However, since these sites are not truly invariant in the corresponding genes, and copy number determinations at individual sites may be subject to error, the likelihood that one of the individual determinations will differ from the true copy number state will increase. To make the final determination, it was necessary that the SMN1 copy number determinations were consistent with each other at 5 or more of the 8 sites (see the Methods section of this example for a complete description of the copy number determination rules). With a posterior probability cutoff of 0.8, most samples had concordant definitions in at least five of the eight sites, and only 1.4% of samples had fewer than 5 concordant sites (Table 10). In 80% of these samples, a reliable copy number determination was made based on the second consensus rule (requiring agreement with the copy number determination made by summing the reads from all 8 sites). "Discordant" sites were more often not identified due to low posterior probability rather than due to inconsistent definitions, and only 15.3% of these were valid definitions that did not agree with the consensus of other sites. Again, a significant proportion of the discrepancies were identified in African populations (Table 10). Using fewer sites for most rules resulted in more misidentifications and misidentifications compared to using eight sites (Table 11).

Checking the SMN Copy Number Determination

[0208] To test this method, 48 samples with known SMN1 and SMN2 copy numbers were sequenced, including 29 SMA probands, 6 SMA carriers, and 13 samples with SMN1 copy numbers greater than 1. SMN1 copy number results were consistent with digital PCR results for all 48 cases , and SMN2 copy number results were consistent for 47 (97.9%) of 48 cases (Tables 6A and 6B). In this single nonmatching case (MB509), the method detected 3 copies of SMN2 , while digital PCR showed 2 copies of SMN2 (Table 12). A more precise study revealed a deletion of 1884 bp. in SMN1 (chr5:70247145-70249029, hg19) in this sample (FIG. 19). The deletion is small (does not significantly change the depth in the 6 kb region used to determine the copy number of intact SMN ) and has not previously been reported (or detected in population data), so the method was not designed to detect it. As a result, this sample was correctly identified as SMA, but the copy number of SMN2 was overestimated by one. The deletion is consistent with copy number determinations made at 8 sites of the SMN1 - SMN2 difference, with the first 2 sites not in the deletion and SMN1 copy number being defined as CN=1, and the next 6 sites being in the deletion and SMN1 copy number being defined as CN=0 .

[0209] The consistency of SMN1/SMN2 / SMN* copy number determinations in 258 trios from the Next Generation Children cohort was analyzed (see Methods section of this example). There was no Mendelian error in any of the determinations (Table 13).

Table 6A. Validation against samples with known SMN1 / SMN2 copy numbers

Number of copies by digital PCR General Agreed Inconsistent Coordination SMN1 0 29 29 0 100.0% 1 6 6 0 100.0% 2 10 10 0 100.0% 3 3 3 0 100.0% General 48 48 0 100.0% SMN2 0 1 1 0 100.0% 1 4 4 0 100.0% 2 29 28 1 96.6% 3 eleven eleven 0 100.0% 4 3 3 0 100.0% General 48 47 1 97.9%

Table 6B. Validation against samples with known SMN1 / SMN2 copy numbers

Number of copies using the orthogonal method General Agreed Inconsistent Coordination SMN1 0 64 64 0 100.0% 1 45 44 1 97.8% 2 897 897 0 100.0% 3 174 174 0 100.0% 4 43 43 0 100.0% 6 1 0 1 0.0% General 1224 1222 2 99.8% SMN2 0 117 117 0 100.0% 1 466 465 1 99.8% 2 541 539 2 99.6% 3 60 60 0 100.0% 4 9 8 1 88.9% General 1193 1189 4 99.7% SMN2Δ7 - 8 0 1089 1089 0 100.0% 1 80 80 0 100.0% 2 4 4 0 100.0% General 1173 1173 0 100.0%

Copy numbers of SMN1 , SMN2 and SMN* by population

[0210] Given the high accuracy demonstrated by validation of digital PCR results, the method was applied to high depth (>30X) WGS data for 12,747 unrelated samples from 1kGP and NIHR BioResource (Table 14). Copy number distribution was analyzed across populations (Europeans, Africans, East Asians, South Asians, and mixed Americans consisting of Colombians, Mexicans, Peruvians, and Puerto Ricans). In FIG. 14B shows a histogram of the number of individuals with different copy numbers of intact SMN1 , intact SMN2, and SMN* . The distributions are similar between the 1kGP and NIHR BioResource samples (FIG. 20). Overall, individuals had more copies of SMN1 than SMN2 . The most common SMN1/SMN2 copy number combinations were 2/2 (44.9%) and 2/1 (33.4%). With the exception of Africans, who showed higher variation in both SMN1 and SMN2 copy numbers, variation in SMN1 copy numbers was much lower than in SMN2 copy numbers. Conversely, 54.7% of Africans had three or more copies of SMN1 , which was more than twice as high as observed in any of the other four populations (FIG. 14B, Table 7). There was an inverse relationship between SMN1 and SMN2 copy number, SMN2 copy number was lower with increasing SMN1 copy number (FIG. 14C, correlation coefficient -0.344, p-value < 2.2e-16). This observation is consistent with a mechanism in which gene conversion occurs between SMN1 and SMN2 . The observed higher copy number of SMN1 compared to that of SMN2 may result from a bias towards conversion of SMN2 to SMN1 or a selection towards low copy number of SMN1 . Africans have significantly lower SMN2 copy numbers than other populations.

[0211] The number of SMA carriers identified in different populations is summarized in Table 7 and Table 15. Of the 12,683 individuals with valid SMN1/SMN2 copy number determinations, Europeans had the highest carrier frequency at 2.2%, followed by Mixed Americans ( 2.05%), East Asians (1.35%), and South Asians (1.67%). Africans had the lowest carriage rate (0.44%). The copy number frequency distributions observed in this example are consistent with previous studies of SMN1/SMN2 copy number distributions in the general population. In addition, the frequency of deletion in exon 7-8 ( SMN* ) was determined in different populations: 21.2% of Europeans and 11.5% of mixed Americans had at least one copy of SMN* , while the frequency was lower in people from the South Asian (3.35%), African (1.1%) and East Asian (0.34%).

[0212] In the Next Generation Children cohort (see Methods section of this example), SMA was identified in two neonatal probands from three analyzes that were independently confirmed. In addition, the copy numbers of SMN1 , SMN2 and SMN* were phase-regulated for each member of the trio (FIG. 14D).

Simulation of determining the number of copies of one site

[0213] The number of reads per site at median sampling depths of 30X, 35X and 40X were modeled based on a Poisson distribution, and reads supporting SMN1 were selected based on a binomial model with all possible combinations of SMN1 and SMN2 copy numbers, with the total SMN copy number being from 2 to 6. Taking into account the number of supporting reads of SMN1 and SMN2 , the posterior probability of the simulated copy number of SMN1 was obtained (see Methods section in this example). The posterior probability was high (greater than 0.9) when at least one copy number value of SMN1 or SMN2 was low (less than or equal to 1) (FIG. 16). When both values were greater than 2, i.e. in the SMN1:SMN2 combinations 2:2, 2:3, 2:4, 3:2, 3:3 and 4:2, the posterior probability often became low and dropped below 0. 9. This is due to greater variability in read depth when the expected copy number is higher. Thus, in these scenarios, generating SMN1 and SMN2 copy number determinations using a single site may be less accurate.

Discrepancies in test samples

[0214] There was one sample of MB509 that showed a discrepancy between our copy number determination and digital PCR results. Upon further analysis, this sample had two copies of SMN2 and one copy of SMN1 with a deletion of 1884 bp. (chr5:70247145-70249029, hg19, FIG. 20). Although read alignment in the SMN1 /2 region is not always precise, careful analysis of the split read results revealed that the reads or their partners overlapped SMN1 -specific bases. Without being limited to any theory, it has been hypothesized that this deletion was correctly placed on SMN1 . The deletion is small (does not significantly change the depth in the 6.3 kb region used to determine copy number of intact SMN) and has not previously been reported (nor found in 1kGP samples, which is a very rare variant), so the method was not designed to detect deletion. As a result, the method determined the total copy number of SMN1 + SMN2 to be 3. The deletion is consistent with the copy number determinations made at the 8 SMN1 - SMN2 difference sites, with the first 2 sites not in the deletion and the SMN1 copy number being determined to be CN=1, and the following 6 sites were in the deletion and the copy number of SMN1 was determined to be CN=0 (FIG. 21A). Based on most rules, the method determined the copy number of SMN1 to be 0, correctly identifying the sample as SMA. The SMN2 copy number was calculated as the total copy number minus the SMN1 copy number, so the method determined the SMN2 copy number to be 3, overestimating it by 1.

[0215] Four other samples, MB231, MB367, MB383 and LP2101748, had discrepancies between the copy number determinations performed and the digital PCR or MLPA results. Read counts and normalized depth values (read counts divided by the depth of haploid samples) at the 8 base difference sites supported our copy number determinations (FIG. 21A), and the discrepancy was likely caused by errors in the orthogonal methods. In two samples, genomic sequencing (GS) detection and digital PCR detection differed by a factor of two (MB231: GS-0.2, PCR-0.4 and MB383: GS-3.1, PCR-6.2). When using digital PCR, a normalization problem may arise, leading to an overestimation of the copy number by a factor of two.

[0216] When comparing MLPA copy number determinations to 1109 1kGP samples, one sample was excluded that did not detect SMN2Δ 7-8 due to the low posterior probability for the total SMN copy number, as well as three samples where the copy number of SMN1 and SMN2 was not determined due to discrepancies in copy number determinations across 8 selected sites that did not follow consensus rules (FIG. 22B).

Detection of silent media

[0217] g.27134T>G SNP may be associated with 2+0 SMA silent carrier status, where one chromosome carries two copies of SMN1 (either by duplicating SMN1 or by gene conversion of SMN2 to SMN1 ) and the other chromosome has no copies of SMN1 . The method of this example can also detect the presence of this SNP and thus can be used to screen for potential silent carriers. This SNP is most closely associated with two-copy SMN1 alleles in Africans, where 84.5% of subjects with three copies of SMN1 and 92.6% of subjects with four copies of SMN1 have the g.27134T>G SNP (Table 7). Identification of this SNP has significantly increased the frequency of carrier detection in Africans, as Africans have a higher frequency of alleles carrying two copies of SMN1 (Table 17 and Table 18). However, 33% of subjects with two copies of SMN1 also had the g.27134T>G SNP, indicating that a significant proportion of SMN1 singlet alleles also carry this SNP. Maximum likelihood estimates were calculated for the percentages of single- and double-copy SMN1 alleles carrying g.27134T>G (Table 17) and residual risks for the combination of copy number and SNP determinations (Table 18). The calculated estimates are similar to previous studies, although there are significant differences in all of these estimates. This variability may be due to population variation: for example, Africans (this example) versus African Americans (previous studies) and Northern Europeans (overrepresented in this example) versus a more diverse sample of Caucasians (previous studies).

Table 7. SMN1 copy number and g.27134T>G frequency depending on the population

Ethnicity General SMN1CN=1 SMN1CN=2 SMN1CN=3 SMN1 CN=4 Quantity g.27134T
>G+ Quantity g.27134T
>G+ Quantity g.27134T
>G+ Quantity g.27134T
>G+ African population 902 4 0 (0.0%) 404 134 (33.17%) 373 315 (84.45%) 121 112 (92.56%) European population 9648 212 0 (0.0%) 8899 4 (0.04%) 524 22 (4.2%) 13 2 (15.38%) South Asians 1199 20 0 (0.0%) 965 1 (0.1%) 195 5 (2.56%) 19 1 (5.26%) East Asians 593 8 0 (0.0%) 552 1 (0.18%) 33 1 (3.03%) 0 0 (no data) Mixed Americans 341 7 0 (0.0%) 296 7 (2.36%) 36 9 (25.0%) 2 1 (50.0%)

Comparison between two converters, BWA and Isaac

[0218] The method in this example over-analyzed both SMN1 and SMN2 reads, and was therefore insensitive to how the resolver distinguishes between the two genes. Thus, using different converters should produce similar results. The BAM data analyzed in this example was obtained using two different converters: BWA for the 1kGP data and various versions of Isaac for the rest. Consistent SMN1 /2 copy number distributions between 1kGP and NIHR samples (Table 19, FIG. 20) indicate that our method is insensitive to the transformer. In addition, the method was tested for consistency by aligning 117 samples with both BWA and Isaac, including 5 SMA samples and 3 media. All 117 samples had exactly the same copy number determinations ( SMN1 / SMN2 / SMN2Δ 7-8) according to the method in this example, and the normalized depths for both exons 1-6 and exons 7-8 were virtually identical (r > 0.999 , FIGURE 22).

Comparison between carrier definitions in this study and Larson et al.

[0219] The carrier determinations obtained from the 1kGP samples in this example (N=37) were compared to those reported by Larson et al. (N=36) and found 26 overlapping definitions (Table 15). It is assumed that the determinations obtained by the method of this example are correct, Larson et al. made 10 false positive (FP) and 11 false negative (FN) determinations. Larson et al. identified carriers by determining whether the proportion of reads supporting SMN1 was less than or equal to 1/3. This study used shallow-depth sequencing data, which is expected to introduce some errors, but more importantly, their approach is prone to errors without determining the total copy number. For example, a sample with one copy of SMN1 and one copy of SMN2 will be determined to be a non-carrier ( SMN1 1/2 fraction), and a sample with two copies of SMN1 and four copies of SMN2 will be determined to be a carrier ( SMN1 1/3 fraction), resulting in false positives and false negatives. results (Table 16).

Additional figures and tables

[0220] In FIG. Figure 15 shows non-limiting example plots, each illustrating the posterior probability distribution for modeling SMN1 copy number using a single site at different read depths and SMN1:SMN2 copy number combinations.

[0221] In FIG. 16 shows a non-limiting example of an IGV snapshot of the SMN2 region in a sample with a deletion in exon 7-8. Horizontal lines pair two reads on the central alignment track. BLAT results for two split reads spanning the breakpoint are shown in the bottom track, showing two segments of the same alignment for the read on each side of the deletion breakpoint.

[0222] In FIG. 17 shows non-limiting examples of graphs illustrating the correlation between raw SMN1 copy numbers at 15 base differences near c840.C>T and raw SMN1 copy numbers at the c840.C>T site. The raw copy number of SMN1 at each site was calculated as the copy number of intact SMN multiplied by the fraction of SMN1 supporting the number of SMN1 + SMN2 reads from the number of reads. Correlation coefficients are given in the header of each graph.

[0223] In FIG. 18 and 18B show non-limiting examples of SMN1 / SMN2 haplotype plots in samples with SMN1 :2 SMN2 :0 and SMN1 :2 SMN2 :1 in 1kGP. The y-axis shows the raw copy numbers of SMN1 as shown in FIG. 16. The X-axis shows the 16 sites, the indexes of which are listed and explained in Table 8. Index No. 13 represents the site c840.C>T. Samples with SMN1: 2 SMN2 :0 are shown together in the top left graph. Samples with SMN1 :2 SMN2 :1 are shown as 5 clusters. FIG. 18A Non-African Population FIG. 18B African population

[0224] In FIG. 19 shows a non-limiting example of an IGV snapshot showing a 1.9 kb deletion. in SMN1 in MB509.

[0225] In FIG. 20 is a non-limiting illustrative graph illustrating SMN1/SMN2 /SMN* copy numbers in the 1kGP and NIHR cohorts.

[0226] In FIG. 21A and 21B show discrepancies and missing determinations in the test samples. In fig. Figure 21 A shows five samples with discrepancies between GS and digital PCR or MLPA determinations. The x-axis shows 16 sites, the indices of which are listed and explained in Table 8. Index #13 represents site c840.C>T. The left y-axis for the bars shows the number of reads supporting SMN1 and SMN2 . The right y-axis for lines shows normalized read depth, a proxy for SMN1 and SMN2 copy numbers (number of reads divided by haploid depth). The header of each panel shows the GS and digital PCR/MLPA determinations for each sample for SMN1 and SMN2 , separated by a comma. In FIG. 21B shows three 1kGP validation samples where the SMN detector did not detect the copy number of SMN1 and SMN2 due to discrepancy between the SMN1 / SMN2 base difference sites. For the consensus rules of the method, eight sites are used: No. 7-8 and No. 10-15. The y-axis shows the raw copy numbers of SMN1 as shown in FIG. 17.

[0227] In FIG. 22 shows copy number determinations obtained from BWA and Isaac BAM.

Table 8. Genomic coordinates of base differences between SMN1 and SMN2

Index Location Selected SMN1 SMN2 Position, hg19 Grounds Position, hg19 Grounds 1 Intron 6 70244142 A 69368717 G 2 Intron 6 70245876 T 69370451 C 3 Intron 6 70246016 G 69370591 A 4 Intron 6 70246019 T 69370594 C 5 Intron 6 70246156 G 69370731 A 6 Intron 6 70246167 T 69370742 C 7 Intron 6 Yes 70246320 G 69370895 A 8 Intron 6 Yes 70246793 G 69371368 A 9 Intron 6 70246919 A 69371499 C 10 Intron 6 Yes 70247219 G 69371799 A eleven Intron 6 Yes 70247290 T 69371870 C 12 Intron 6 Yes 70247724 G 69372304 A 13 Exon 7
(c.840 C>T) Yes 70247773 C 69372353 T 14 Intron 7 Yes 70247921 A 69372501 G 15 Intron 7 Yes 70248036 A 69372616 G 16 Exon 8 70248501 G 69373081 A

Table 9. Frequencies of SMN1 haplotypes with the SMN2 allele and SMN2 haplotypes with the SMN1 allele in two simple copy number states ( SMN1 =CN2 and SMN2 =CN0 or SMN1 =CN2 and SMN2 =CN1)

(Numbers in parentheses indicate the contribution of representatives of the African population.)

Site index Number of SMN1 haplotypes with reliable copy number determination Number of SMN1 haplotypes with SMN2 allele Percentage Number of SMN2 haplotypes with reliable copy number determination Number of SMN2 haplotypes with SMN1 allele Percentage 1 12292 490 (71) 4 5041 101 (34) 2 2 9372 542 (79) 5.8 3669 46 (0) 1.3 3 11784 187 (48) 1.6 4788 48 (1) 1 4 11056 205 (51) 1.9 4428 43 (1) 1 5 10212 312 (51) 3.1 4087 34 (1) 0.8 6 9974 1787 (111) 17.9 3946 28 (1) 0.7 7 11956 58 (0) 0.5 4874 45 (3) 0.9 8 12218 15 (1) 0.1 5005 8 (0) 0.2 9 11872 79 (47) 0.7 4831 56 (35) 1.2 10 12484 20) 0 5137 39 (29) 0.8 eleven 11964 19 (5) 0.2 4880 10) 0 12 12506 eleven) 0 5148 0 (0) 0 13 12836 0 (0) 0 5313 0 (0) 0 14 12386 9 (6) 0.1 5088 0 (0) 0 15 12544 9 (4) 0.1 5167 33 (24) 0.6 16 12336 12 (3) 0.1 5063 76 (41) 1.5

Table 10. Number of samples with different numbers of consensus sites at 8 SNP sites

(Numbers in parentheses indicate the contribution of representatives of the African population.)

SNP mapping SMN1CN=1 CN=2 CN=3 CN=4 CN=no definition General Proportion of discordant sites 8 163 6325 594 111 0 7193 (475) 0 (0) 7 52 3141 285 28 0 3506 (199) 11.3 (1.6) 6 25 1197 150 9 0 1381 (137) 16.3 (6) 5 9 356 86 6 1 458 (74) 21.1 (10) <5 2* 92* 44* 1* 36 175 (26) 19.6 (6.9)

*Determinations are made in these samples based on the second majority rule (see "Methods").

Table 11. Number of undeterminations due to nonconformity and nonconformity determinations made with a reduced number of sites.

Number of sites for majority rule 8
(Requires 5 to match) 6
(4) 4
(3) 2
(2) 1 (c.840C)
(1) Number of uncertainties due to inconsistency 175 298 766 1149 700 Number of definitions different from those obtained using 8 sites 0 0 1 6 41

Table 12. Test samples

Sample ID SMN Copy Number Determinant Digital PCR intact
SMN1 copy number intact
SMN2 copy number SMN1 copy number SMN2 copy number NA03813 0 3 0 3 NA09677 0 3 0 3 NA23689 0 3 0 3 NA00232 0 2 0 2 NA10684 0 2 0 2 NA23687 1 2 1 2 NA23688 1 2 1 2 NA03815 1 1 1 1 MB122 2 0 2 0 MB226 2 1 2 1 MB119 3 1 3 1 MB370 3 1 3 1 MB489 0 2 0 2 MB364 0 2 0 2 MB691 0 2 0 2 MB488 0 2 0 2 MB219 0 2 0 2 MB228 0 2 0 2 MB501 0 2 0 2 MB362 0 2 0 2 MB692 0 2 0 2 MB234 0 2 0 2 MB693 0 2 0 2 MB510 0 2 0 2 MB114 0 2 0 2 MB116 1 2 1 2 MB115 1 2 1 2 MB104 2 2 2 2 MB384 2 2 2 2 MB338 2 2 2 2 MB344 2 2 2 2 MB345 2 2 2 2 MB349 2 2 2 2 MB113 2 2 2 2 MB366 2 2 2 2 MB351 3 2 3 2 MB355 0 3 0 3 MB361 0 3 0 3 MB378 0 3 0 3 MB232 0 3 0 3 MB106 0 3 0 3 MB222 0 3 0 3 MB509 0 3 0 2 MB112 0 3 0 3 MB339 1 3 1 3 MB377 0 4 0 4 MB356 0 4 0 4 MB503 0 4 0 4

Table 13. Copy numbers of SMN1 , SMN2 and SMN* were determined in 258 trios in the Next Generation Children cohort

SMN1 SMN2 SMN* Number of families Father Mother Proband 1 Proband 2 Number of families Father Mother Proband 1 Proband 2 Number of families Father Mother Proband 1 Proband 2 207 2 2 2 53 2 2 2 174 0 0 0 8 2 2 2 2 29 2 1 1 20 0 1 0 8 2 3 3 27 1 2 2 15 0 1 1 8 3 2 2 23 1 2 1 15 1 0 0 7 3 2 3 23 2 1 2 9 1 0 1 4 2 3 2 17 1 1 1 6 0 0 0 0 3 1 2 1 12 2 0 1 4 1 1 1 3 1 2 2 eleven 1 1 2 3 1 1 0 2 1 1 0 9 1 1 0 2 0 2 1 2 2 2 1 7 0 1 1 2 1 0 1 0 2 2 3 2 3 6 0 2 1 2 1 0 1 1 2 3 3 3 4 1 0 1 2 2 1 1 1 2 1 1 3 0 0 0 1 0 2 2 1 2 2 3 3 2 2 1 1 1 1 2 2 1 2 1 1 1 2 0 1 2 1 2 2 2 1 3 0 2 2 1 3 1 2 2 1 1 2 2 2 1 3 2 2 2 1 3 2 2 2 2 2 2 2 2 3 2 2 3 3 2 3 2 3 1 0 1 0 1 1 0 0 1 1 3 2 1 1 4 3 1 2 3 2 1 2 4 4 1 3 0 1 1 3 1 2 1 3 2 2 1 3 2 4 1 4 1 2

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

Ethnicity 1kGP NIHR BioResource (including unrelated NGCs) NIHR BioResource (including NGC, total) African population 661 253 295 European population 503 9186 11652 South Asians 489 713 1012 East Asians 504 91 97 Mixed Americans 347 0 0 Other 0 0 1127 General 2504 10243 14183

Table 15. Copy number frequencies of SMN1, SMN2 and SMN* by population

Ethnicity General SMN1 SMN2 SMN* 1 2 3 4 0 1 2 3 4 0 1 2 African population 902 4
(0.44%) 404
(44.79%) 373
(41.35%) 121
(13.41%) 226
(25.06%) 449
(49.78%) 214
(23.73%) 13
(1.44%) 0
(0.0%) 892
(98.89%) 9
(1.0%) 1
(0.11%) European population 9648 212
(2.2%) 8899
(92.24%) 524
(5.43%) 13
(0.13%) 833
(8.63%) 3850
(39.9%) 4667
(48.37%) 279
(2.89%) 19
(0.2%) 7591
(78.74%) 1912
(19.83%) 137
(1.42%) South Asians 1199 20
(1.67%) 965
(80.48%) 195
(16.26%) 19
(1.58%) 78
(6.51%) 400
(33.39%) 686
(57.26%) 29
(2.42%) 5
(0.42%) 1155
(96.65%) 40
(3.35%) 0
(0.0%) East Asians 593 8
(1.35%) 552
(93.09%) 33
(5.56%) 0
(0.0%) 28
(4.72%) 211
(35.58%) 340
(57.34%) 12
(2.02%) 2
(0.34%) 591
(99.66%) 2
(0.34%) 0
(0.0%) Mixed Americans 341 7
(2.05%) 296
(86.8%) 36
(10.56%) 2
(0.59%) thirty
(8.8%) 136
(39.88%) 162
(47.51%) eleven
(3.23%) 2
(0.59%) 302
(88.56%) 37
(10.85%) 2
(0.59%)

Table 16. Comparison of carrier determinations obtained in 1kGP samples in this example and in Larson et al.

Sample ID Ethnicity SMN1 copy number SMN2 copy number Number of SMN copies* Defined as a carrier by Larson et al. Carrier probability adjusted by Larson et al. MLPA Verified GS Definition HG03583 AFR 1 1 0 Yes 0.645 Yes HG01205 AMR 1 1 0 Yes 0.756 HG01892 AMR 1 1 0 Yes 0.902 Yes HG01801 EAS 1 1 0 Yes 0.541 NA11932 EUR 1 1 0 Yes 0.716 NA20760 EUR 1 1 0 Yes 0.638 Yes NA20896 SAS 1 1 0 Yes 0.514 Yes HG01948 AMR 1 2 0 Yes 0.678 Yes HG02265 AMR 1 2 0 Yes 0.982 HG01085 AMR 1 2 0 Yes 1 NA20812 EUR 1 2 0 Yes 0.999 Yes NA20764 EUR 1 2 0 Yes 0.982 Yes HG00324 EUR 1 2 0 Yes 0.997 Yes NA12383 EUR 1 2 0 Yes 1 HG03953 SAS 1 2 0 Yes 0.972 HG02771 AFR 1 3 0 Yes 0.997 HG01893 AMR 1 3 0 Yes 1 HG02079 EAS 1 3 0 Yes 0.976 NA20814 EUR 1 3 0 Yes 1 HG00281 EUR 1 3 0 Yes 1 Yes HG00346 EUR 1 3 0 Yes 1 Yes HG03740 SAS 1 3 0 Yes 0.874 HG02087 EAS 1 4 0 Yes 1 HG02134 EAS 1 4 0 Yes 1 NA12778 EUR 1 4 0 Yes 1 HG01773 EUR 1 4 0 Yes 1 Yes HG01492 AMR 2 2 0 Yes 0.914 NA19723 AMR 2 2 0 Yes 0.681 NA18542 EAS 2 2 0 Yes 0.633 HG00525 EAS 2 2 0 Yes 0.763 Yes NA20792 EUR 2 2 0 Yes 0.671 Yes NA11843 EUR 2 2 0 Yes 0.509 NA19711 AFR 2 3 0 Yes 0.943 NA19346 AFR 2 3 0 Yes 0.52 Yes HG01248 AMR 2 4 0 Yes 0.935 HG01094 AMR 2 4 0 Yes 0.738 HG02156 EAS 1 0 0 No 2.36E-33 HG02180 EAS 1 1 0 No 7.26E-05 NA20790 EUR 1 1 0 No 0.489 Yes NA20787 EUR 1 1 1 No 0.322 Yes HG01686 EUR 1 1 1 No 0.00119 Yes NA19456 AFR 1 2 0 No 0.278 HG01455 AMR 1 2 0 No 0.176 HG01863 EAS 1 2 0 No 0.42 HG01612 EUR 1 2 0 No 1.20E-07 Yes NA20845 SAS 1 2 0 No 0.398 HG03928 SAS 1 2 0 No 0.442 Yes

Table 17. Maximum likelihood estimates for the percentage of single- and double-copy SMN1 alleles carrying g.27134T>G

Ethnicity Single copy allele SMN1 Double copy allele SMN1 African population 18.4% 78.5% European population 0.02%
* (1kGP Europeans: 0.11%) 4.35%
* (1kGP Europeans: 10.0%) South Asians 0.05% 2.54% East Asians 0.09% 2.94% Mixed Americans 1.2% 24.5%

*The NIHR BioResource cohort, which occupies the majority of the European population analyzed in this example due to its large sample size, includes samples from Northern Europe that carry a lower frequency of the g.27134T>G SNP than the more diverse European samples from the 1000 Genomes Project .

Table 18. SMA carrier detection and residual risk estimates from this example

Ethnicity Carriage frequency ^a Detection rate (copy number) ^a Residual risk (CN=2) Detection rate (CN+SNP) This example Luo et al. ^c Feng et al. ^d Alias et al. ^e Residual risk (CN=2, SNP -) Residual risk (CN=2, SNP+) Residual risk (CN=2, SNP -) Residual risk (CN=2, SNP+) Residual risk (CN=2, SNP -) Residual risk (CN=2, SNP+) Residual risk (CN=2, SNP -) Residual risk (CN=2, SNP+) African population 1 of 72 70.5% 1 of 129 91.8% 1 of 346 1 of 58 1 of 396
(African American) 1 of 34 1 of 375
(African American) 1 of 39 no data no data European population 1 of 47 94.8% 1 of 790 95.0% 1 of 814
(1kGP Europeans
1 of 846 1 of 12
(1kGP Europeans 1 of 27) 1 of 770 1 of 29 1 of 921 1 of 69 1 of 888
(Spaniards) ~1 Asian ^b 1 of 59 93.3% 1 of 767 93.4% 1 of 779 1 of 57 1 of 702 ~1 1 of 907 1 of 61 no data no data Mixed Americans 1 of 68 90.0% 1 of 559 91.9% 1 of 674 1 of 71 1 of 1762
(Spaniards) 1 of 140 1 of 906
(Spaniards) 1 of 99 no data no data

a Numbers and frequencies of SMN1 alleles for residual risk calculations, taken from Sugarman et al. Panethnic carrier screening and prenatal diagnosis of spinal muscular atrophy: clinical laboratory analysis of >72,400 specimens. Eur J Hum Genet. 2012;20(1):27-32. doi:10.1038/ejhg.2011.134;

b Includes people from South and East Asia;

with Luo et al. An Ashkenazi Jewish SMN1 haplotype specific to duplication alleles improves pan-ethnic carrier screening for spinal muscular atrophy. Genet Med Off J Am Coll Med Genet. 2014;16(2):149-156. doi:10.1038/gim.2013.84;

d 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;19(8):936-944. doi:10.1038/gim.2016.215;

e Alias et al. Utility of two SMN1 variants to improve spinal muscular atrophy carrier diagnosis and genetic counseling. Eur J Hum Genet. 2018;26(10):1554. doi:10.1038/s41431-018-0193-4.

Table 19. SMN1/SMN2/SMN2 Δ 7-8 copy numbers in 1kGP and NIHR cohorts

1kGP NIHR Kolmogorov-Smirnov criterion General SMN CN=1 2 3 4 General SMN1CN=1 2 3 4 EUR 503 15 463 25 0 9145 197 8436 499 13 1 2.98% 92.05% 4.97% 0.00% 2.15% 92.25% 5.46% 0.14% EAS 502 7 470 25 0 91 1 82 8 0 0.9999 1.39% 93.63% 4.98% 0.00% 1.10% 90.11% 8.79% 0.00% AFR 653 3 293 261 96 249 1 111 112 25 0.8284 0.46% 44.87% 39.97% 14.70% 0.40% 44.58% 44.98% 10.04% SAS 489 5 397 77 10 710 15 568 118 9 1 1.02% 81.19% 15.75% 2.04% 2.11% 80.00% 16.62% 1.27%

Discussion

[0228] Due to the high sequence homology between SMN1 and SMN2 , the SMN region is difficult to determine by both short and long read sequencing, and until now this important region has been excluded from standard WGS analysis. This example demonstrates a method that can independently determine the copy number of SMN1 and SMN2 using short-read WGS data, filling an important gap in SMA diagnosis and carrier screening for precision medicine initiatives. Accurate measurement of SMN1 and SMN2 copy numbers is not only important for the diagnosis of SMA, but is also a prognostic indicator and the basis of therapeutic options. SMN2 copy number has been used as a criterion in many clinical trials of SMA, including Nusinersen and Zolgensma.

[0229] To demonstrate this method, the copy number of SMN1 and SMN2 was determined using sequencing data from 12,747 samples covering five different subpopulations. The following were identified: 251 samples with loss of the entire gene (less than two copies) and 1317 with gain of the entire gene (more than two copies) SMN1 ; 6241 samples with loss of the entire gene and 1274 with gain of the entire SMN2 gene; 2144 samples carrying one or more copies of the truncated form SMN* . It is not possible to determine precisely the role played by deletions, duplications, or gene conversion leading to copy number changes in a given region. Evidence supporting all three mechanisms includes: 1) 3853 samples with a total copy number of ( SMN1 + SMN2 ) CN<4 (deletions), 2) 670 samples with a total CN>4 (duplicates), and 3) a strong inverse correlation between the number copies of SMN1 and SMN2 (gene conversion, FIG. 14C). In addition, a carrier frequency of 1:42 to 1:101 was determined, depending on the ancestral population (Table 7). Copy number frequencies varied widely across the population, and the population results in this example are consistent with previous population-based studies. Although the alignment provides qualitative support for the accuracy of the method, the accuracy of the method was directly assessed by comparing copy number determinations obtained by the method with digital PCR results. In this direct comparison, all copy number determinations (48/48) of SMN1 and 98% (47/48) of SMN2 were consistent with the digital PCR-based results. One of the discrepancies was due to a 2 kb deletion that was not targeted by the method, and importantly, the method correctly identified the SMA status of this sample.

[0230] In this example, copy number determination was optimized for individuals of all backgrounds and thus limited differentiation of SMN1/2 to a functionally significant splice variant and seven high-match splice variant sites in all populations (FIG. 14A). By quantifying the correspondence between all reference differences and the splice variant, the method allowed us to identify variations in these fixed differences that, if not adequately accounted for (e.g., removed from our analysis), could lead to errors in our copy number determinations. The lack of fixed differences would be particularly problematic when analyzing Africans, since they have more diverse haplotypes. Population genetic studies, for example, including those using long-length sequencing, can help more directly profile haplotype diversity in populations and identify new site variants that can further improve the accuracy of SMN1 / SMN2 differentiation.

[0231] One type of silent carrier occurs when a subject has two copies of the SMN1 gene, but they both belong to the same haplotype. The SNP (g.24134T>G) was used to identify subjects who are at increased risk of becoming carriers if the SMN1 copy number value is two, but the risk associated with this SNP can vary widely between studies and populations (Table 17). If a subject has only one copy of SMN1 , then the subject can be definitely identified as a carrier, but this variant only indicates a 2-8% chance of being a carrier if the copy number of SMN1 is two. Thanks to WGS, there are multiple variants that occur with different combinations of SMN1 and SMN2 copy numbers, and the identification of additional markers that can be used to improve our ability to identify these silent carriers. Additionally, loss of the c.840C>T splice variant currently explains approximately 95% of SMA cases, with the remaining cases involving other pathogenic variants. These other pathogenic variants represent another type of silent carrier. The method can directly genotype these other pathogenic variants as part of the testing process, further improving the ability to identify SMA carriers and cases.

[0232] Although there are complex regions of the genome in which normal WGS pipelines do not identify variants, this example demonstrates the ability to use WGS in combination with a targeted bioinformatics approach to identify one such complex region. This targeted strategy (WGS + custom bioinformatics) can be applied to a number of complex variants, such as repeat expansions and CYP2D6 , described herein. Traditionally, performing all known genetic tests and carrier screening on each individual has been cost-effective, so candidates for a particular genetic test have been identified using information such as carrier frequency and family history. However, this process meant that many people without a family history who might benefit from knowledge of SMA status did not have access to this data routinely. Once WGS analysis can accurately detect all SNVs and CNVs in all clinically relevant genes, then a more general and population-based genetic testing strategy with a single test will be possible. Improving WGS as a replacement for a single existing genetic test will help facilitate the integration of more genetic tests and screening tests into WGS, providing more general population-wide access to genetic testing. WGS provides a valuable opportunity to evaluate the entire genome for genetic variation, and continued development of more targeted bioinformatics solutions for complex areas with WGS data will help bring the promise of personalized medicine one step closer to reality.

Example 2

Precise genotyping of CYP2D6 using whole gene sequencing data

[0233] This example and Appendix A describe genotyping of CYP2D6 using whole genome sequencing data. The contents of Appendix A are incorporated herein by reference in their entirety.

[0234] CYP2D6 is involved in the metabolism of 25% of all drugs and is a key target for personalized medicine. Genotyping of CYP2D6 is challenging due to its high polymorphism, the presence of common structural variants (SVs), and high sequence similarity to the pseudogene paralog CYP2D7 . Described herein is a bioinformatics method, also referred to herein as Cyrius, that can accurately genotype CYP2D6 using whole genome sequencing (WGS) data. This method showed superior performance (97.9% true) compared to other methods (85.6-88.8%) in 138 samples with consensus determinations by GeT-RM and in 50 additional samples by Pacific Biosciences, California , Inc. (Menlo Park, California), also known as PacBio, sequencing data. A specific distinguishing feature of the method is the ability to determine structural variants of star alleles. The method correctly identified 97.2% (70/72) of structural variants of star alleles compared to 77.8-88.9% (56/72 and 64/72) for other methods. When applying the method to 2504 samples from the 1000 Genomes Project (1kGP), CYP2D6 star alleles involving SV were found to be 32.2% more common than previously reported for some populations. This example presents a benchmarking analysis with the largest validation data set. In some embodiments, the method is a useful tool for pharmacogenomic applications with WGS. The method could help bring the promise of precision medicine a step closer to reality.

Introduction

[0235] There is significant variation in the response of subjects to large numbers of clinically prescribed medications. A significant factor influencing this differential response to drugs is the genetic composition of drug-metabolizing genes. Precision medicine requires genotyping of pharmacogenes 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 defined by the Pharmacogene Variation (PharmVar) Consortium ( pharmvar.org/gene/CYP2D6 ). CYP2D6 star alleles are copies of the CYP2D6 gene, defined by a combination of small variants (such as single nucleotide variants (SNVs) and insertions/deletions (indels)) and structural variants (SV), and correspond to different levels of CYP2D6 enzymatic activity, such as weak, intermediate, normal or ultra-fast metabolism.

[0236] Genotyping of CYP2D6 is complicated by the presence of a non-functional paralogue, CYP2D7 , which is located upstream of CYP2D6 and shares 94% sequence similarity with several nearly identical regions. Deletions and duplications of CYP2D6 and fusions of CYP2D6 with its pseudogenic paralog CYP2D7 are common. Traditionally, CYP2D6 genotyping has been performed using chips or polymerase chain reaction (PCR)-based methods such as TaqMan assays, digital droplet PCR (ddPCR), and long-fragment PCR. These assays differ in the number of star alleles (variants) they target, resulting in variability in genotyping results across assays. A common limitation of these methods is that 1) the wild-type allele *1 is often determined by default when none of the target variants are found, or 2) a parental allele such as *2, which determines the true star allele, is not tested. These assays have low throughput and often make it difficult to detect structural variants.

[0237] Whole genome profiling at high throughput and in a clinically relevant time frame is possible with next generation sequencing (NGS). Large-scale population sequencing efforts are underway, and pharmacogenomic testing may be a desirable goal. Genotyping of CYP2D6 by NGS is particularly challenging due to common gene conversions between CYP2D6 and CYP2D7 (hereafter referred to as CYP2D6 /7), common SVs (gene deletions, duplications, and CYP2D6 /7 fusions), and the sequence similarity of CYP2D/7 , resulting in ambiguous read alignment for any of the genes. Some existing detectors cannot detect complex structural variants and have been shown to have low performance. Other existing determinants, such as Aldy (Numanagic et al. Allelic decomposition and exact genotyping of highly polymorphic and structurally variant genes. Nat Commun. 2018;9(1):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. Genet Med. 2019;21(2):361. Doi: 10.1038/s41436-018-0054 -0), relies on precise read alignment of CYP2D6 sequence reads to detect depth-based SVs and derive haplotype configurations based on observed minor variants and SVs. However, accurate reading of CYP2D6 sequences is often not possible at many positions throughout the gene because the sequence is very similar to CYP2D7 or even indistinguishable due to gene conversion. As a result, depth patterns may be ambiguous and determiners may identify false positive/negative minor variants. Some definers do not support hg38, so many studies will require realignment to hg37 to use these tools.

[0238] Availability of a panel of reference samples through the CDC Genetic Testing Reference Material 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), where consensus genotypes of major pharmacogenetic genes were obtained using multiple genotyping platforms, allowed assessment of genotyping accuracy for newly developed methods. GeT-RM covers 43 of the 106 star alleles of CYP2D6 . In addition, many of the methods used for these consensus genotypes may be prone to error, leading to conflicts between methods. The availability of high-quality long reads allows a complete picture of CYP2D6 to be obtained for better screening of complex variants and haplotypes. This paper describes Cyrius, a WGS-based genotyping method for CYP2D6 that overcomes problems with CYP2D6 and CYP2D7 (herein referred to as CYP2D6/7 ). Cyrius has superior genotyping accuracy compared to Aldy and Stargazer in 138 GeT-RM reference samples and 50 samples with PacBio whole genome sequencing data, covering 41 of the 106 known star alleles. The method was applied to high-depth sequencing data for 2504 unrelated samples from the 1000 Genomes Project (1kGP) to report the distribution of star alleles in five ethnic populations. This analysis demonstrates differences with PharmGKB frequencies, highlighting the potential pitfalls associated with combining limited star allele definitions generated using different technologies designed to identify specific subsets of known star alleles. This analysis expands the current understanding of the genetic diversity of CYP2D6 , especially complex star alleles with SV.

Materials and methods

Samples

[0239] The following were analyzed: WGS data for 138 GeT-Rm reference samples, including 96 samples that were genotyped in the initial GeT-RM study and updated in the latest version of GeT-RM, as well as 42 additional samples that were recently added to latest version of GeT-RM. For the first batch of 96 samples, WGS was performed using the TruSeq DNA PCR-free sample preparation kit using 150 bp paired-end reads sequenced on HiSeq X instruments from Illumina, Inc. (San Diego, CA, USA). The GRCh37 genome assembly was used for read alignment. Sequencing data for 70 of these samples was downloaded from ebi.ac.uk/ena/data/view/PRJEB19931. WGS data for a second batch of 42 samples was downloaded from the NYGC as part of the 1000 Genomes Project (see below).

[0240] For the population studies, data from the 1000 Genomes Project (1kGP) was used, for which WGS BAMs for 2504 samples were downloaded from ncbi.nlm.nih.gov/bioproject/PRJEB31736/. These BAM files were obtained by sequencing 2 x 150 bp reads. on Illumina NovaSeq 6000 instruments from libraries without PCR and alignment with the human reference sample hs38DH. WGS data for 70 GeT-RM samples were downloaded from ebi.ac.uk/ena/data/view/PRJEB19931.

PacBio sequencing

[0241] gDNA samples were purchased from the Coriell Institute for Medical Research (Coryell, NJ, USA). The quality of gDNA samples was assessed using Nanodrop (ThermoFisher, MA, USA). The A280/A260 ratio should be between 1.8 and 2.0, and the A260/230 ratio should be ≥ 2.0. The molecular weight of gDNA was assessed using the Femto Pulse system (Agilent CA, USA). Most DNA fragments should be >40 kb in size. If the quality of the gDNA sample from Coriell was lower than the protocol requirement, fresh DNA was extracted from B cells (Coryell, NJ, USA) using a Qiagen DNA extraction kit (Qiagen, CA, USA).

[0242] 10 μg of gDNA was fragmented to 15 kb. using Covaris g-tubes according to the manufacturer's instructions (Covaris, MA, USA). DNA was purified using 0.45x AMPure XP beads (Beckman Coulter, IN, USA) according to the manufacturer's instructions. The size of the cut DNA was confirmed using the Femto Pulse System (Agilent, California, USA).

[0243] Libraries were constructed in accordance with the PacBio protocol “Preparing HiFi SMRTbell® Libraries using SMRTbell Template Prep Kit 1.0” or “HiFi SMRTbell® Libraries using SMRTbell Express Template Prep Kit 2.0” (PacBio, California, USA). The library size was chosen to be 15~20 kb. using a Sage Elf instrument with 0.75% agarose (Sage Science, MA, USA). Quality control of all libraries was performed using Qubit (Life Technologies, CA, USA) and Femto Pulse (Agilent, CA, USA).

The PacBio Sequel II sequencing platform was used for sequencing. WGS data with 20× coverage are essentially obtained from 2~3 SMRT cells (Pacific Biosciences, CA, USA). CYP2D6 genotyping method

[0244] The method described in this example by Cyrius first determines the sum of the copy number (CN) of CYP2D6 /7, according to a method similar to that described in example 1. The number of reads was calculated directly from the BAM file, aligned to WGS based on all reads, mapped to CYP2D6 or CYP2D7 , including reads with zero mapping quality to account for regions of high sequence homology. The total number of reads was normalized by region length. GC correction was then performed against 3000 genome-wide preselected 2-kb regions. These 3000 normalization regions were randomly selected from the genome to provide stable coverage of population samples to reveal sequencing depth and capture GC bias. Population normalized depth values were modeled using a univariate mixture of 11 distributions that are centered around each integer copy number value representing copy number states ranging from 0 to 10. CYP2D6 + CYP2D7 CNS was defined using a normal distribution mixture model (GMM) with a threshold posterior probability 0.95. The same approach was used to determine the copy number of the 1.5-kb spacer region. between the REP7 repeat and CYP2D7 to reveal the copy number of fusion genes containing REP7 (FIG. 23).

[0245] In FIG. 23 is a non-limiting example plot showing the quality of WGS data in the CYP2D6 /7 region. The average mapping quality for the 1kGP samples is plotted for each position in the CYP2D6 /7 region. In a window 200 bp long. a median filter is used. REP6, REP7 and CYP2D6 /7 exons 9 are shown as boxes on the left ( CYP2D6 ) and right ( CYP2D7 ). Two 2.8 kb repeat regions. downstream, CYP2D6 (REP6) and CYP2D7 (REP7) are identical and essentially incompatible. The dotted frame indicates the spacer region between CYP2D7 and REP7. The two major homologous regions in the genes are shaded.

[0246] The method identified 118 differentiating bases of CYP2D6 / CYP2D7 (see supporting information for this example, FIG. 26). At each of these differentiating base positions, Cyrius determined the number of chromosomes carrying CYP2D6 and the number carrying CYP2D7 by combining the total copy number of CYP2D6 + CYP2D7 with the number of reads supporting each of the gene-specific bases. Based on the determined total copy number, Cyrius iterated through all possible combinations of CYP2D6 and CYP2D7 copy numbers and obtained the combination that gave the highest posterior probability for the observed number of CYP2D6 and CYP2D7 reads. Gene fusions were determined by identifying bases when the copy number of CYP2D6 changed (FIG. 27).

[0247] Cyrius analyzed read alignments to identify small variants that define star alleles. Variants of interest were divided into those that fell within homologous regions of CYP2D6 / CYP2D7 (i.e., low-quality mapping regions in FIG. 23) and variants that occurred in unique regions of CYP2D6 . For the first variant, Cyrius looked for the CYP2D6 read variant and its corresponding CYP2D7 site. In the latter case, Cyrius used reads aligned to CYP2D6 . The number of copies determined in the area was also taken into account when determining small variants. For example, where a duplication fusion *68 has been identified, one haplotype should have an intact copy of CYP2D6 plus a copy of *68, and the other haplotype should have an intact copy of CYP2D6 and therefore the CYP2D6 copy number should be 3 upstream of exon 2 and 2 downstream of exon 2. exon 2.

[0248] Finally, Cyrius mapped specific structural variants and minor variants according to the definition of star alleles (downloaded and analyzed from PharmVar, pharmvar.org/gene/ CYP2D6 , last accessed March 2019) to designate star alleles, which further grouped into haplotypes when, for example, there were more than two copies of CYP2D6. To do this, preliminary information was included to determine the exact haplotypes, for example *68 was in the same haplotype as *4, and *36 was in the same haplotype as *10). These preliminary values were generated from the tandem arrangement patterns described in PharmVar and are also supported by our hard data (12/12 for *68 and 25/25 for *36). An option was available to map only certain structural variants and minor variants to star alleles with known functions.

[0249] Of the 131 star alleles identified in PharmVar (last accessed March 2020), 25 are still awaiting curing, so the example excludes them and focuses on the 106 selected alleles (Cyrius has another option to include these uncured ones alleles). Of these 106 star alleles, four were removed from our target list, none of which were in GeT-RM. Deleted star alleles include *61 and *63 (both of unknown function), which are CYP2D6 /7 fusion genes very similar to *36, with the fusion breakpoint slightly higher. Since it was not possible to distinguish the exon 7-exon 8 region between CYP2D6 /7 (FIG. 26), these two star alleles cannot be distinguished from *36 and will be designated by Cyrius as *36. In addition, *27 (normal function) and *32 (unknown function) were removed; *27 and *32 share g.42126938C>T, a gene conversion variant in a highly homologous region (ideally the read would align with CYP2D7 . When counting reads supporting CYP2D6 and CYP2D7 at the same site, it may be difficult to accurately distinguish 1 copy of CYP2D6 and 3 copies of CYP2D7 from 2 copies each of 20. Thus, *27 would be called *1, and *32 would be called *41.

Checking the validity of GeT-RM results and long reading

[0250] When comparing the CYP2D6 definitions made by Cyrius, Aldy and Stargazer with the consensus genotypes provided by GeT-RM, the genotype was considered a match provided that all star alleles in the true genotype were present, the haplotype assignment was different. An example of this occurs in several samples listed in GeT-RM as *1/*10+*36+*36, but referred to by Aldy as *1+*36/*10+*36.

[0251] When validating genotype definitions from PacBio, PacBio reads covering the entire CYP2D6 gene were analyzed to identify small variants known to define star alleles. Long (~10 kb) reads allow given variants to be fully stepwise converted into haplotypes, and these haplotypes are compared to a star allele table to determine which star allele is read each time. reads carrying structural variations were identified by aligning reads to a set of control contigs that were constructed to represent known structural variants (*5/*13/*36/*68/duplications).

Application of Aldy and Stargazer

[0252] Aldy v2.2.5 was run using the command " aldy genotype -p lumina -g CYP2D6 ".

[0253] Stargazer v1.0.7 was run on the CYP2D6 genotype using VDR as the reference gene, with GDF and VCF files as input.

[0254] Since Aldy and Stargazer only support GRCh37, 1kGP samples were initially aligned to hs38DH using Isaac, and reorientation relative to GRCh37 was performed using Isaac.

results

Checking and comparing characteristics

[0255] CYP2D6 determinations by Cyrius, Aldy and Stargazer, based on 188 samples for which reliable, high-quality information was obtained. Whole-genome sequencing results of these 188 samples were compared, including 138 GeT-RM samples and 50 PacBio-confidence samples (Table 20, Table 21). PacBio CCS data allowed the detection and visualization of breakpoints of common and rare structural variants in the region (FIG. 24) and thus served as a valuable resource for the study of complex star alleles and confirmed the phasing of star allele variants. On short reads, these SV samples exhibited different depth signals that allowed accurate SV detection (FIG. 27).

Table 20. Summary of reliability comparison results

Determinant Total number of GeTRMs Total number of PacBio General Full approval Sensitivity Deletion
N=18 Duplication
N=14 Merger
N=40 No SV
N=116 Sensitivity, samples with SV Sensitivity, samples without SV Cyrius 138 50 188 184* 97.9% 17 14 39 114 97.2% 98.3% Aldy 167 88.8% 16 12 36 103 88.9% 88.8% Stargazer 161 85.6% 17 eleven 28 105 77.8% 90.5%

* After searching for three divergent samples, improvements were made to Cyrius, and Cyrius was able to accurately identify 187 of these 188 samples.

Table 21. Cyrius/Aldy/Stargazer results regarding the reliability of GeT-RM and PacBio data

Sample ID Credibility Cyrius Aldy Stargazer Reliable source 1kGP Type HG00276 *4/*5 *4/*5 *4/*5 *4/*5 GeT-RM Deletion NA10831 *4/*5 *4/*5 *4/*5 *4/*5 GeT-RM Deletion NA12873 *1/*5 *1/*5 exit with error *1/*5 GeT-RM Deletion NA17235 *1/*5 *1/*5 *1/*5 *1/*5 GeT-RM Deletion NA18855 *1/*5 *1/*5 *1/*5 *1/*5 GeT-RM Deletion NA18868 *2/*5 *2/*5 *2/*5 *2/*5 GeT-RM Deletion HG01706 *1/*5 *1/*5 *1/*5 *1/*5 PacBio Deletion HG00615 *2/*5 *2/*5 *2/*5 *2/*5 PacBio Deletion HG02523 *1/*5 *1/*5 *1/*5 *1/*5 PacBio Deletion NA18992 *1/*5 *1/*5 *1/*5 *1/*5 GeT-RM Deletion NA18861 *5/*29 *29/*5 *29/*5 *5/*29 GeT-RM Deletion NA19035 *2/*5 *2/*5 *2/*5 *2/*5 GeT-RM 1kGP Deletion NA18945 *1/*5 *1/*5 *1/*5 *1/*5 GeT-RM 1kGP Deletion HG03225 *5/*56 *10/*5
(updated to *56/*5) *5/*56 *5/*56 GeT-RM 1kGP Deletion HG03259 *5/*106 *106/*5 *106/*5 *5/*106 GeT-RM 1kGP Deletion HG03246 *5/*43 *43/*5 exit with error *5/*43 GeT-RM 1kGP Deletion NA19317 *5/*5 *5/*5 *5/*5 *2/*2 GeT-RM 1kGP Deletion NA18873 *5/*17 *17/*5 *17/*5 *5/*17 GeT-RM 1kGP Deletion HG00436 *2x2/*71 *2x2/*71 *2+*2/*71 *2x2/*71 GeT-RM duplication NA07439 *4x2/*41 *41/*4x2 *4+*4/*41 *4x2/*41 GeT-RM duplication NA17244 *2x2/*4x2
(+hybrid) *2x2/*4x2 *4+*4/*63+*78
+*2 *2x2/*4x2 GeT-RM duplication NA19226 *2/*2x2 *2/*2x2 *2/*2+*2 *2/*2x2 GeT-RM duplication NA24027 *2x2/*6 *2x2/*6 *2+*2/*6 *1x2/*6 GeT-RM duplication NA19920 *1/*4x2 *1/*4x2 *1/*4+*4 *1/*4x2 GeT-RM duplication NA19819 *2/*4x2 *2/*4x2 *2/*4+*4 *2/*4x2 GeT-RM duplication NA19207 *2x2/*10 *10/*2x2 *10/*2+*2 *2Dx2/*10 GeT-RM duplication NA17454 *1x2/*2x2 *1x2/*2x2 *1+*34/*2+*2 *1x2/*2x2 GeT-RM duplication NA19109 *2x2/*29 *2x2/*29 *2+*2/*29 *2x2/*29 GeT-RM duplication NA15245 *4x2/*4 *4/*4x2 *4/*4+*4 *4/*4x2 GeT-RM duplication HG00337 *2x2/*22 *2x2/*22 *2+*2/*22 *1/*2x2 GeT-RM 1kGP duplication HG01622 *1/*2x2 *1/*2x2 *1/*2+*2 *1/*34x2 PacBio duplication HG03131 *17/*2x2 *17/*2x2 *17/*2+*2 *2Dx2/*17 PacBio duplication HG01190 *68+*4/*5 *5/*4+*68 *4/*68 *5/*68+*4 GeT-RM merger NA12878 *3/*4+*68 *3/*4+*68 *3/*4+*68 *3/*4+*68 PacBio merger NA12877 *4/*4+*68 *4/*4+*68 *4/*4+*68 *4/*4+*68 PacBio merger NA21781 *2x2/*68+*4 *2x2/*4+*68 *2+*2/*68+*4 *2x2/*68+*4 GeT-RM merger HG01772 *4+*68/*4+*68 *4+*68/*4+*68 *68+*4/*68+*4 ungenotyped PacBio merger NA11832 *1/(*68)+*4 *1/*4+*68 *1/*68+*4 *1/*68+*4 GeT-RM merger NA12878 *3/(*68)+*4 *3/*4+*68 *3/*68+*4 *3/*68+*4 GeT-RM 1kGP merger NA12154 (*68)+*4/*33 *33/*4+*68 *33/*68+*4 *33x2/*68+*4 GeT-RM 1kGP merger HG00731 *4/*4+*68 *4/*4+*68 *4/*68+*4 *4/*68+*4 PacBio 1kGP merger HG00553 *29/*4+*68 *29/*4+*68 *29/*68+*4 *29/*68+*4 PacBio 1kGP merger NA23874 *4/*4 *4/*4+*68 *4/*68+*4 *4/*68+*4 GeT-RM merger NA24008 *1/*4 *1/*4+*68 *1/*68+*4 *1/*68+*4 GeT-RM merger NA18524 *1/*36x2+*10 *1/*10+*36
+*36 *1+*36/*36
+*10 *36+*10/*36
+*10 GeT-RM merger NA18526 *1/*36x2+*10 *1/*10+*36
+*36 *1+*36/*36
+*10 *36+*10/*36
+*10 GeT-RM merger NA18540 (*36+)10/*41 *41/*10+*36
+*36 *36+*10/*61
+*69 *36+*10/*36
+*10 GeT-RM merger NA18564 *2A/*36+*10 *2/*10+*36 *2/*36+*10 *2/*36+*10 GeT-RM merger NA18565 *10/*36x2 *10/*10+*36 *10/*36+*10 *10/*36+*10 GeT-RM merger NA18617 *36+*10/*36
+*10 *36+*10/*36
+*10 *36+*10/*36
+*10 *36+*10/*36
+*10 GeT-RM merger NA18959 *2/*36+*10 *2/*10+*36 *2/*36+*10 *2/*36+*10 GeT-RM merger NA23246 *10x2/*36
+*10 *10x2/*10
+*36 *10+*10/*36
+*10 *10x2/*36+*10 GeT-RM merger NA18980 *2/*36+*10 *2/*10+*36 *2/*36+*10 *2/*36+*10 GeT-RM merger NA18642 *36+*10/*1
+*90 *36+*10/*1
+*90 *1+*90/*36
+*10 *1x2/*36+*10 GeT-RM 1kGP merger HG00463 *36+*10/*36
+*10 *36+*10/*36
+*10 *36+*10/*36
+*10 *36+*10/*36
+*10 GeT-RM 1kGP merger HG02373 *14/*36+*10 *14/*10+*36 *14/*36+*10 *14/*36+*10 GeT-RM 1kGP merger NA18572 *36+*10/*41 *41/*10+*36 *10+*83/*69 *36+*10/*41 GeT-RM 1kGP merger NA18632 *36x2+*10/
*52 *52/*10+*36
+*36 *36+*10/*36
+*52 *36+*10/*36
+*10 GeT-RM 1kGP merger NA18563 *1/*36+*10 *1/*10+*36 *1/*36+*10 *1/*36+*10 GeT-RM 1kGP merger NA18545 *5/*36x2
+*10x2 *36+*10/*36
+*10 *36+*10/*36
+*10 *36+*10/*36
+*10 GeT-RM 1kGP merger HG02068 *10/*10+*36 *10/*10+*36 *10/*36+*10 *10/*36+*10 PacBio merger HG00612 *10/*10+*36 *10/*10+*36 *10/*36+*10 *10/*36+*10 PacBio merger HG00597 *49/*10+*36
+*36+*83 *49/*10+*36
+*36+*83 *36+*10+*49
/*36+*39 *36+*10/*39x3 PacBio merger HG02015 *10/*10+*36 *10/*10+*36 *10/*36+*10 *10/*36+*10 PacBio merger HG02071 *2/*10+*36 *2/*10+*36 *2/*36+*10 *2/*36+*10 PacBio merger HG02129 *1/*10+*36 *1/*10+*36 *1/*36+*10 *1/*36+*10 PacBio merger HG02074 *10+*36/*10
+*36 *36+*10/*36
+*10 *36+*10/*36
+*10 *36+*10/*36
+*10 PacBio merger HG00844 *10+*36/*10
+*36 *36+*10/*36
+*10 *36+*10/*36
+*10 *36+*10/*36
+*10 PacBio 1kGP merger HG005
(NA24631) *49/*10+*36 *49/*10+*36 *49/*10+*36 *10/*36+*10 PacBio merger NA19785 *1/*13+*2 *2+*13/*1 *1/*79+*2 ungenotyped GeT-RM 1kGP merger HG00290 *2+*13/*1 *2+*13/*1 *1/*79+*2 ungenotyped PacBio 1kGP merger HG00421 *2/*10x2 *10+*36/*2
(updated to
*10x2/*2) *2/*36+*10 *2/*36+*10 GeT-RM, PacBio 1kGP fusion,* 10D represents fusion, see FIG. 27 HG00589 *1/*21 *1/*21 *1/*21 *1/*21 GeT-RM No SV NA06991 *1/*4 *1/*4 *1/*4 *1/*4 GeT-RM No SV NA07000 *2 (*35)/*9 *35/*9 *35/*9 *9/*35 GeT-RM No SV NA07019 *1/*4 *1/*4 *1/*4 *1/*4 GeT-RM No SV NA07029 *1/*35 *1/*35 *1/*35 *1/*35 GeT-RM No SV NA07055 *4/*4 *4/*4 *4/*4 *4/*4 GeT-RM No SV NA07056 *2/*4 *2/*4 *2/*4 *2/*4 GeT-RM No SV NA07348 *1/*6 *1/*6 *1/*6 *1/*6 GeT-RM No SV NA07357 *1/*6 *1/*6 *1/*6 *1/*6 GeT-RM No SV NA10847 *1/*41 *1/*41 *1/*41 *1/*41 GeT-RM No SV NA10851 *1/*4 *1/*4 *39/*4 *1/*4 GeT-RM No SV NA10854 *1/*4 *1/*4 *1/*4 *1/*4 GeT-RM No SV NA11839 *1/*2 *1/*2 *1/*2 *1/*2 GeT-RM No SV NA11993 *1/*9 *1/*9 *1/*9 *1/*9 GeT-RM No SV NA12003 *4/*35 *35/*4 *35/*4 *4/*35 GeT-RM No SV NA12006 *4/*41 *41/*4 *4/*41 *4/*41 GeT-RM No SV NA12145 *1/*4 *1/*4 *1/*4 *1/*4 GeT-RM No SV NA12156 *1/*4 *1/*4 *1/*4 *1/*4 GeT-RM No SV NA12236 *1/*4 *1/*4 *1/*4 *1/*4 GeT-RM No SV NA12717 *1/*1 *1/*1 *1/*1 *1/*1 GeT-RM No SV NA12813 *2/*4 *2/*4 *2/*4 *2/*4 GeT-RM No SV NA17074 *1/*2 *1/*2 *1/*2 *1/*2 GeT-RM No SV NA17102 *1/*40 *1/*40 *1/*40 *1/*40 GeT-RM No SV NA17204 *1/*35 *1/*35 *1/*35 *1/*35 GeT-RM No SV NA17227 *1/*9 *1/*9 *1/*9 *1/*9 GeT-RM No SV NA17234 *1/*41 *1/*41 *1/*41 *1/*41 GeT-RM No SV NA17448 *1/*28 *1/*28 *1/*28 *1/*28 GeT-RM No SV NA17641 *2/*35 *2/*35 *2/*35 *2/*35 GeT-RM No SV NA17642 *1/*1 *1/*1 *1/*61 *1/*1 GeT-RM No SV NA17657 *4/*9 *4/*9 *4/*9 *4/*9 GeT-RM No SV NA17658 *1/*2 *1/*2 *1/*2 *1/*2 GeT-RM No SV NA17660 *1/*2 *1/*2 *1/*2 *1/*2 GeT-RM No SV NA17673 *1/*4 *1/*4 *1/*4 *1/*4 GeT-RM No SV NA17679 *1/*4 *1/*4 *1/*4 *1/*4 GeT-RM No SV NA17702 *1/*35 *1/*35 *1/*35 *1/*35 GeT-RM No SV NA18484 *1/*17 *1/*17 *61-like/*78 *1/*17 GeT-RM No SV NA18509 *2/*17 *17/*2 *17/*2 *2/*17 GeT-RM No SV NA18518 *17/*29 *17/*29 *17/*29 *17/*29 GeT-RM No SV NA18519 *1/*29 *106/*29 *106/*29 *29/*106 GeT-RM No SV NA18544 *10/*41 *10/*41 *10/*41 *10/*41 GeT-RM No SV NA18552 *1/*14 *1/*14 *1/*14 *1/*14 GeT-RM No SV NA18942 *2/*2 *2/*2 *2/*2 *2/*2 GeT-RM No SV NA18952 *2/*2 *2/*2 *2/*2 *2/*2 GeT-RM No SV NA18966 *1/*2 *1/*2 *1/*2 *1/*2 GeT-RM No SV NA18973 *1/*21 *1/*21 *1/*21 *2D/*21 GeT-RM No SV NA19003 *1/*1 *1/*1 *1/*1 *1/*1 GeT-RM No SV NA19007 *1/*1 *1/*1 *1/*1 *1/*1 GeT-RM No SV NA19095 *1/*29 *1/*29 *1/*29 *1/*29 GeT-RM No SV NA19122 *2/*17 *17/*2 *17/*2 *2/*17 GeT-RM No SV NA19143 *2 (*45)/*10 *10/*45 *10/*45 *2/*10 GeT-RM No SV NA19147 *17/*29 *17/*29 *17/*29 *17/*29 GeT-RM No SV NA19174 *4/*40 *40/*4 *4/*40 *4/*40 GeT-RM No SV NA19176 *1/*2 *1/*2 *1/*2 *1/*2 GeT-RM No SV NA19178 *1/*1 *1/*1 *1/*1 *1/*1 GeT-RM No SV NA19213 *1/*1 *1/*1 *1/*1 *1/*1 GeT-RM No SV NA19239 *15/*17 *15/*17 *15/*17 *2/*15 GeT-RM No SV NA19789 *1/*1 *1/*1 *1/*61 *1/*1 GeT-RM No SV NA19908 *1/*46 *1/*46;*43/*45 *43/*45 *1/*46 GeT-RM No SV NA19917 *1/*40 *1/*40 *1/*40 *1/*40 GeT-RM No SV NA20296 *1/*2 *1/*2 *1/*2 *1/*2D GeT-RM No SV NA20509 *4/*35 *35/*4 *35/*4 *4/*35 GeT-RM No SV NA23275 *1/*40 *1/*17
(updated to *1/*40) *1/*40 *1/*40 GeT-RM No SV NA23348 *7/*35 *35/*7 *35/*7 *7/*35 GeT-RM No SV HG03882 *1/*112 *1/*112 *1/ 1-like+*61 *1/*1 GeT-RM 1kGP No SV HG03780 *1/*112 *1/*112 *1/*1-like *1/*1 GeT-RM 1kGP No SV NA19238 *1/*17 *1/*17 *1/*17 *1/*17 GeT-RM, PacBio 1kGP No SV NA20803 *2/*22 *2/*22 *2/*22 *1/*2 GeT-RM 1kGP No SV HG04206 *2/*113 *113/*2 *1/*2 *2/*113 GeT-RM 1kGP No SV HG01108 *2/*106 *106/*2 *106/*2 *2/*106 GeT-RM 1kGP No SV NA20875 *1/*111 *1/*111 *111/*2 *1/*2 GeT-RM 1kGP No SV HG01094 *1/*31 *1/*31 *1/*31 *1/*31 GeT-RM 1kGP No SV HG01086 *1/*31 *1/*31 *1/*31 *1/*31 GeT-RM 1kGP No SV NA07048 *1/*4 *1/*4 *10/*74-like *1/*4 GeT-RM 1kGP No SV HG03703 *1/*99 *1/*99 *1/*10 *1/*99 GeT-RM 1kGP No SV NA20289 *6/*11 *11/*6 *11/*6 *6/*11 GeT-RM 1kGP No SV NA19700 *4/*29 *29/*4 *29/*4 *4/*29 GeT-RM 1kGP No SV HG00373 *2/*2 *2/*2 *2/*2 *2/*2 GeT-RM 1kGP No SV NA21105 *3/*111 *111/*3 *2/*3-like *2/*3 GeT-RM 1kGP No SV NA11881 *2/*3 *2/*3 *2/*3 *2/*3 GeT-RM 1kGP No SV NA12815 *2/*41 *2/*41 *2/*41 *2/*41 GeT-RM 1kGP No SV HG01680 *28/*59 *28/*59 *28/*59 *28/*59 GeT-RM 1kGP No SV HG03643 *2/*7 *2/*7 *2/*7 *2/*7 GeT-RM 1kGP No SV HG03781 *2/*99 *2/*99 *10/*2 *2/*99 GeT-RM 1kGP No SV HG00111 *3/*3 *3/*3 *3/*3 *3/*3 GeT-RM 1kGP No SV NA06989 *9/*9 *9/*9 *9/*9-like *9/*9 GeT-RM 1kGP No SV NA19777 *1/*82 *1/*82 *1/*82 *1/*1 GeT-RM 1kGP No SV HG02723 *17/*2 *17/*2 *17/*2 *2/*17 PacBio No SV HG03522 *1/*1 *1/*1 *1/*1 *1/*1 PacBio No SV HG00450 *10/*41 *10/*41 *10/*41 *10/*41 PacBio No SV HG03453 *1/*29 *1/*29 *1/*29 *1/*29 PacBio No SV HG01687 *1/*6 *1/*6 *1/*6 *1/*6 PacBio No SV HG02984 *2/*4 *2/*4 *2/*4 *2/*4 PacBio No SV HG01763 *1/*1 *1/*1 *1/*1 *1/*1 PacBio No SV HG03098 *2/*29 *2/*29 *2/*29 *2/*29 PacBio No SV HG03041 *17/*29 *17/*29 *17/*29 *17/*29 PacBio No SV HG02622 *17/*46 *17/*46 *61-like/*78 *17/*46 PacBio No SV HG01621 *2/*33 *2/*33 *2/*33 *2/*33 PacBio No SV HG03579 *1/*2 *1/*2 *1/*2 *1/*2 PacBio No SV HG02975 *17/*29 *17/*29 *17/*29 *17/*29 PacBio No SV HG03101 *1/*17 *1/*17 *1/*17 *1/*17 PacBio No SV HG03065 *106/*29 *106/*29 *106/*29 *29/*106 PacBio No SV HG03486 *1/*17 *1/*17 *1/*17 *1/*17 PacBio No SV HG03308 *1/*29 *1/*29 *1/*29 *1/*29 PacBio No SV HG00513 *10/*10 *10/*10 *10/*10 *10/*10 PacBio 1kGP No SV HG00143 *1/*4 *1/*4 *1/*4 *1/*4 PacBio 1kGP No SV NA20527 *2/*4 *2/*4 *2/*4 *2/*4 PacBio 1kGP No SV HG00732 *41/*9 *41/*9 *41/*9 *9/*41 PacBio 1kGP No SV HG01119 *1/*4 *1/*4 *1/*4 *1/*4 PacBio 1kGP No SV HG00554 *4/*4 *4/*4 *4/*4 *4/*4 PacBio 1kGP No SV HG01254 *2/*41 *2/*41 *2/*41 *2/*41 PacBio 1kGP No SV HG00186 *2/*4 *2/*4 *2/*4 *2/*4 PacBio 1kGP No SV HG00263 *1/*35 *1/*35 *1/*35 *1/*2 PacBio 1kGP No SV NA19239 *15/*17 *15/*17 *15/*17 *1/*15 PacBio 1kGP No SV NA19437 *17/*2 *17/*2 *17/*2 *2/*17 PacBio 1kGP No SV NA19449 *1/*17 *1/*17 *1/*17 *1/*17 PacBio 1kGP No SV HG002
(NA24385) *2/*4 *2/*4 *2/*4 *2/*4 PacBio No SV

[0256] Through comparison with GeT-RM samples, three samples were found in which the definitions of all three determinants were consistent, but not consistent with the GeT-RM consensus. PacBio whole genome sequencing confirmed that the three definitions were correct and the GeT-RM consensus needed to be updated (FIG. 24).

[0257] In FIG. Figure 24 shows structural variants confirmed by PacBio CCS reads. PacBio read confirms deletion (*5), duplication, and fusion (*36, *68, and *13). Graphs were generated using SV-VIZ2 (zotero.org/google-docs/?xAunA6 ⁾ . In the case of deletions and duplications, due to identical sequence in the REP6/7 regions, the exact positions of the breakpoints within REP6/7 were not available. The break points in A and B are for illustration purposes only. The genotypes of the samples in panel AE are *2/*5, *17/*2x2, *10/*10+*36, *29/*4+*68 and *1/*2+*13, respectively.

[0258] Cyrius initially made four determinations diverging from the GeT-RM confidence, showing a sensitivity of 97.9%. Included among these discrepancies was specimen NA19908 (GeT-RM defined *1/*46), in which Cyrius identified 1/*46 and *43/*45 as two possible diplotypes. Both of these two star allele combinations produce the same set of options. Neither read phase analysis nor population frequency analysis can rule out any combination of genotypes. The genotyping results from the various analyzes that produced the GeT-RM consensus for this sample also showed a discrepancy between *1/*46 and* *43/*45, indicating the complexity of these combinations (Table 22). Future sequencing of more samples of either diplotype may help identify new variants that differentiate them.

Table 22. GeT-RM results for sample NA19908

Consensus genotype TaqMan
+CNV
+XL-PCR PharmacoScan
v.r6 PharmacoScan
custom
v.r6+20180103 iPLEX
CYP2D6 V1.1 iPLEX
CYP2D6
V1.1
+individual panel and VeriDose Sanger, NGS or SMRT sequencing *1/*46 *1/*45 *1/*2 *1/*46, *43/*45 *1/*46 N/A *46 (ASXL-PCR) Sanger; NGS

[0259] In the remaining three samples where Cyrius was inconsistent, errors were identified and Cyrius was improved to identify the correct genotypes. First, in NA23275 (*1/*40), the 18-bp insertion defining *40 was initially missed because reads containing the insertion were often aligned not as having the insertion, but as soft stitches. The qualifier has been improved to take into account soft crosslinks when searching for a variant. Second, in HG03225 (*5/*56), reads derived from CYP2D7 were aligned to CYP2D6 , preventing detection of the defining variant*56. The determinant was improved to be more sensitive to variant reads in a given region. Finally, in HG00421 (*10x2/*2), the merge was erroneously designated as *36, as were the other two qualifiers. Closer examination of this sample using PacBio data revealed a different fusion, *10D, with the fusion breakpoint located downstream of exon 9 (FIG. 28). This fusion performs the same function as *10 (reduced function), while *36 is nonfunctional due to exon 9 derived from CYP2D7 . The determinant has been improved to be able to determine *10D. Although these three samples were treated as misidentifications in this example, improvements made to Cyrius following the search for these three samples allowed 187 of the 188 samples to be accurately identified, highlighting how more hard data and more population data can reveal limitations. which can provide improvements to the determinant for subsequent samples.

[0260] On the contrary, both other determinantsCYP2D6 had a sensitivity of less than 90% compared to these samples. Aldy had a sensitivity of 88.8%. Specifically, it redefined mergesCYP2D6/CYP2D7, such as *61, *63, *78 and *83 (identified 8 of 21 divergent samples, Table 21). The fusion identified by Aldy may be refuted by the PacBio data in FIG. 29. Stargazer had a sensitivity of 85.6% and was most error prone in the presence of SV. The sensitivity of samples with SV was only 77.8%, and 16 of the 27 divergent determinations were in samples with structural variants. It should be noted that he incorrectly identified NA19317 (*5/*5) as *2/*2, while the double deletion was completely absent. Stargazer was unable to genotype two samples using the *13 fusion (Table 21). In addition, Stargazer showed a high error rate when merging *36 (7 incorrect identifications out of 25 of all samples with *36). Specifically, Stargazer incorrectly identified all 5 samples that had more than one copy of *36 in the same haplotype.

[0261] Collectively, the 188 validation samples used in this example confirmed the accuracy of Cyrius CYP2D6 detection in 48 different haplotypes (Table 23), including 41 star alleles, as well as several common and rare SV structures such as duplications, *2+* 13, *4+*68, *10+*36, *10+*36+*36 and *10+*36+*36+*83 (new haplotype not previously reported, see FIG. 30A and 30B). These 41 star alleles that were tested in data screening represent 38.7% of the 106 curated star alleles currently listed in PharmVar and 53.4% (31 of 58) of those with known function. They overlap 96.4% of the haplotypes identified by Cyrius in the 1kGP samples (Table 23, also see next section).

Table 23. Haplotypes confirmed in this example and their frequency in 1kGP

Haplotype Pan-
ethnicity European population Mixed Americans People from
East Asia African population South Asians Confirmed in this example In the complete GeT-RM set Function *1 33.43 35.79 45.97 26.19 26.25 39.26 x x Normal *2 14.86 16.2 18.44 7.74 13.24 20.45 x x Normal *3 0.54 1.79 0.58 0 0.23 0.2 x x No *4 5.83 11.83 8.79 0.2 2.34 8.08 x x No *5 3.49 2.39 2.02 3.47 5.82 2.56 x x No *6 0.5 2.09 0.29 0 0.08 0.1 x x No *7 0.18 0 0 0 0 0.92 x x No *9 0.7 2.49 1.3 0 0.08 0 x x Reduced *10 5.41 1.39 1.44 15.08 4.39 3.78 x x Reduced *eleven 0.02 0 0 0 0.08 0 x x No *13 0.1 0.2 0.14 0 0.08 0.1 x x No *14 0.18 0 0 0.89 0 0 x x Reduced *15 0.06 0 0 0 0.23 0 x x No *17 5.25 0.2 0.86 0 19.29 0 x x Reduced *21 0.1 0 0 0.5 0 0 x x No *22 0.06 0.3 0 0 0 0 x x Unknown *28 0.12 0.5 0.14 0 0 0 x x Unknown *29 2.64 0 0.29 0 9.83 0 x x Reduced *31 0.12 0.2 0.58 0 0 0 x x No *33 0.18 0.6 0.29 0 0 0.1 x x Normal *34 0.02 0 0 0 0.08 0 Normal *35 1.48 4.77 2.45 0 0.23 0.61 x x Normal *36 0.1 0 0 0.1 0.3 0 No *39 0.08 0 0.14 0 0.08 0.2 x Normal *40 0.24 0 0 0 0.91 0 x x No *41 6.15 9.05 6.05 3.77 1.59 11.86 x x Reduced *43 0.5 0.1 0 0 1.06 1.02 x x Unknown *45 0.88 0 0.29 0 3.18 0 x x Normal *46 0.14 0 0.14 0 0.45 0 x x Normal *49 0.1 0 0 0.5 0 0 x Reduced *52 0.02 0 0 0.1 0 0 x x Unknown *56 0.02 0 0 0 0.08 0 x x No *59 0.06 0.2 0.14 0 0 0 x x Reduced *71 0.12 0 0 0.6 0 0 x x Unknown *82 0.06 0 0.43 0 0 0 x x Unknown *84 0.02 0 0 0 0.08 0 Reduced *86 0.44 0 0 0 0 2.25 Unknown *99 0.04 0 0 0 0 0.2 x x No *106 0.32 0 0.14 0 1.13 0 x x Unknown *108 0.06 0.3 0 0 0 0 x Unknown *111 0.16 0 0 0 0 0.82 x x Unknown *112 0.04 0 0 0 0 0.2 x x Unknown *113 0.16 0 0 0 0 0.82 x x Unknown *1x2 0.5 0.5 1.15 0.1 0.45 0.51 x x Increased *1x3 0.02 0 0 0 0.08 0 Increased *2x2 1.14 1.49 0.58 0.6 2.12 0.41 x x Increased *2x3 0.04 0.1 0 0 0.08 0 Increased *4x2 0.84 0.3 0.14 0 2.87 0 x x No *4x3 0.04 0 0 0 0.15 0 No *9x2 0.02 0.1 0 0 0 0 Normal *10x2 0.06 0 0 0.3 0 0 x x Reduced *17x2 0.02 0 0 0 0.08 0 x Normal *29x2 0.1 0 0 0 0.38 0 Normal *35x2 0.02 0 0.14 0 0 0 Increased *43x2 0.04 0 0.14 0 0.08 0 Unknown *45x3 0.02 0 0 0 0.08 0 Increased *10+*36 7.23 0 0.14 34.62 0.08 1.12 x x Reduced *4+*68 1.94 5.57 2.45 0 0.23 2.15 x x No *4+*68+*68 0.08 0.1 0.43 0 0 0 No *10+*36+*36 0.32 0 0 1.59 0 0 x x Reduced *10+*36+*36+*36 0.02 0 0 0.1 0 0 x x Reduced *2+*13 0.06 0.2 0.14 0 0 0 x x Normal *4+*4N 0.14 0.7 0 0 0 0 x No *1+*90 0.02 0 0 0.1 0 0 x x Unknown *10+*36+*36+*83 0.02 0 0 0.1 0 0 x Reduced Unknown 2.36 0.6 3.75 3.37 2.27 2.25 % of haplotypes overlapping the validation set 96.4 98.1 95.4 96.5 96.3 95.3

Frequency of CYP2D6 haplotypes in five ethnic populations

[0262] Given the high accuracy reported in the previous section, Cyrius was used outside of the validation sets to study CYP2D6 in the general population. The distribution of haplotypes by population (European, African, East Asian, South Asian and mixed American consisting of Colombian, Mexican, Peruvian and Puerto Rican) in 2504 1kGP samples was analyzed (FIG. 25, Table 23). Cyrius finally identified diplotypes in 2445 (97.6%) of 2504 samples, including 46 distinct star alleles, with 41 star alleles overlapping with those included in the validation data. These 41 previously confirmed star alleles represent 96.5% of all star alleles that were detected in the 1kGP samples (Table 23).

[0263] In FIG. 25 shows a non-limiting example plot showing CYP2D6 allele frequencies for five ethnic populations for the ten most common CYP2D6 altered function haplotypes. One haplotype (*2x2) is characterized by increased function, two haplotypes (*4 and *4 + *68) by lack of function, and the remaining haplotypes by decreased function.

[0264] In the 59 samples where Cyrius did not make a definitive diplotype determination, 10 samples had an inconclusive SV determination, 30 samples had variant determinations that did not match any known star allele, four samples had the same ambiguity between *1/*46 and * 43/*45 as described above for validation sample NA19908, and 15 samples had final star allele determinations that Cyrius could not unambiguously convert to diplotypes.

[0265] In most cases, the haplotype frequencies are consistent with pharmGKB (FIGS. 31A and 31B, Table 24). For example, Africans have a high frequency of *17 (~20%) and *29 (~10%), South Asians have a high frequency of *41 (~12%), Europeans have a high frequency of *4 (18-20%, including *4+*68), and East Asians have a high frequency of *10 (40-50%, including *10+*36). The increased sensitivity of Cyrius to structural variants provides a more complete picture of the frequencies of structural variants in populations. Among these, a haplotype containing the *10 + *36 fusion is very common in East Asians (>30% compared to 1-2% reported in PharmGKB, FIGS. 31A and 31B), as well as another haplotype containing merge *4 +. *68 is also quite common in Europeans (>5%, data not available in PharmGKB, FIGS. 31A and 31B). Taken together, SV-adjusted haplotype frequencies were 32.2%, 5.57%, 1.47%, 1.34%, and 0.45% greater than those reported in PharmGKB in East Asians, Europeans, Americans, and Africans and South Asia, respectively (overall frequencies in PharmGKB are 7.48%, 5.33%, 5.17%, 9.9% and 6.19%, respectively).

[0266] There are several other haplotypes for which lower frequencies were observed than in PharmGKB (FIGS. 31A and 31B), highlighting the difficulty of combining data from multiple studies using different technologies. These include *2 in Africans and South Asians. Since *2 is the default assignment, if some other star alleles are not tested, its frequency may be overestimated in PharmGKB. A lower frequency of *41 was found in Africans. According to PharmGKB, *41 has not always been consistently identified by the detector SNP across studies, leading to an overestimation of the frequency of *41, especially in individuals of African descent. The much higher frequency of *29 in South Asians in PharmGKB (6% versus 0% estimated in this example) was due to an error in PharmGKB: 0.2% in the publication of Sistonen et al. ( CYP2D6 worldwide genetic variation shows high frequency of altered activity variants and no continental structure. Pharmacogenet Genomics. 2007;17(2):93-101. doi:10.1097/01.fc.0000239974.69464.f2) was erroneously included in PharmGKB as 20%. Among Europeans, the frequency of *34 and *39 was much lower. *34 and *39 are defined by one of two variants that define 2, so both of these two variants had to be tested in any study describing CYP2D6 . *34 and *39 are present at levels >1% in only 3 of 91 studies of Europeans in PharmGKB, among which Wesmiller et al. (The Association of CYP2D6 Genotype and Postoperative Nausea and Vomiting in Orthopedic Trauma Patients. Biol Res Nurs. 2013; 15 (4): 382-389. doi: 10.1177/1099800412449181), who reported only *39 and had 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; 17(2): 5-14. doi: 10.2478/BJMG-2014-0069) reported both *34 and *39 and were specific to the country of 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; 9. doi: 10.3389/far.2018.00305), did not report *34 or *39, but PharmGKB may have mistakenly assumed the frequency reported for * 35, for frequency *34.

CYP2D6 / CYP2D7 Differentiating Base Analysis

[0267] A total of 208 single nucleotide differences between CYP2D6 /7 were extracted from the reference genome. In 1kGP samples, where the total copy number of CYP2D6 + CYP2D7 is 4, i.e. no structural changes were detected, the percentage of samples in which the CYP2D6 base copy number was determined to be 2 at 208 sites was queried (FIG. 26). A small percentage of samples with two copies of the CYP2D6 base were observed at many sites, suggesting that the CYP2D6 / CYP2D7 base difference is not fixed in the population, so base differences cannot be used to distinguish between the two genes. Using site-by-site read alignment can introduce significant noise when differentiating two genes. A total of 118 highly stable sites were selected where >98% of samples showed two copies of CYP2D6 bases for CYP2D6 / CYP2D7 differentiation, resulting in a clean signal for SV detection.

Additional figures and tables

[0268] In FIG. 26 shows that sites of CYP2D6 / CYP2D7 base differences are highly variable in the population. The y-axis shows the frequency of samples in which the CN for the CYP2D6 base is identified in 2 of all samples having a total CYP2D6 + CYP2D7 CN of 4. The x-axis shows the genome coordinates in hg38. CYP2D6 exons are shown as gray boxes above the graph. The black horizontal line indicates the 98% cutoff.

[0269] In FIG. Figure 27 shows unprocessed CYP2D6 CNs at CYP2D6/7 differentiation sites in SV examples. Raw CYP2D6 CNs were calculated as the total CYP2D6+CYP2D7 CNs multiplied by the ratio of CYP2D6 supporting reads from CYP2D6 to CYP2D7 supporting reads. The large diamond indicates the copy number of genes derived from CYP2D6 at the end of the gene (can be a complete CYP2D6 gene or a fusion gene ending with CYP2D6 ), calculated as the total CN of CYP2D6+CYP2 D7 minus the CN of the CYP2D7 spacer region (see FIG. 23). To detect SV, the CN of CYP2D6 was determined at each site, and a change in the CN of CYP2D6 in the gene indicated the presence of SV. For example, in HG01161, CN CYP2D6 changed from 2 to 1 between exon 7 and exon 9, indicating a CYP2D7-CYP2D6 hybrid gene. In HG00553, CN CYP2D6 changed from 2 to 3 between exon 1 and exon 2, indicating a CYP2D6-CYP2D7 hybrid gene.

[0270] In FIG. 28 shows that PacBio data support the *10D fusion in HG00421. For comparison, a sample with *36 (HG00612) is shown. PacBio reads containing fusions are reads with shaded bases, bases that represent a software stitch made by the aligner and were derived from part of the CYP2D7 fusion. The breakpoints of the fusions are close to each other, but the breakpoint for *36 is located upstream of the base differences in exon 9 (located within the black block), and the breakpoint for *10D is located downstream, leaving the CYP2D6 gene intact.

[0271] In FIG. 29 shows that the PacBio data had a false *61 ( CYP2D6 / CYP2D7 hybrid) produced 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). The PacBio data showed no structural variant in this region (each read was fully aligned, with no soft stitches indicating unaligned parts).

[0272] In FIG. 30A and 30B show the new haplotype *10+*36+*36+*83 in HG00597. FIG. 30A The depth chart shown in FIG. 27 shows that HG00597 had three copies of *36-like fusions, all of which had a breakpoint in the homologous region between exon 7 and exon 9. FIG. 30B Screenshot of IGV data from PacBio, showing all reads containing fusions, that is, those aligned by soft stitching. One copy of the fusion gene lacked g.42130692G>A, a SNP that was located at *36 but not at *83, as shown in the region flanked by the two black vertical lines. This copy was *83, and unlike what was reported in PharmVar, it was a fusion gene with REP7 rather than REP6, otherwise the copy number of the region downstream of exon 9 would have been 3 instead of 2 in FIG. 30A.

[0273] In FIG. 31A and 31B compared the frequencies of 1kGP and pharmGKB. Each dot represents a haplotype with a frequency >=0.5% for 1kGP or pharmGKB. SV-associated haplotypes were noted, including the two haplotypes with the greatest deviation (*10+*36 in East Asians and *4+*68 in Europeans). Other haplotypes with deviated values are labeled (*2, *41, *34, *39, *2, and *29). A diagonal line is drawn for each panel. Correlation coefficients are given for each population (*10+*36 excluded for East Asians and *4+*68 excluded for Europeans for calculation). In FIG. 31B shows values in the lower range of values (<5%).

[0274] In FIG. 32 shows a non-limiting example of an IGV snapshot showing de novo assembly of PacBio reads into HG00733 not including the *68 fusion.

Table 24. Comparison of the haplotype frequencies of certain Cyrius and the frequencies of pharmGKB

Haplotype Ethnicity pharmGKB 1kGP Function *1 African population 9.53 26,786 Normal function *2 African population 18.81 13,509 Normal function *3 African population 0.15 0.233 No function *4 African population 3.33 2,407 No function *6 African population 0 0.078 No function *9 African population 0 0.078 Decreased function *10 African population 6.71 4,814 Decreased function *15 African population 0.57 0.233 No function *17 African population 19.58 19,798 Decreased function *29 African population 10.73 10,093 Decreased function *35 African population 0 0.233 Normal function *39 African population 0 0.078 Normal function *40 African population 1.31 0.932 No function *41 African population 11.47 1.553 Decreased function *43 African population 0.96 1,087 Undefined function *45 African population 5.77 3.261 Normal function *46 African population 0 0.466 Normal function *1x2 African population 1.12 0.466 Enhanced function *2x2 African population 1.73 2,174 Enhanced function *4x2 African population 1.53 2.95 No function *5 African population 5.52 5,978 No function *106 African population 0 1.165 Undefined function *1 Americans 51.05 47,598 Normal function *2 Americans 22.09 19,219 Normal function *3 Americans 0.02 0.601 No function *4 Americans 10.25 9.159 No function *6 Americans 0.25 0.3 No function *7 Americans 0.5 0 No function *9 Americans 0.45 1.351 Decreased function *10 Americans 1.44 1,502 Decreased function *12 Americans 1.7 0 No function *17 Americans 0.48 0.901 Decreased function *28 Americans 0.09 0.15 Undefined function *29 Americans 0.19 0.3 Decreased function *33 Americans 0.17 0.3 Normal function *35 Americans 0.97 2,553 Normal function *41 Americans 2.33 6,306 Decreased function *82 Americans 2.5 0.45 Unknown function *1x2 Americans 2.86 1.201 Enhanced function *2x2 Americans 0.61 0.601 Enhanced function *4x2 Americans 0.11 0.15 No function *35x2 Americans 0 0.15 Enhanced function *5 Americans 1.59 2.102 No function *4+*68 Americans 0 2,553 No function *1 East Asians 24.74 27,216 Normal function *2 East Asians 12.09 8,041 Normal function *4 East Asians 0.54 0.206 No function *10 East Asians 43.56 15,464 Decreased function *14 East Asians 0.29 0.928 Decreased function *21 East Asians 0.35 0.515 No function *34 East Asians 1.02 0 Normal function *39 East Asians 0.59 0 Normal function *41 East Asians 2.27 3,918 Decreased function *49 East Asians 1.05 0.515 Decreased function *52 East Asians 0.18 0.103 Undefined function *65 East Asians 2.95 0 Undefined function *69 East Asians 1.17 0 No function *71 East Asians 0.12 0.619 Undefined function *1x2 East Asians 0.34 0.103 Enhanced function *2x2 East Asians 0.45 0.619 Enhanced function *10x2 East Asians 0.61 0.309 Decreased function *5 East Asians 4.84 3,608 No function *10+*36 East Asians 1.24 35,979 Decreased function *10+*36+*36 East Asians 0.45 1,649 Decreased function *1 European population 23.79 36,593 Normal function *2 European population 18.52 16.23 Normal function *3 European population 1.58 1,714 No function *4 European population 18.56 11,895 No function *6 European population 1.11 1.915 No function *9 European population 2.75 2.52 Decreased function *10 European population 1.58 1.411 Decreased function *17 European population 0.36 0.202 Decreased function *28 European population 0 0.504 Undefined function *31 European population 0.12 0.202 No function *33 European population 1.9 0.605 Normal function *34 European population 5.54 0 Normal function *35 European population 4.64 4.839 Normal function *39 European population 3.16 0 Normal function *41 European population 9.23 9.173 Decreased function *43 European population 0 0.101 Undefined function *59 European population 0.65 0.202 Decreased function *1x2 European population 0.83 0.504 Enhanced function *2x2 European population 0.84 1.512 Enhanced function *4x2 European population 0.66 0.302 No function *9x2 European population 0.01 0.101 Normal function *5 European population 2.99 2,419 No function *4+*68 European population 0 5,444 No function *4+*4N European population 0 0.706 No function *1 South Asians 24.9 40.147 Normal function *2 South Asians 29.3 20,964 Normal function *3 South Asians 0.11 0.21 No function *4 South Asians 9.13 8,281 No function *6 South Asians 0 0.105 No function *7 South Asians 0.41 0.943 No function *10 South Asians 8.84 3.878 Decreased function *29 South Asians 6.08 0 Decreased function *35 South Asians 1.1 0.629 Normal function *39 South Asians 0.2 0.21 Normal function *41 South Asians 12.29 12,159 Decreased function *1x2 South Asians 0.56 0.419 Enhanced function *2x2 South Asians 0.95 0.419 Enhanced function *5 South Asians 4.68 2.621 No function *10+*36 South Asians 0 1.153 Decreased function *4+*68 South Asians 0 2,201 No function *86 South Asians 0 2,306 Unknown function *111 South Asians 0 0.839 Unknown function *113 South Asians 0 0.839 Unknown function

Discussion

[0275] This example describes Cyrius, a method that allows the complex region of CYP2D6 to be accurately diplotyped. A unique feature of this example is that long reads of the data were used to confirm both haplotypes and SVs. Long reads provide a unique opportunity to confirm breakpoint regions of common SVs ( CYP2D6 deletions and duplications, and CYP2D6 /7 fusion genes) and confirm CYP2D6 gene phasing. Using 188 samples, including 50 with long-read validation data, as an orthogonal validation data set, Cyrius was shown to outperform other CYP2D6 genotypes, achieving 97.9% accuracy compared to 88.8% for Aldy and 85.6% for Stargazer. In particular, compared to these existing CYP2D6 identifiers, Cyrius enabled reads in regions where CYP2D6/7 have high similarity. Ambiguous read alignments in these regions can lead to incorrect copy number estimates and errors in identifying small variants. By accounting for possible unaligned reads and selecting a set of robust CYP2D6 /7 differentiation sites, Cyrius is able to detect star alleles much better using SV, achieving 97.2% accuracy compared to 88.9% for Aldy and 77.8% for Stargazer.

[0276] A total of 41 different star alleles were confirmed across all 188 validation samples, which represented (38.7%) of all star alleles listed in PharmGKB, including 53.4% of alleles with known functional status. Although based on analysis of the 1kGP samples in this example, the validation set included only 38.7% of all known star alleles, they were estimated to represent 96.5% of star alleles in the pangenomic population. Generally, allele frequencies calculated for 2504 1kGP samples from five ethnic populations were consistent with previous studies of simple star alleles. Conversely, some star alleles that were driven by the presence of SV showed very different frequencies, likely because many of the star alleles affected by SV are difficult to detect using conventional assays. This highlights the inherent pitfalls of pooling results from studies that have used many different CYP2D6 assays, some of which may be designed to simply detect a subset of star alleles. For example, of the 5 assays used to generate consensus GeT-RM genotypes, individual accuracies ranged from 47.1% to 75.2% relative to consensus (Table 25). A single method that can identify all known star alleles in a single analysis is the best choice for creating a population-level database.

Table 25: Accuracy of individual GeT-RM assays

TaqMan
+CNV
+XL-PCR PharmacoScan
v.r6 PharmacoScan
custom
v.r6+20180103 iPLEX
CYP2D6 V1.1 iPLEX
CYP2D6
V1.1
+individual panel and VeriDose Samples not listed as n/a 60.9% 78.2% 78.2% 100.0% 29.1% Accuracy 75.2% 47.1% 60.7% 59.8% 69.2%

[0277] Additionally, Cyrius was used to analyze 2504 1kGP samples from five ethnic populations to determine star allele frequencies. The calculated allele frequencies are consistent with previous studies of simple star alleles, and Cyrius has significantly improved estimates of star allele frequencies involving structural variants that may be difficult to detect using conventional methods.

[0278] Some existing methods rely on precise read alignment to distinguish between CYP2D6 and CYP2D7 , which may be error-prone due to several regions of high sequence similarity between the two genes, particularly between intron 1-exon 2 and exon 7-exon 9 ,Ambiguous alignment can introduce noise in depth profiles, ,leading to false CNV detections. Additionally, erroneous read alignments can lead to false positive or false negative variant determinations. In contrast, Cyrius first determined the total copy number of CYP2D6 + CYP2D7 by counting all reads that match either gene, and a total copy number of less than 4 clearly indicates the presence of SV. Not all differences based on the reference genome were used to determine the exact position of the SV. Many CYP2D6 / CYP2D7 base differences are not fixed, so not all of these positions can be used to reliably distinguish CYP2D6 from CYP2D7 (FIG. 26). Cyrius used 118 CYP2D6 / CYP2D7 differentiation positions selected to determine the exact position of the SV. By first determining the total copy number and then differentiating them using a subset of suitable differentiating bases, Cyrius was able to achieve more accurate SV determinations. To identify small variants, Cyrius overcomes the dependence on single-unit alignments by searching for variant reads at both the CYP2D6 position and the corresponding CYP2D7 position, thereby obtaining the most accurate small variant definitions.

[0279] In the example, long reads of the data were used to confirm both haplotypes and SV definitions. The PacBio data in this example provide a clear picture of the CYP2D6 - CYP2D7 region with high quality long reads (10 - 20 kb). In particular, PacBio data allows the identification of breakpoint regions for common structural variants ( CYP2D6 deletions and duplications, and CYP2D6 - CYP2D7 fusion genes). Even with PacBio reads, genotyping CYP2D6 can be challenging and may require targeted analysis, especially for structural variants involving duplication ( CYP2D6 and CYP2D6 - CYP2D7 duplication ) where the duplicated region contains >10 kb. For example, the de novo assembly approach failed to capture the *68 fusion in sample HG00733 (Figures 31A and 31B). Additionally, PacBio reads are not long enough to cover more than one copy of the duplicated sequence, and PacBio reads are too long to detect copy number in read counting (for short reads), so copy number calculation is difficult. Whole genome short-read sequencing provides the most accurate solution for CYP2D6 genotyping.

[0280] When analyzing 1kGP samples, Cyrius was able to determine the final genotype of over 97.6% of samples. In some embodiments, Cyrius was able to detect the remaining 2.4% of samples. For example, in samples where multiple haplotype configurations are possible, it may be useful to use a probabilistic approach to obtain the most likely genotype given the observed variations. In addition, continued sequencing and testing of more samples will help confirm the ability to genotype rare star alleles, as well as identify new variants that can be used to distinguish ambiguous diplotypes. This process was demonstrated in this example, where improvements were made to better identify three star alleles that were initially misidentified in 188 validation samples. The improvements were useful for genotyping at the population level, as three star alleles were found in almost 1% (23 of 2504) of 1kGP samples.

[0281] As new star alleles are identified, new star alleles may be added to the Cyrius database. One consideration when adding new star alleles that are defined by new variants is that the variants are unlikely to be considered in previous star allele definitions. As a result, there may be novel combinations of new and existing variants that cannot match any known combination, resulting in a lack of definitions. For example, Cyrius includes the ability to genotype for 25 new star alleles added in PharmVar v4 (not included in GeT-RM, Aldy or Stargazer). However, five (*119, *122, *135, *136, *139) of the 25 new star alleles have new variants that, when included, resulted in missing definitions in samples that could have been previously defined, indicating the existence of common new star alleles with a combination of variants not reflected in PharmVar. As a result, these five star alleles were removed along with two others (*127, with a gene conversion variant in the homologous region and *131 variants at the interference site) leaving the remaining 18. New star alleles may be possible as new variants/stars are identified alleles. Publicly available WGS datasets, such as the 2504 1kGP samples analyzed here, can be an important component of integrating new variants into star allele definitions, as these data allow for the rapid assessment of variants in many samples with different genotypes.

[0282] WGS provides a valuable opportunity to profile all genetic variations throughout the genome, but many of the clinically important regions/variants are beyond the capabilities of most secondary analysis pipelines. CYP2D6 is among complex genomic regions that are both clinically important and also require targeted bioinformatics solutions in addition to normal WGS pipelines. Such targeted methods have already been successfully applied to some challenging areas, such as the SMN1 gene responsible for spinal muscular atrophy, as shown in Example 1. More targeted methods, such as Cyrius, can accelerate the development of pharmacogenomics and enable personalized medicine.

Additional Considerations

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

[0284] One skilled in the art will appreciate that for this and other processes and methods described herein, the functions performed in the processes and methods may be implemented in a different order. In addition, the described steps and operations are provided as examples only, and some steps and operations may be optional, combined into fewer steps and operations, or expanded into additional steps and operations without detracting from the spirit of the described embodiments.

[0285] With respect to the application of substantially any plural and/or singular terms herein, those skilled in the art may change plural to singular and/or singular to plural as the context and/or scope of application requires. Throughout this document, various singular/plural combinations may be explicitly stated for clarity. As used in this specification and the appended claims, the singular forms include references to the plural unless the context clearly indicates otherwise. Accordingly, phrases such as “device configured for” are intended to include one or more of the aforementioned devices. Such one or more of the above-mentioned devices may also be collectively configured to perform the above-mentioned enumerations. For example, "a processor configured to perform statements A, B, and C may include a first processor configured to perform statements A operating in conjunction with a second processor configured to perform statements B and C. Any reference to "or ” as used herein implies the inclusion of “and/or” unless otherwise noted.

[0286] Those skilled in the art will appreciate that, in general, terms used herein, and particularly in the appended claims (e.g., the body of the appended claims), are essentially considered to be "non-limiting" terms ( e.g. , the term "including " should be interpreted as "including without limitation", the term "having" should be interpreted as "having at least", the term "includes" should be interpreted as "includes without limitation", etc.). In addition, those skilled in the art will appreciate that if reference to a presented claim is specifically intended, then such intent will be clearly stated in the claim, and in the absence of such reference, there is no such intent. For example, to facilitate understanding, the following appended claims may contain the introductory phrases “at least one” and “one or more” to introduce the enumeration of claims. However, the use of such phrases should not imply that the introduction of an enumeration in a claim using the singular form limits any particular claim containing such introduced enumeration in the claim to embodiments containing only one such enumeration, even if the same claim the invention includes the introductory phrases "one or more" or "at least one" and indefinite articles such as "some" ( for example , the singular form should be interpreted as "at least one" or "one or more"); the same is true with respect to the use of definite articles used to introduce the enumeration of the claims. Moreover, even if a specific number of claims presented is explicitly stated, those skilled in the art will understand that such listing should be interpreted to mean at least the stated amount ( for example, simply quoting the phrase "two statements" without other modifiers means at least two statements or two or more statements). In addition, in cases where a rule similar to the rule "at least one of A, B, C, etc." is used, such a construction is generally intended for one skilled in the art to understand the rule ( eg , "A system containing at least one of A, B, and C" will include, without limitation, systems that have only A, only B, only C, both A and B, both A and C, both B and C, and/or both A , B and C, etc.). In cases where a rule similar to "at least one of A, B, C, etc." is used, such a construction is generally intended for one skilled in the art to understand the rule ( e.g. , "A system containing at least one of A, B or C" will include, without limitation, systems that have only A, only B, only C, both A and B, both A and C, both B and C, and/or both A, B and C, etc.). In addition, those skilled in the art will appreciate that virtually any dividing word and/or phrase presenting two or more alternative terms, whether in the description, claims, or drawings, should be understood to imply the possibility of including one of the terms, any of terms or both terms. For example, the phrase "A or B" should be understood to include the possibilities "A" or "B" or "A and B."

[0287] Moreover, if features or aspects of the disclosure are described in terms of Markush groups, those skilled in the art will appreciate that the description is thus also described in terms of any individual member or subset of members of the Markush group.

[0288] As one skilled in the art will appreciate, for any and all purposes, such as in terms of providing a written description, all ranges described herein also include all possible sub-ranges and combinations of sub-bands thereof. Any range listed can easily be recognized as sufficiently descriptive to allow 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 described herein can be easily broken down into a lower third, a middle third, and an upper third, etc. As will also be appreciated by one skilled in the art, all expressions such as “up to”, “at least”, “greater than”, “less than”, etc. include the number listed and refer to ranges, which can subsequently be broken down into sub-bands as discussed above. Finally, as one skilled in the art will appreciate, the range includes each individual element. Thus, for example, a group having 1-3 points refers to groups having 1, 2 or 3 points. Similarly, a group containing 1-5 items refers to groups containing 1, 2, 3, 4 or 5 items, etc.

[0289] It should be understood that various embodiments of this specification are described herein for purposes of illustration and that various modifications may be made without departing from the scope and spirit of this specification. Accordingly, the various embodiments described herein are not intended to be limiting, and the scope and spirit are set forth in the following claims.

[0290] It should be understood that not all objects or benefits may be achieved in accordance with any specific embodiment described herein. Thus, for example, those skilled in the art will appreciate that certain embodiments may be configured to operate in a manner that provides or optimizes one benefit or group of benefits described herein without the need to achieve other objects or benefits that may be described or suggested in this document.

[0291] All processes described herein can be implemented and fully automated using program code modules executed by a computer system that includes one or more computers or processors. The program modules may be stored on any type of non-volatile computer readable medium or other computer storage device. Some or all of the methods can be implemented on specialized computer hardware.

[0292] Many other variations than those described herein will be apparent from this description. For example, depending on the 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 removed altogether (for example, not all actions or events described are necessary for the practical implementation of the algorithms ). Moreover, in certain embodiments, actions or events may be performed simultaneously, for example, through multi-threaded processing, interrupt processing, or multiple processors or processor cores, or other parallel architectures, rather than sequentially. In addition, different tasks or processes may be performed by different machines and/or computing systems that may operate together.

[0293] Various illustrative logic blocks and modules described in connection with the embodiments disclosed herein may be implemented or executed by a machine, such as a data processing unit or digital signal processor (DSP), application specific integrated circuit (ASIC) , field programmable gate array (FPGA) or other programmable logic device, discrete logic gate or transistor logic circuit, discrete hardware component, or any combination thereof, designed to perform the functions described in this document. The processor may be a microprocessor, but in an alternative embodiment, the processor may be a controller, microcontroller, or state machine, combinations thereof, or the like. The processor may include electrical circuitry configured to process computer-executable instructions. In another embodiment, the processor includes an FPGA or other programmable device that performs logical operations without processing computer-executable instructions. The processor may also be implemented as a combination of computing devices, such as a combination of a DSP and a microprocessor, a plurality of microprocessors, one or more microprocessors coupled to a DSP core, or any other such configuration. Although described herein primarily in terms of digital technology, 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 mixed analog and digital circuitry. The computing environment may include any type of computer system, including, without limitation, a computer system based on a microprocessor, a system computer, a digital signal processor, a portable computing device, a device controller, or an in-device computing engine, and many others.

[0294] Any process descriptions, elements, or blocks in the flowcharts described herein and/or shown in the accompanying figures should be understood to potentially represent modules, segments, or portions of code that include one or more executable instructions for implementing specific logic functions or elements in a process. Alternative embodiments are included within the scope of the embodiments described herein, in which elements or functions may be removed, performed in an order different from those shown or described, including substantially simultaneously or in reverse order, depending on the functionality involved, as will understandable to those skilled in the art.

[0295] It should be emphasized that many variations and modifications can be made to the embodiments described above, elements of which should be understood to refer to other suitable examples. All such modifications and variations are intended to be included within the scope of this specification and are protected by the following claims.

Claims (63)

1. A method for determining the copy number of the survival motor neuron 1 ( SMN1 ) gene, comprising: controlled by hardware processor: obtaining sequencing data containing multiple sequence reads obtained from a sample from a subject aligned to the survival motor neuron gene 1 ( SMN1 ) or survival motor neuron gene 2 ( SMN2 ); determining (i) a first number of sequence reads from the plurality of sequence reads aligned to a first region of SMN1 or SMN2 containing at least 1 of the first 6 exons of the SMN1 gene or the SMN2 gene, respectively, and (ii) a second number of sequence reads from the plurality of sequence reads, aligned to a second region of SMN1 or SMN2 containing at least exon 7 or exon 8 of the SMN1 gene or the SMN2 gene, respectively; determining (i) a first normalized number of sequence reads aligned to a first region of the SMN1 or SMN2 gene, and (ii) a second normalized number of sequence reads aligned to a second region of SMN1 or SMN2 , using (i) the length of the first SMN1 or SMN2 region and ( ii) the length of the second region SMN1 or SMN2 , respectively; determining (i) the copy number of the total number of survival motor neuron (SMN) genes, each of which is an intact SMN1 gene, an intact SMN2 gene, a truncated SMN1 gene, or a truncated SMN2 gene; and (ii) the copy number of any intact SMN genes, each of which is an intact SMN1 gene or an intact SMN2 gene, using a mixture normal distribution model containing a plurality of normal distributions, each representing a different integer copy number, given (i) the first normalized the 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; for one set of SMN1 gene-specific bases associated with an intact SMN1 gene, determining the most likely combination from a set of possible combinations, each containing a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene, which add up to the number of copies of any specific intact genes SMN taking into account (a) the number of sequence reads from the set of sequence reads with bases that contain a base specific to the SMN1 gene, and (b) the number of sequence reads from the set of sequence reads with bases that contain the SMN2 gene-specific base of the SMN2 gene corresponding a base specific to the SMN1 gene; And determining the copy number of the SMN1 gene using the most likely combination of the possible copy number of the SMN1 gene and the possible number of copies of the SMN2 gene, determined for a specific base for the SMN1 gene. 2. The method of claim 1, wherein the sequencing data comprises whole genome sequencing (WGS) data or short-read WGS data. 3. Method according to paragraphs. 1, 2, wherein the subject is a fetal subject, a neonatal subject, a pediatric subject, an adolescent subject, or an adult subject. 4. Method according to any one of paragraphs. 1-3, in which the sample contains cells or extracellular DNA. 5. Method according to any one of paragraphs. 1-4, in which the sample contains fetal cells or extracellular fetal DNA. 6. Method according to any one of paragraphs. 1-5, wherein a sequence read from a plurality of sequence reads is aligned to a first SMN1 or SMN2 region or a second SMN1 or SMN2 region with an alignment quality score of approximately zero. 7. Method according to any one of paragraphs. 1-6, in which the first region SMN1 or SMN2 contains exon 1 to 6 of the SMN1 gene or the SMN2 gene, respectively, and has a length of approximately 22.2 kb, and the second region SMN1 or SMN2 contains exon 7 and exon 8 of the gene SMN1 or SMN2 gene, respectively, and is approximately 6 kb in length. 8. Method according to any one of paragraphs. 1-7, wherein determining (i) a first normalized number of sequence reads aligned to a first region SMN1 or SMN2 and (ii) a second normalized number of sequence reads aligned to a second region comprises: determining (i) a first normalized number of sequence reads , aligned to the first SMN1 or SMN2 region, and (ii) a second normalized number of sequence reads aligned to the 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, and (iii) sequence read depths of the subject's genome region, excluding genetic loci containing the SMN1 gene and the SMN2 gene, as determined by sequencing data. 9. The method of claim 8, wherein determining (i) a first normalized number of sequence reads aligned to the first SMN1 or SMN2 region, and (ii) a second normalized number of sequence reads aligned to the second SMN1 or SMN2 region, comprises: determining (i) the length-normalized number of sequence reads aligned to the first SMN1 or SMN2 region, and (ii ) the length-normalized number of sequence reads aligned to the second SMN1 or SMN2 region, using ( i ) the length of the first region SMN1 or SMN2 and (ii) the length of the second region SMN1 or SMN2 , respectively; And determining (i) a first normalized read depth of sequences aligned to the first SMN1 or SMN2 region, and (ii) a second normalized read depth of sequences aligned to a second SMN1 or SMN2 region, based on (i) a first length-normalized count of the SMN1 or SMN2 region and (ii) a second length-normalized count of the SMN1 or SMN2 region, respectively, using the depth of sequence reads of a region of the subject's genome other than the genetic loci containing the SMN1 gene and the SMN2 gene; a first normalized number of sequence reads aligned to the first SMN1 or SMN2 region, and a second normalized number of sequence reads aligned to the second SMN1 or SMN2 region, which represent the first normalized depth and the second normalized depth, respectively. 10. Method according to any one of paragraphs. 1-9, wherein determining (i) a first normalized number of sequence reads aligned to a first region SMN1 or SMN2 and (ii) a second normalized number of sequence reads aligned to a second region comprises: determining (i) a first normalized number of sequence reads , aligned to the first SMN1 or SMN2 region, and (ii) a second normalized number of sequence reads aligned to the second SMN1 or SMN2 region, using (i) the GC content of the first SMN1 or SMN2 region and (ii) the GC content of the second SMN1 region or SMN2 , respectively, and (iii) the depth of sequence reads of a region of the subject's genome other than the genetic loci containing the SMN1 gene and the SMN2 gene as determined by sequencing and (iv) the GC content of the genomic region. 11. Method according to any one of paragraphs. 8 - 10, wherein the read depth of the region includes the average depth or median read depth of the sequence of a region of the subject's genome other than the genetic loci containing the SMN1 gene and the SMN2 gene, as determined by sequencing. 12. The method of claim 11, wherein the region includes approximately 3000 preselected regions of approximately 2 kb in length. each in the subject's genome. 13. Method according to any one of paragraphs. 1-12, wherein (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 . 14. Method according to any one of paragraphs. 1-13, in which the normal distribution model includes a one-dimensional normal distribution model. 15. Method according to any one of paragraphs. 1-14, in which the set of normal distributions of the mixture of normal distributions model represents an integer number of copies from 0 to 10. 16. Method according to any one of paragraphs. 1-15, wherein the mean for each of a plurality of normal distributions is the integer number of copies represented by the normal distribution. 17. Method according to any one of paragraphs. 1-16, wherein determining (i) the copy number of all SMN genes and (ii) the copy number of any intact SMN genes includes determining (i) the copy number of all SMN genes and (ii) the copy number of any intact SMN genes using a normal mixture model distributions, and a first predetermined posterior probability threshold given (i) a first normalized number of sequence reads aligned to the first SMN1 or SMN2 region, and (ii) a second normalized number of sequence reads aligned to the second SMN1 or SMN2 region, respectively. 18. The method of claim 17, wherein the first determined posterior probability threshold is 0.95. 19. Method according to any one of paragraphs. 1-18, comprising determining the copy number of truncated SMN genes using (i) a certain copy number of total SMN genes and (ii) a certain copy number of intact SMN genes. 20. The method of claim 19, wherein the copy number of the truncated SMN genes is the difference between (i) the determined total copy number of the SMN genes and (ii) the determined copy number of the intact SMN genes. 21. Method according to any one of paragraphs. 1-20, in which the SMN1 gene-specific base is a splicing enhancer. 22. Method according to any one of paragraphs. 1-21, wherein the SMN1 gene-specific base is the base at c.840 of the SMN1 gene. 23. Method according to any one of paragraphs. 1-22, in which the most likely combination of a possible SMN1 gene copy number and a possible SMN2 gene copy number is associated with the highest posterior probability, compared to other combinations of a plurality of combinations given (a) the number of sequence reads from the plurality of sequence reads with bases that contain an SMN1 gene-specific base, and (b) the number of sequence reads from a plurality of sequence reads with bases that contain the corresponding SMN2 gene-specific base. 24. Method according to any one of paragraphs. 1-23, in which determining the most likely combination of a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene includes: determining the most likely combination from a plurality of possible combinations, each of which contains a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene, which give in total, a certain number of copies of any intact SMN genes, taking into account the ratio of (a) the number of sequence reads from the set of sequence reads with bases that contain the SMN1 gene-specific base, and (b) the number of sequence reads from the set of sequence reads with bases that contain the specific for the SMN2 gene, the base of the SMN2 gene corresponding to the base specific for the SMN1 gene. 25. Method according to any one of paragraphs. 1-24, in which determining the most likely combination of a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene includes: determining (a) the number of sequence reads from the plurality of sequence reads with bases that contain the SMN1 gene-specific base, and (b) the number of sequence reads from the plurality of sequence reads with bases that contain the SMN2 gene-specific base corresponding to the gene-specific SMN1 base; determining the ratio of (a) the number of sequence reads from a plurality of sequence reads with bases that contain an SMN1 gene-specific base, and (b) the number of sequence reads from a plurality of sequence reads with bases that contain an SMN2 gene-specific base corresponding to the SMN2 gene-specific SMN1 gene base; And determining the most likely combination from a set of possible combinations, each containing a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene, which add up to a certain number of copies of any intact SMN gene based on the ratio of (a) the number of sequence reads from the set of sequence reads with bases that contain the SMN1 gene-specific base, and (b) the number of sequence reads from the plurality of sequence reads with bases that contain the SMN2 gene-specific base of the SMN2 gene corresponding to the SMN1 gene-specific base. 26. Method according to any one of paragraphs. 1-25, wherein determining the most likely combination of a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene includes: for each of a plurality of SMN1 gene-specific bases, determining the most likely combination from a plurality of possible combinations, each of which contains a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene that add up to a certain number of copies of any intact SMN genes that is associated with the largest posterior probability given (a) the number of sequence reads from the set of sequence reads with bases that contain the SMN1 gene-specific base, and (b) the number sequence reads from a plurality of sequence reads with bases that contain an SMN2 gene-specific SMN2 gene base corresponding to an SMN1 gene-specific base, 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, 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 a plurality of SMN1 gene-specific bases. 27. The method of claim 26, wherein the SMN1 gene-specific base has a match with each of the plurality of SMN1 gene-specific bases other than the SMN1 gene-specific base above a predetermined match threshold. 28. The method according to claim 27, in which the compliance threshold is 97%. 29. Method according to any one of paragraphs. 26-28, in which the plurality of SMN1 gene-specific bases includes 8 SMN1 gene-specific bases. 30. Method according to any one of paragraphs. 26-29, wherein each of a plurality of SMN1 gene-specific bases is located in intron 6, exon 7, intron 7, or exon 8 of the SMN1 gene. 31. Method according to any one of paragraphs. 26-30, wherein a plurality of bases specific to the SMN1 gene if the subject belongs to a first race, a plurality of bases specific to the SMN1 gene if the subject belongs to a second race, and a plurality of bases specific to the SMN1 gene if the subject belongs to an unknown race, differ. 32. Method according to any one of paragraphs. 26-31, in which the race of the subject is unknown and yet many bases specific to the SMN1 gene are not race specific. 33. Method according to any one of paragraphs. 26-31, wherein the race of the subject is known and wherein a plurality of bases specific to the SMN1 gene are specific to the race of the subject. 34. Method according to any one of paragraphs. 26-33, additionally including: obtaining information about the subject's race; And based on the obtained race information, a plurality of bases specific for the SMN1 gene are selected from the plurality of bases specific for the SMN1 gene. 35. Method according to any one of paragraphs. 1-34, in which 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 number of copies of the SMN1 gene and the possible number of copies of the SMN2 gene determined for each of a plurality of bases specific for the gene SMN1 . 36. Method according to any one of paragraphs. 1-35, wherein determining the SMN1 gene copy number includes: determining the SMN1 gene copy number using the most likely combination of a possible SMN1 gene copy number and a possible SMN2 gene copy number determined for a base specific to the SMN1 gene, and a second predetermined posterior threshold probabilities for the combination of the possible number of copies of the SMN1 gene and the possible number of copies of the SMN2 gene. 37. The method of claim 36, wherein the second predetermined posterior probability threshold is 0.6 or 0.8. 38. Method according to any one of paragraphs. 26-37, in which the majority of possible SMN1 gene copy number values are agreed upon, and wherein the determined SMN1 gene copy number represents the agreed upon possible SMN1 gene copy number. 39. The method according to claim 38, including: determining a possible combination containing a possible number of copies of the SMN1 gene and a possible number of copies of the SMN2 gene that add up to a certain number of copies of any intact SMN gene, taking into account (a) the number of sequence reads from a plurality of sequence reads with bases that contain any of the plurality of bases, specific for the SMN1 gene, and (b) the number of sequence reads from a plurality of sequence reads with bases that contain any of a plurality of corresponding bases specific for the SMN2 gene; And determining the possible copy number of a possible combination is the consensus possible copy number of the SMN1 gene. 40. Method according to any one of paragraphs. 1-39, wherein determining the copy number of the SMN1 gene includes determining the copy number of the SMN1 gene to be zero, one, or greater than one. 41. Method according to any one of paragraphs. 1-40, including determining the spinal muscular atrophy (SMA) status of a subject based on the copy number of the SMN1 gene. 42. The method of claim 41, wherein the subject's SMA status includes SMA, SMA carrier/non-SMA carrier, and non-SMA carrier. 43. Method according to any one of paragraphs. 1-42, comprising identifying a subject as a silent carrier of SMA using a set of sequence reads from a plurality of sequence reads aligned to g.27134 of the SMN1 gene and base sequence reads aligned to g.27134 of the SMN1 gene. 44. Method according to one of paragraphs. 1-43, including determining a treatment recommendation for a subject based on a certain number of copies of the SMN1 gene. 45. The method of claim 44, wherein the treatment recommendation includes administering Nusinersen and/or Zolgensma to the subject.