A method and apparatus for predicting a target gene copy number type

Abstract

A method and device for predicting a genotype of a target gene, the method comprising: a baseline correction step; a feature extraction step, comprising obtaining a cluster center of a specific region of the target gene according to a corrected depth of each capture site, and taking the cluster center as a feature value of the corresponding region; and a prediction step, comprising predicting a genotype of the target gene in sequencing data of a sample to be tested according to the feature value. The method has the advantages of fast analysis speed, high accuracy, wide application range, no need for reference set correction, avoidance of the influence of sequencing fluctuation on the analysis result, ease of implementation of high-throughput screening typing at a population level, and applicability to individualized medication guidance for more populations and acceleration of large-scale research and development of pharmacogenomics.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of bioinformatics, and particularly relates to a method and device for predicting a target gene copy number type. BACKGROUND

[0002] Cytochrome P450 isoenzyme 2D6 (CYP2D6) is a liver enzyme of great importance to pharmacogenetics, accounting for 2% of total liver CYP450 enzyme protein, and is involved in the metabolism of up to 25% of clinically commonly used drugs. The gene encoding the CYP2D6 metabolic enzyme is CYP2D6, located on chromosome 22q13.1, and is part of the cytochrome P450 gene family, which encodes metabolic enzymes involved in phase I metabolism and clearance of many endogenous substrates (such as steroids, fatty acids, biogenic amines, etc.) and a variety of commonly used drugs (antidepressants, antipsychotics, antiarrhythmic drugs, opioid drugs, and beta-receptor blockers, etc.). Moreover, among the cytochrome P450 gene family, CYP2D6 is the only non-inducible enzyme, which leads to a very high correlation between CYP2D6 gene variation and individual differences in enzyme activity. According to the ability to metabolize CYP2D6 substrates (i.e., enzyme activity), subjects can be divided into four categories, from strong to weak enzyme activity: ultrarapid metabolizer (UM), extensive metabolizer (EM), intermediate metabolizer (IM), and poor metabolizer (PM). The CYP2D6 gene is highly polymorphic, with more than 130 reported allelic variations (https: / / www.pharmvar.org / gene / CYP2D6), including point mutations, whole gene deletions, and complex structural variations. Due to the presence of abundant polymorphic sites and complex site combination forms, PharmVar has standardized the nomenclature of gene polymorphisms in the cytochrome P450 gene family, including the star nomenclature of the CYP2D6 gene, such as CYP2D6*3, CYP2D6*4, CYP2D6*5, CYP2D6*10, etc.

[0003] CYP2D6*5 represents a whole gene deletion, which is a copy number variation, corresponding to a non-functional metabolic enzyme, and is the second largest CYP2D6 polymorphic variation in East Asian population, with a frequency as high as 5%. If the CYP2D6 genotype of a subject is *5 / *5, the corresponding CYP2D6 enzyme activity is almost completely lost (i.e. PM metabolic type), which will lead to a large accumulation of drugs metabolized by CYP2D6 enzymes in the body, and the concentration of active products is too low, resulting in the risk of serious adverse reactions or insufficient efficacy. For example, codeine needs to be metabolized by CYP2D6 enzyme to active metabolite morphine to exert analgesic effect, and CYP2D6 PM patients have poor metabolic capacity for codeine, which cannot make the morphine in the blood reach the target concentration, thereby leading to insufficient efficacy, and also causing codeine to accumulate in the body, resulting in the risk of adverse reactions. This insufficient efficacy cannot be improved by increasing the dose of the drug, and other drugs that are not metabolized by CYP2D6 enzymes need to be considered. Therefore, accurate genotyping of CYP2D6 gene copy number has very important guiding significance for clinical drug use.

[0004] Due to the presence of pseudogenes CYP2D8 (92% homologous) and CYP2D7 (97% homologous) of CYP2D6, and the fusion gene of CYP2D6 and CYP2D7, the accurate detection of CYP2D6 gene copy number becomes very complex. There is no approved kit for CYP2D6*5 detection on the market. However, there are some kits or detection techniques for research, such as PCR electrophoresis method (long fragment PCR + agarose gel electrophoresis method), fluorescent quantitative PCR method (such as: ThermoFisher CYP2D6 TaqMan MGB probe), mass spectrometry method (such as: Agena Veridose CYP2D6), comparative genomic hybridization method (such as: Affymetrix DMET plus microarray), high-throughput sequencing method (Panel, WGS) and the like. The PCR electrophoresis method has high accuracy, but it needs to be followed by agarose gel electrophoresis after long-time PCR (about 3 hours), which is complicated, low throughput and time-consuming. Although the fluorescent quantitative PCR method is absolute quantification, it has poor repeatability. The mass spectrometry method and the comparative genomic hybridization method have complex experimental steps, high environmental requirements, great technical difficulty, high cost, and each experiment needs a normal copy number control sample as a reference set. SUMMARY

[0005] According to the first aspect, in an embodiment, a method for predicting a target gene copy number type is provided, comprising:

[0006] The baseline correction step includes using the median and / or average sequencing depth of each capture site of the target gene in the sequencing data of samples with known target gene copy number type as the target gene whole gene deletion type as the calibration baseline, calculating the baseline depth of each capture site of the target gene, and subtracting the baseline depth from the standardized depth of each capture site of the target gene in the sequencing data of the sample to be tested to obtain the corrected depth.

[0007] The feature extraction step includes obtaining the cluster center of the specific region of the target gene based on the corrected depth of each capture site, and using the cluster center as the feature value of the corresponding region.

[0008] The prediction step includes predicting the copy number type of the CYP2D6 gene in the sequencing data of the sample to be tested based on the feature value.

[0009] According to a second aspect, in one embodiment, an apparatus for predicting the copy number type of a target gene is provided, comprising:

[0010] The baseline correction module includes using the median and / or average sequencing depth of each capture site of the target gene in the sample sequencing data with the known target gene copy number type as the target gene whole gene deletion type as the calibration baseline, calculating the baseline depth of each capture site of the target gene, and subtracting the baseline depth from the standardized depth of each capture site of the target gene in the sequencing data of the sample to be tested to obtain the corrected depth;

[0011] The feature extraction module includes obtaining the cluster center of the specific region of the target gene based on the corrected depth of each capture site, and using the cluster center as the feature value of the corresponding region.

[0012] The prediction module includes predicting the copy number type of the target gene in the sequencing data of the sample to be tested based on the feature value.

[0013] According to a third aspect, in one embodiment, an apparatus is provided, comprising:

[0014] Memory, used to store programs;

[0015] A processor for implementing the method as described in the first aspect by executing a program stored in the memory.

[0016] According to a fourth aspect, in one embodiment, a computer-readable storage medium is provided, the medium storing a program that can be executed by a processor to implement the method as described in the first aspect.

[0017] According to the above embodiments, a method and apparatus for predicting the copy number type of a target gene are provided. This method has fast analysis speed, high accuracy, wide applicability, and does not require reference set correction, thus avoiding the impact of sequencing fluctuations on the analysis results. It is also easy to achieve high-throughput screening and typing at the population level, which can be used for personalized medication guidance for more populations and accelerate the large-scale research and development of pharmacogenomics. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the CYP2D6 gene copy number in one embodiment;

[0019] Figure 2 This is a flowchart of CYP2D6 gene copy number analysis in one embodiment;

[0020] Figure 3 This is an example diagram showing the normalized depth of different regions of the CYP2D6 gene in 0N, 1N, and 2N type samples according to one embodiment. Detailed Implementation

[0021] The present invention will now be described in further detail with reference to specific embodiments and accompanying drawings. Similar elements in different embodiments are referred to by associated similar element reference numerals. In the following embodiments, many details are described to facilitate a better understanding of this application. However, those skilled in the art will readily recognize that some features may be omitted in different situations, or may be replaced by other elements, materials, or methods. In some cases, certain operations related to this application are not shown or described in the specification. This is to avoid obscuring the core parts of this application with excessive description. For those skilled in the art, detailed description of these related operations is not necessary; they can fully understand the related operations based on the description in the specification and general technical knowledge in the art.

[0022] Furthermore, the features, operations, or characteristics described in the specification can be combined in any suitable manner to form various embodiments. At the same time, the steps or actions in the method description can be rearranged or adjusted in a manner obvious to those skilled in the art. Therefore, the various orders in the specification and drawings are only for the clear description of a particular embodiment and do not imply a necessary order, unless otherwise stated that a particular order must be followed.

[0023] The serial numbers assigned to components in this article, such as "first" and "second", are used only to distinguish the objects being described and have no sequential or technical meaning.

[0024] According to a first aspect, in one embodiment, a method for predicting the copy number type of a target gene is provided, comprising:

[0025] The baseline correction step includes using the median and / or average sequencing depth of each capture site of the target gene in the sequencing data of samples with known target gene copy number type as the target gene whole gene deletion type as the calibration baseline, calculating the baseline depth of each capture site of the target gene, and subtracting the baseline depth from the standardized depth of each capture site of the target gene in the sequencing data of the sample to be tested to obtain the corrected depth.

[0026] The feature extraction step includes obtaining the cluster center of the specific region of the target gene based on the corrected depth of each capture site, and using the cluster center as the feature value of the corresponding region.

[0027] The prediction step includes predicting the copy number type of the CYP2D6 gene in the sequencing data of the sample to be tested based on the feature value.

[0028] In one embodiment, during the baseline correction step, the target gene contains a pseudogene. The activity of the product expressed by the pseudogene is almost completely lost.

[0029] In one embodiment, the target gene in the baseline correction step includes, but is not limited to, the CYP2D6 gene.

[0030] In one embodiment, in the baseline correction step, the target gene whole gene deletion type means that both alleles of the target gene are whole gene deletion types, which can be represented by target gene *5 / *5.

[0031] In one embodiment, the target gene deletion type in the baseline correction step includes, but is not limited to, the CYP2D6 gene *5 / *5.

[0032] In one embodiment, the method of the present invention can be extended to other genes; for genes with pseudogenes, it is necessary to find regions that are different from pseudogenes, and the relationship between normalization depth and known copy number types can be analyzed to find specific regions that can reflect copy number, thereby enabling the prediction of copy number types.

[0033] In one embodiment, in the feature extraction step, the specific region includes at least one region among the regions of exons 1, 3, 5, and 6 and introns 2, 5, and 6 of the target gene.

[0034] In one embodiment, in the feature extraction step, the specific region includes all regions of exons 1, 3, 5, and 6 and introns 2, 5, and 6 of the target gene.

[0035] In one embodiment, the sequencing depth of each capture site in the baseline correction step can be a normalized depth. The purpose of normalization is to reduce data fluctuations caused by sequencing type, different experimental conditions, sample type, etc.; without normalization, the data will be affected by the aforementioned heterogeneous factors.

[0036] In one embodiment, in the baseline correction step, the median and average of the sequencing depth (preferably normalized depth) of each capture site can be used as the calibration baseline. However, the average is easily affected by outliers, while the median better reflects the location of baseline clustering and is not sensitive to outliers. Therefore, the median is preferred as the calibration baseline.

[0037] In one embodiment, the formula for calculating the standardized depth in the baseline correction step is as follows:

[0038] (normalized_depth) ij ) m×n =(depth) ij / mean_depth i ) m×n ;

[0039] Where, depth ij The mean_depth represents the depth of the j-th capture site of the target gene in the sequencing data of the i-th sample. i Normalized_depth represents the average depth of the i-th sample to be tested. ij denoted as the normalized depth of the j-th capture site of the target gene in the sequencing data of the i-th sample to be tested, m represents the number of samples, and n represents the length of the target gene. The step of calculating the normalized depth can eliminate the influence of different experimental conditions and sequencing data volume of samples within the batch.

[0040] In one embodiment, n is 4312.

[0041] In one embodiment, the baseline depth is calculated in the baseline correction step according to the following formula:

[0042] (baseline j ) 1×n =(median(normalized_depth) ·j )) 1×n ;

[0043] Among them, baseline j This represents the baseline depth of the j-th capture site of the target gene, specifically the median of the normalized depths of the j-th capture site across all 0N samples. 0N samples are target gene *5 / *5 type samples, meaning samples where both alleles are completely deleted.

[0044] In one embodiment, the corrected depth is calculated in the baseline correction step according to the following formula:

[0045] (fixed_depth ij ) m×n =(normalized_depth ij -baseline j ) m×n ;

[0046] Among them, fixed_depth ij The normalized_depth represents the baseline-corrected depth of the j-th capture site of the target gene in the sequencing data of the i-th sample. ij The baseline represents the normalized depth of the j-th capture site of the target gene in the sequencing data of the i-th sample. j This represents the baseline depth of the j-th capture site. The corrected depth is used for subsequent gene copy number analysis. If the target gene is the CYP2D6 gene, this step can eliminate the influence of factors such as CYP2D7 and CYP2D8 homologs, and different sequencing chip platforms between batches.

[0047] In one embodiment, the feature extraction step includes:

[0048] The steps for calculating the sum of squared distances are as follows: For at least one region among exons 1, 3, 5, 6 and introns 2, 5, 6 of the target gene, firstly, the baseline-corrected depth of a capture site within the region is randomly selected as the initial cluster center, and then the sum of squared distances between the baseline-corrected depths of other capture sites within the region and the cluster center is calculated.

[0049] Repeat the steps, moving the cluster center capture site, and repeat the calculation of the sum of squared distances according to the steps of calculating the sum of squared distances until the sum of squared distances reaches its minimum value. Use the cluster center at this point as the feature value of the region.

[0050] In one embodiment, in the step of calculating the sum of squares of distances, for each region of exons 1, 3, 5, 6 and introns 2, 5, 6 of the target gene, the baseline-corrected depth of a capture site within the region is first randomly selected as the initial cluster center, and then the sum of squares of the distances between the baseline-corrected depths of other capture sites within the region and the cluster center is calculated.

[0051] In one embodiment, the sequencing data is capture sequencing data that captures the target gene or whole genome sequencing (WGS).

[0052] In one embodiment, the sequencing data is sequencing data aligned to a reference genome sequence.

[0053] In one embodiment, the reference genome sequence is derived from a reference group.

[0054] In one embodiment, the reference genome sequence includes a common sequence from a reference group.

[0055] In one embodiment, the reference genome sequence includes, but is not limited to, at least a portion of the hg19 human genome, hg38 genome, hg18 genome, hg17 genome, or hg16 genome.

[0056] In one embodiment, the sequencing data is sequencing data generated based on next-generation sequencing (NGS).

[0057] In one embodiment, in the prediction step, the copy number type includes at least one of 0N, 1N, and 2N, where 0N represents the target gene *5 / *5 type, 1N represents the target gene *1 / *5 type, and 2N represents the target gene *1 / *1 type.

[0058] In one embodiment, the prediction step includes using a model to predict the copy number type of the target gene in the sequencing data of the sample to be tested, based on the feature value.

[0059] In one embodiment, the algorithm of the model includes, but is not limited to, ensemble learning algorithms.

[0060] In one embodiment, the ensemble learning algorithm includes, but is not limited to, at least one of Bagging-based algorithms and Boosting-based algorithms.

[0061] In one embodiment, the Bagging-based algorithm includes, but is not limited to, random forest.

[0062] In one embodiment, the Boosting-based algorithm includes, but is not limited to, at least one of Adaboost, GBDT, and XGBOOST.

[0063] According to a second aspect, in one embodiment, an apparatus for predicting the copy number type of a target gene is provided, comprising:

[0064] The baseline correction module includes using the median and / or average sequencing depth of each capture site of the target gene in the sample sequencing data with the known target gene copy number type as the target gene whole gene deletion type as the calibration baseline, calculating the baseline depth of each capture site of the target gene, and subtracting the baseline depth from the standardized depth of each capture site of the target gene in the sequencing data of the sample to be tested to obtain the corrected depth;

[0065] The feature extraction module includes obtaining the cluster center of the specific region of the target gene based on the corrected depth of each capture site, and using the cluster center as the feature value of the corresponding region.

[0066] The prediction module includes predicting the copy number type of the target gene in the sequencing data of the sample to be tested based on the feature value.

[0067] According to a third aspect, in one embodiment, an apparatus is provided, comprising:

[0068] Memory, used to store programs;

[0069] A processor for implementing the method as described in the first aspect by executing a program stored in the memory.

[0070] According to a fourth aspect, in one embodiment, a computer-readable storage medium is provided, the medium storing a program that can be executed by a processor to implement the method as described in the first aspect.

[0071] Currently available methods for detecting CYP2D6 gene copy number have certain shortcomings, such as the cumbersome and low-throughput nature of PCR electrophoresis, the high cost of mass spectrometry and comparative genomic hybridization, and the poor reproducibility of fluorescence PCR. With the popularization of sequencing technology, some analysis software based on high-throughput sequencing data for detecting CYP2D6 gene copy number has emerged; however, they all have limitations. For example, Stargazer v1.0.7 and Aldy v2.2.3 only support variant detection results based on the hg19 reference genome, and for captured sequencing data, they only support the PGRNseq capture sequencing panel. Although Cyrius v1.1.1 supports alignment results based on the hg19 / hg38 reference genome, it is only applicable to high-throughput sequencing data from WGS (Whole Genome Sequencing). To address the limitations of the aforementioned detection / analysis methods, in one embodiment, this invention provides a novel method for analyzing the copy number of the CYP2D6 gene based on high-throughput sequencing data. This method offers fast analysis speed, high accuracy, and wide applicability (supporting hg19 / hg38 reference genomes and suitable for panel and WGS sequencing data that have captured the CYP2D6 gene). Furthermore, it eliminates the need for reference set correction, avoiding the impact of sequencing fluctuations on the analysis results. It also facilitates high-throughput screening and genotyping at the population level, enabling personalized medication guidance for a wider range of populations and accelerating large-scale pharmacogenomics research and development.

[0072] In one embodiment, the present invention provides a method for using high-throughput sequencing data to achieve copy number analysis of the CYP2D6 gene (see schematic diagram). Figure 1This invention relates to a method and apparatus for accurate analysis of *5 / *5 types (abbreviated as "0N"), *1 / *5 types (abbreviated as "1N"), and *1 / *1 types (abbreviated as "2N"), applicable to panel and WGS sequencing data capturing the CYP2D6 gene. The apparatus only requires inputting the alignment result file (bam) of the sample to be analyzed, and automatically outputs the copy number of the CYP2D6 gene in the sample. The flowchart is shown below. Figure 2 .

[0073] In one embodiment, the star-shaped nomenclature for the CYP2D6 gene copy number type is defined according to the PharmVar database (https: / / www.pharmvar.org / gene / CYP2D6). This database provides standardized nomenclature for the CYP family and gene loci, and is an authoritative database for pharmacogenomics research.

[0074] In one embodiment, this invention does not involve CYP2D6*3, CYP2D6*4, CYP2D6*7, and CYP2D6*10, as these are single SNP variants / indel variants, which are relatively easy to detect using other techniques or NGS. Therefore, this invention does not integrate these variants and focuses on analyzing CYP2D6 copy number variants (*5 whole gene deletion).

[0075] In one embodiment, CYP2D6*5 represents the CYP2D6 full gene deletion, specifically referring to the standard nomenclature database https: / / www.pharmvar.org / gene / CYP2D6.

[0076] In one embodiment, *1 is generally used to indicate that the target variant was not detected.

[0077] Example 1

[0078] The main steps of this embodiment are as follows:

[0079] 11. Depth of Standardization

[0080] In this embodiment, the known CYP2D6 genotype of the samples was obtained by long-PCR followed by agarose gel electrophoresis. This validation method synthesizes information from multiple studies and designs specific primers for specific regions.

[0081] The principle is as follows: Three pairs of primers are designed upstream, downstream and inside the CYP2D6 gene, respectively. The PCR products are approximately 2.5Kb, 3.5Kb and 5.1Kb. The copy number type of the sample is determined based on the band distribution of the PCR products after gel running.

[0082] Reference information:

[0083] 1) Selection and optimization of large-fragment PCR method for CYP2D6*5 genotyping (https: / / www.doc88.com / p-7428798324659.htm), Table 1 of the article mentions traditional large-fragment PCR and Figure 1 2) PMID18957039: {Development of a PCR-based strategy for CYP2D6 genotyping, including gene multiplication of worldwide potential use}. See the section on "Detection of the CYP2D6*5Allele" in that paper for details.

[0084] This embodiment uses conventionally constructed capture sequencing data (also known as panel sequencing data) on an MGISEQ-2000, PE100 sequencing platform, with insert lengths of approximately 250–300 bp. First, the sequencing depth of 800 samples (whole blood genomic DNA) with known CYP2D6 gene copy number types was calculated at each location of the gene (i.e., each capture site). Simultaneously, the average depth of each sample within the capture interval was calculated. The normalized depth of each location on the CYP2D6 gene was obtained by dividing the sequencing depth at each location by the average depth of the samples, as shown in the following formula:

[0085] (normalized_depth) ij ) m×n =(depth) ij / mean_depth i ) m×n , where depth ij The mean_depth represents the depth of the j-th position of the CYP2D6 gene in the i-th sample. i Normalized_depth represents the average depth of the i-th sample. ij The normalized depth of the CYP2D6 gene at the j-th position in the i-th sample is represented by m, where m represents the number of samples (800 in this example) and n represents the length of the CYP2D6 gene (4312 in this example). This step can eliminate the influence of different experimental conditions and sequencing data volume of samples within the batch.

[0086] 12. Baseline Correction

[0087] In this embodiment, the median normalized depth of the CYP2D6 gene at each location in 22 samples with a known CYP2D6 gene copy number of 0N was used as the calibration baseline: (baseline)j ) 1×n =(median(normalized_depth) ·j )) 1×n baseline j Let `j` be the baseline depth of the CYP2D6 gene at position j. Then, subtract the baseline from the normalized depth of each position of this gene in the sample to obtain the corrected depth: `(fixed_depth)`. ij ) m×n =(normalized_depth ij -baseline j ) m×n This is used for subsequent gene copy number analysis, where fixed_depth ij This represents the baseline correction depth of the CYP2D6 gene at the j-th position in the i-th sample. This step can eliminate the influence of factors such as CYP2D7 and CYP2D8 homologous genes and different sequencing chip platforms between batches.

[0088] 13. Feature Extraction

[0089] This embodiment, based on the alignment results of the actual sequencing reads of samples with a CYP2D6 gene copy number type of 0N, found that exons 1, 3, 5, 6 and introns 2, 5, 6 (reference transcript NM_000106.6) can clearly reflect the CYP2D6 copy number type of the sample (e.g., Figure 3 As shown in the figure, these regions are specifically aligned with the sequencing read lengths of homologous genes for CYP2D6. Furthermore, the alignment results of the actual sequencing read lengths in these regions for 1N and 2N samples also reflect the copy number genotype of CYP2D6 in the samples well. Therefore, this embodiment uses the cluster centers at the baseline-corrected depth of the above-mentioned regions for copy number analysis.

[0090] For each region of CYP2D6, firstly, the baseline-corrected depth of a location within the region is randomly selected as the initial cluster center. Then, the sum of squares of the distances between the baseline-corrected depths of other locations within the region and the cluster center is calculated. The cluster center is moved and the second step above is repeated until the sum of squares of the distances reaches its minimum value. The cluster center at this point is used as the feature value of the region.

[0091] 14. Model Training

[0092] This embodiment uses the Random Forest ensemble learning algorithm to train a model on the feature values ​​of the above seven regions from 800 samples with known CYP2D6 gene copy numbers. The core idea of ​​the Random Forest algorithm is to select a strong classifier by voting on the classification results of several weak classifiers.

[0093] This embodiment uses CART decision trees as weak classifiers, employing the Gini coefficient to select the optimal feature and determine its optimal binary split point. Assuming the training set size is M, for each decision tree, M training samples are first randomly and with replacement drawn from the training set. Each sample has K features; k (k <= K) features are randomly selected from all features, and the optimal feature is chosen as the node to build the CART decision tree. The value of k remains constant during the tree's growth. Each decision tree is generated to the maximum extent possible without pruning. The above steps are repeated until the number of trained decision trees reaches a preset value. Finally, each decision tree votes on the classification of the input sample's features; the category with the most votes is the classification result of the strong classifier, thus achieving accurate prediction of the CYP2D6 gene copy number.

[0094] The model trained in this embodiment achieves 100% prediction accuracy on 800 samples with known copy numbers. To further demonstrate the reliability of the parameters used in the model, 10-fold cross-validation was performed. First, all samples were randomly divided into 10 equal-sized subsets. Then, these 10 subsets were iterated sequentially, with the current subset serving as the validation set and the remaining 9 subsets as the training set for model training and prediction accuracy evaluation. The model used in this embodiment achieved 100% prediction accuracy in 10-fold cross-validation.

[0095] The training and copy number prediction methods for the CYP2D6 gene copy number prediction model in this embodiment are as follows:

[0096] 21. Obtaining the sample BAM file

[0097] First, the raw FASTQ data from high-throughput sequencing is preprocessed to remove adapter sequences, low-quality sequences, and low-quality portions of the sequences. Software used includes any one or a combination of at least two of SOAPnuke, FastP, or Trimmomatic (SOAPnuke is used in this embodiment). Then, the preprocessed clean sequence data is aligned to the hg19 or hg38 human reference genome (hg19 is used in this embodiment) using BWA software. The aligned results are then sorted using any one of samtools, Picard, or GATK (Picard is used in this embodiment). Finally, GATK software is used to mark duplicates in the sorted results and perform base quality correction to obtain the input BAM file required in this embodiment.

[0098] 22. Feature Extraction

[0099] A large number of samples with known CYP2D6 gene copy numbers were selected as the training set (800 samples were used in this embodiment, and the number of samples with different copy numbers was kept as equal as possible). The pysam module was used to calculate the average depth and the sequencing depth of each position of the CYP2D6 gene for each sample, and then standardized them. Then, the median of the standardized depth of the CYP2D6 gene at each position in the ON samples was calculated as the calibration baseline, and baseline correction was performed on the standardized depth of all samples. Finally, the cluster centers of the standardized baseline-corrected depths of exons 1, 3, 5, and 6 and introns 2, 5, and 6 of the CYP2D6 gene were calculated as feature values.

[0100] 23. Training of the CYP2D6 gene copy number prediction model

[0101] Based on the description above, feature values ​​of exons 1, 3, 5, and 6, and introns 2, 5, and 6 of the CYP2D6 gene for each sample were extracted. Then, the `sklearn.ensemble.RandomForestClassifier` function (with default parameters) from the scikit-learn module was used to train the model on all samples using these seven feature values. In this embodiment, the model's prediction accuracy on the training set reached 100%. Simultaneously, 10-fold cross-validation was used to demonstrate the reliability of the model parameters. An accuracy of over 99% across 10 iterations indicates that the model parameters are optimal; otherwise, increasing the sample size or adjusting the model parameters (increasing the number of base decision trees) should be considered.

[0102] 24. Validate the model's performance using the test set.

[0103] The trained model was used to perform feature extraction and copy number prediction on 50 samples. The feature value information, detection results, and consistency judgment results are as follows:

[0104] Table 1

[0105]

[0106]

[0107] The results in the table above show that all results are consistent with the expected results, with an accuracy rate of 100%.

[0108] In one embodiment, based on high-throughput sequencing data and a random forest ensemble learning algorithm, the cluster centers of the normalized baseline-corrected depth of exons 1, 3, 5, and 6 and introns 2, 5, and 6 of the CYP2D6 gene are used as feature values. This invention enables high-throughput detection of the CYP2D6 gene copy number. Traditional PCR electrophoresis cannot achieve high-throughput detection and is cumbersome. Comparative genomic hybridization and mass spectrometry require control samples for genotyping and are greatly affected by sample quality, placing high demands on the experimental environment, technical platform, and personnel. This invention can achieve high-throughput genotyping while avoiding complex experimental procedures.

[0109] In one embodiment, the present invention also effectively avoids the limitations of CYP2D6 copy number typing software such as Stargazer, Aldy, and Cyrius, which are only applicable to specific panels / WGS or specified reference genomes, and does not require reference set correction, thus avoiding the impact of sequencing fluctuations on the analysis results.

[0110] In one embodiment, the analysis method of the present invention can be applied to high-throughput sequencing data such as whole-genome sequencing and targeted capture sequencing.

[0111] With the development of high-throughput sequencing technology and the decrease in cost, each individual can obtain a large amount of sequencing data. However, most testing methods only perform targeted analysis for the purpose of the test, resulting in a significant waste of data. In one embodiment, the analysis method of the present invention can reanalyze existing high-throughput sequencing data to obtain CYP2D6 gene copy number genotyping results. Combined with conventional SNP locus genotyping, the CYP2D6 enzyme metabolotype can be comprehensively determined to guide personalized medication, providing secondary benefits to the tested individuals. Furthermore, because this analysis method is easy to operate and highly efficient, it can achieve large-scale population-level CYP2D6 copy number analysis, promoting in-depth research in pharmacogenomics and guiding personalized clinical medication.

[0112] In one embodiment, the present invention creatively uses the cluster centers of the normalized baseline-corrected depth of the regions of exons 1, 3, 5, 6 and introns 2, 5, 6 of the CYP2D6 gene as feature values, and combines them with the random forest ensemble learning algorithm to achieve fast and accurate prediction of the CYP2D6 copy number.

[0113] In one embodiment, the present invention detects the copy number of the cytochrome P450 isoenzyme 2D6 gene (CYP2D6), and its application scope is all genomic products or services involving CYP2D6 genotyping. The detection subjects are healthy people or people who need personalized medication guidance, and the market application prospects are very broad.

[0114] Those skilled in the art will understand that all or part of the functions of the various methods in the above embodiments can be implemented by hardware or by computer programs. When all or part of the functions in the above embodiments are implemented by computer programs, the program can be stored in a computer-readable storage medium, which may include: read-only memory, random access memory, disk, optical disk, hard disk, etc., and the program is executed by a computer to achieve the above functions. For example, the program can be stored in the memory of a device, and when the program in the memory is executed by the processor, all or part of the above functions can be achieved. In addition, when all or part of the functions in the above embodiments are implemented by computer programs, the program can also be stored in a server, another computer, disk, optical disk, flash drive, or external hard drive, etc., and can be downloaded or copied to the memory of a local device, or the system of the local device can be updated. When the program in the memory is executed by the processor, all or part of the functions in the above embodiments can be achieved.

[0115] The above examples illustrate the present invention only to aid in understanding it and are not intended to limit the scope of the invention. Those skilled in the art can make various simple deductions, modifications, or substitutions based on the principles of this invention.

Claims

1. A method for predicting the copy number type of a target gene, characterized in that, include: The baseline correction step includes using the median and / or average sequencing depth of each capture site of the target gene in the sequencing data of samples with known target gene copy number typology as the target gene whole-genome deletion genotype as the calibration baseline; calculating the baseline depth of each capture site of the target gene; and subtracting the baseline depth from the standardized depth of each capture site of the target gene in the sequencing data of the sample to be tested to obtain the corrected depth; the target gene includes the CYP2D6 gene; the target gene whole-genome deletion genotype refers to both alleles of the target gene being whole-genome deletion genotype, i.e., target gene *5 / *5; the copy number typology includes at least one of 0N, 1N, and 2N, where 0N represents the target gene *5 / *5 typology, 1N represents the target gene *1 / *5 typology, and 2N represents the target gene *1 / *1 typology; The feature extraction step includes obtaining the cluster center of the specific region of the target gene based on the corrected depth of each capture site, and using the cluster center as the feature value of the corresponding region. The specific region includes at least one region from exons 1, 3, 5, 6 or introns 2, 5, 6 of the target gene; The feature extraction step includes: a distance sum of squares calculation step. For at least one region among exons 1, 3, 5, 6 and introns 2, 5, 6 of the target gene, firstly, the baseline-corrected depth of a capture site within the region is randomly selected as the initial cluster center. Then, the sum of squares of the distances between the baseline-corrected depths of other capture sites within the region and the cluster center is calculated. This step is repeated, moving the cluster center capture site and repeatedly calculating the sum of squares of the distances according to the distance sum of squares calculation step, until the sum of squares of the distances reaches its minimum value. The cluster center at this point is then used as the feature value of the region. The prediction step includes predicting the copy number type of the target gene in the sequencing data of the sample to be tested based on the feature value.

2. The method as described in claim 1, characterized in that, In the feature extraction step, the specific region includes all regions of exons 1, 3, 5, and 6 and introns 2, 5, and 6 of the target gene.

3. The method as described in claim 1, characterized in that, In the baseline correction step, the sequencing depth of each capture site is a normalized depth.

4. The method as described in claim 1, characterized in that, In the baseline correction step, the baseline depth is calculated using the following formula: ； in, The baseline depth is the j-th capture site of the target gene. Specifically, it refers to the median of the normalized depth of the j-th capture site in all 0N samples. 0N samples are target gene *5 / *5 type samples, that is, target gene samples in which both alleles are whole gene deletions.

5. The method as described in claim 1, characterized in that, In the baseline correction step, the formula for calculating the standardized depth is as follows: ； in, This represents the depth of the j-th capture site of the target gene in the sequencing data of the i-th sample. This represents the average depth of the i-th sample to be tested. denoted by , m represents the normalized depth of the j-th capture site of the target gene in the sequencing data of the i-th sample to be tested, m represents the number of samples, and n represents the length of the target gene.

6. The method as described in claim 1, characterized in that, In the baseline correction step, the corrected depth is calculated using the following formula: ； in, This represents the baseline-corrected depth of the j-th capture site of the target gene in the sequencing data of the i-th sample. This represents the normalized depth of the j-th capture site of the target gene in the sequencing data of the i-th sample. This represents the baseline depth of the j-th capture site.

7. The method as described in claim 1, characterized in that, The sequencing data refers to either captured sequencing data or whole-genome sequencing data that captures the target gene.

8. The method as described in claim 1, characterized in that, The sequencing data is sequencing data aligned to a reference genome sequence.

9. The method as described in claim 8, characterized in that, The reference genome sequence was derived from a reference group.

10. The method as described in claim 9, characterized in that, The reference genome sequence contains a common sequence from the reference group.

11. The method as described in claim 10, characterized in that, The reference genome sequence contains at least a portion of the hg19 human genome, hg38 genome, hg18 genome, hg17 genome, or hg16 genome.

12. The method as described in claim 1, characterized in that, The sequencing data is generated based on next-generation sequencing.

13. The method as described in claim 1, characterized in that, The prediction step includes using a model to predict the copy number type of the target gene in the sequencing data of the sample to be tested, based on the feature values.

14. The method as described in claim 1, characterized in that, The algorithm of the model includes ensemble learning algorithms.

15. The method as described in claim 14, characterized in that, The ensemble learning algorithm includes at least one of Bagging-based algorithms and Boosting-based algorithms.

16. The method as described in claim 15, characterized in that, The Bagging-based algorithm includes Random Forest.

17. The method as described in claim 16, characterized in that, The Boosting-based algorithm includes at least one of Adaboost, GBDT, and XGBOOST.

18. A device for predicting the copy number type of a target gene, characterized in that, include: The baseline correction module includes using the median and / or average sequencing depth of each capture site of the target gene in the sequencing data of samples with known target gene copy number typology as the target gene whole-genome deletion genotype as the calibration baseline. It calculates the baseline depth of each capture site of the target gene, subtracts the baseline depth from the standardized depth of each capture site of the target gene in the sequencing data of the sample to be tested, and obtains the corrected depth. The target gene includes the CYP2D6 gene. The target gene whole-genome deletion genotype refers to both alleles of the target gene being whole-genome deletion, i.e., target gene *5 / *5. The copy number typology includes at least one of 0N, 1N, and 2N, where 0N represents the target gene *5 / *5 typology, 1N represents the target gene *1 / *5 typology, and 2N represents the target gene *1 / *1 typology. The feature extraction module includes obtaining the cluster center of the specific region of the target gene based on the corrected depth of each capture site, and using the cluster center as the feature value of the corresponding region. The specific region includes at least one region from exons 1, 3, 5, 6 or introns 2, 5, 6 of the target gene; The feature extraction step includes: a distance sum of squares calculation step. For at least one region among exons 1, 3, 5, 6 and introns 2, 5, 6 of the target gene, firstly, the baseline-corrected depth of a capture site within the region is randomly selected as the initial cluster center. Then, the sum of squares of the distances between the baseline-corrected depths of other capture sites within the region and the cluster center is calculated. This step is repeated, moving the cluster center capture site and repeatedly calculating the sum of squares of the distances according to the distance sum of squares calculation step, until the sum of squares of the distances reaches its minimum value. The cluster center at this point is then used as the feature value of the region. The prediction module includes predicting the copy number type of the target gene in the sequencing data of the sample to be tested based on the feature value.

19. An apparatus, characterized in that, include: Memory, used to store programs; A processor for implementing the method as described in any one of claims 1 to 17 by executing a program stored in the memory.

20. A computer-readable storage medium storing a program that can be executed by a processor to implement the method as claimed in any one of claims 1 to 17.

Images