Method and device for detecting copy number variation

Abstract

The application provides a copy number variation detection method and device, and belongs to the technical field of bioinformatics. The method comprises the following steps: calculating the logR value of each bin; correcting the logR value of each bin to obtain a corrected Z-score; counting the corrected Z-score value according to the chromosome arm; determining the CNV state of each chromosome arm according to the comparison between the corrected Z-score value and a threshold in each chromosome arm; calculating the noise-reduced Z-score of each bin; screening significant windows; and merging adjacent significant windows into candidate peaks. The logR value of each bin is clustered to obtain five components. The posterior probability of each candidate peak is calculated based on the five components, and the local CNV state of the candidate peak is determined according to the posterior probability. The application realizes accurate identification of the key CNV state of a tumor through robust baseline estimation, arm-level decision-making of spatial autocorrelation correction and the like.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics technology, specifically relating to a method and apparatus for detecting copy number variations. Background Technology

[0002] Copy number variations (CNVs) are of paramount importance in tumor detection, forming a core component of tumor genomic heterogeneity. In recent years, molecular pathology has revealed a series of copy number variations (CNVs) closely related to tumor development and progression, such as +7 / -10, 1p / 19q co-deletion, EGFR amplification, CDKN2A / B homozygous deletion, PDGFRA amplification, MDM2 amplification, and PTEN deletion. These CNV events are not only key driver variants but also directly influence targeted therapy selection (such as EGFR inhibitors, CDK4 / 6 inhibitors, PDGFRA inhibitors, and PRMT5 inhibitors) and prognostic assessment.

[0003] Currently, the main methods for detecting CNVs in clinical practice include fluorescence in situ hybridization (FISH), microarray comparative genomic hybridization (aCGH), and bioinformatics analysis based on next-generation sequencing (NGS). NGS-based CNV detection, due to its advantages such as high throughput, adjustable resolution, and the ability to simultaneously detect other variant types, has gradually become the mainstream method for scientific research and clinical auxiliary diagnosis. However, existing CNV analysis algorithms for tumor samples still have the following significant shortcomings in practical applications:

[0004] The distinction between local events and chromosome arm-level events is blurred: whether amplification originates from focal high-copy amplification (such as extrachromosomal DNA, ecDNA) or from chromosome arm-level amplification has entirely different clinical significance. The former often indicates high invasiveness and drug resistance, while the latter may be associated with chromosomal instability. However, many algorithms only output a single level of CNV status or use a fixed-width threshold for rigid classification, leading to misclassification of event types. Summary of the Invention

[0005] To address the aforementioned issues, this invention provides a method and apparatus for detecting copy number variations (CNVs). Through robust baseline estimation, arm-level decision-making with spatial autocorrelation correction, local event detection with adaptive noise correction, and empirical Bayesian classification, it achieves accurate identification of key CNV events in tumors, thereby providing more reliable technical support for precision clinical diagnosis and treatment.

[0006] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0007] On one hand, the present invention provides a method for detecting copy number variations, comprising: acquiring sequencing data of a sample to be tested, and calculating the logR value of each bin based on the sequencing data; fitting the logR distribution using a Gaussian mixture model to estimate the normal diploid baseline logR; correcting the logR value of each bin based on the normal diploid baseline logR to obtain a corrected Z-score; statistically analyzing the corrected Z-score values ​​by chromosome arm, and determining the CNV status of each chromosome arm based on the comparison between the corrected Z-score value and a threshold within each chromosome arm; correcting the logR of each bin based on the noise of the chromosome arm and the whole genome to obtain a noise-reduced Z-score for each bin; screening significant windows and merging adjacent significant windows into candidate peaks; clustering the logR values ​​of each bin to obtain five components; calculating the posterior probability of each candidate peak belonging to each group based on the five components, and determining the local CNV status of the candidate peak based on the posterior probability; and outputting the CNV status of the chromosome arm and / or the local CNV status.

[0008] Furthermore, the normal diploid baseline is obtained as follows: a Gaussian mixture model is used to fit the logR distribution, and the number of components G∈{1,2,3,4,5} is set. The optimal model is selected by Bayesian information criterion. The mean and standard deviation of logR of each component are extracted from the fitting results, and the component with the mean closest to 0 is selected as the normal diploid baseline logR. If the Gaussian mixture model fails or the absolute value of the mean is >0.5, the model is regressed to robust estimation, and the normal diploid baseline is calculated based on the median of logR of each component.

[0009] Furthermore, the robust estimate is calculated as follows:

[0010] ;

[0011] ;

[0012] in, The mean of the normal diploid baseline logR; The standard deviation of the normal diploid baseline logR; Let R be the logR of the i-th bin in component G. Component G has n bins. Calculate the mean and standard deviation of logR for each of the five components.

[0013] Furthermore, the logR values ​​of each bin are corrected based on the normal diploid baseline logR to obtain the corrected Z-score values:

[0014] ;

[0015] Zi This represents the Z-score value of the i-th bin. Let R be the logR of the i-th bin. The mean of the normal diploid baseline logR; Let logR be the standard deviation of the normal diploid baseline; then, the corrected Z-score values ​​for all bins are fitted using a Gaussian mixture model to obtain the mean Z-score of the normal diploid baseline. and standard deviation .

[0016] Furthermore, the CNV status of each chromosome arm is determined by comparing the corrected Z-score value within each chromosome arm with a threshold, including:

[0017] ;

[0018] ;

[0019] Where T C This represents the number of bins contained in chromosome arm C. Z-score of chromosome arm C The number of bins, Z-score of chromosome arm C The number of bins, The mean of the Z-score is... The standard deviation of the Z-score;

[0020] The method for determining the CNV state S of a chromosome arm is as follows:

[0021] ;

[0022] In this context, Gain represents amplification, Loss represents deletion, and None represents normal.

[0023] Furthermore, the Z-score after denoising for each bin is calculated, and significant windows are selected, including:

[0024]

[0025] in:

[0026] ;

[0027] ;

[0028] ;

[0029] Suppose that chromosome arm C contains a total of T bins. Let logR be the value of the i-th bin in chromosome arm C. This represents the median absolute deviation of logR in chromosome arm C; Let be the variance of logR for all bins in chromosome arm i; The standard deviation among all bins after noise reduction of chromosome arm C; The corrected Z-score was recalculated using chromosome arms and genome-wide background noise.

[0030] ;

[0031] in, Let R be the Z-score value of the i-th bin of chromosome arm C after noise reduction, logR i μ is the coverage value of the i-th bin. b The mean of logR at the baseline of normal diploid individuals; calculated using a standard normal distribution. The two-sided p-values ​​were obtained, and FDR correction was performed on the p-values ​​of all bins to obtain q. i If q i If the value is less than 0.05, the bin is marked as significant; other chromosome arms are calculated using the algorithm described above.

[0032] Furthermore, adjacent significant windows are merged into candidate peaks, including: merging significant bins according to chromosome, arm, and positional continuity. If adjacent significant bins have the same chromosome, the same arm, and the difference between the start position and the end position of the previous window is ≤1bp, they are classified into the same candidate peak.

[0033] Furthermore, regarding the corrected Clustering was performed using k-means clustering, which automatically identified five components: normal N, high-copy amplification H, medium-copy amplification G, heterozygous deletion L, and homozygous deletion M.

[0034] Furthermore, the calculation of the posterior probability includes: using the proportion of each component after Laplace smoothing as the prior probability; calculating the likelihood of the candidate peak belonging to each component based on the logR mean of the candidate peak; multiplying the likelihood by the prior probability and normalizing it to obtain the posterior probability of the candidate peak belonging to each component; and determining whether the candidate peak is a high-copy amplification, medium-copy amplification, heterozygous deletion, or homozygous deletion based on the comparison result of the posterior probability with the preset threshold θ, and the relationship between the logR mean of the candidate peak and the mean and standard deviation of each component.

[0035] This invention also provides a device for detecting copy number variations, comprising: a calculation module for acquiring sequencing data of a sample to be tested and calculating the logR value of each bin based on the sequencing data; a CNV status discrimination module at the chromosome arm level for fitting the logR distribution using a Gaussian mixture model to estimate the normal diploid baseline logR; correcting the logR value of each bin based on the normal diploid baseline logR to obtain a corrected Z-score; statistically analyzing the corrected Z-score values ​​by chromosome arm; and determining the Z-score value within each chromosome arm relative to a threshold. The comparison determines the CNV status of each chromosome arm; the local CNV status discrimination module corrects the logR of each bin based on the noise of the chromosome arm and the whole genome to obtain the noise-reduced Z-score of each bin, filters significant windows, and merges adjacent significant windows into candidate peaks; the logR values ​​of each bin are clustered to obtain five components; the posterior probability of each candidate peak belonging to each group is calculated based on the five components, and the local CNV status of the candidate peak is determined based on the posterior probability; the output module outputs the CNV status of the chromosome arm and / or the local CNV status.

[0036] The beneficial effects of the technical solutions provided by the embodiments of the present invention include:

[0037] (1) In chromosome arm level detection, the GMM model is used to evaluate the baseline and perform secondary Z-score correction. The CNV status of each chromosome arm is determined by comparing the corrected Z-score value with the threshold. In local event detection, a dynamic k-means five-category classification (homozygous deletion, heterozygous deletion, normal, local intermediate copy number amplification, and local high copy number amplification) is pioneered. Combined with adaptive Z-score and empirical Bayes, the posterior probability integrates component distribution information, achieving accurate discrimination of the five types of events. At the same time, the fine distinction between local high copy number amplification and intermediate copy number amplification events can provide more refined guidance for treatment and prognosis. For example, in glioma, local high copy number amplification of EGFR suggests that the gene may exist in the form of ecDNA, resulting in stronger treatment resistance and a worse prognosis.

[0038] (2) Baseline estimation is performed using GMM and k-means unsupervised learning methods to automatically identify diploid peaks from tumor data. Adaptive heteroscedasticity correction is used to eliminate noise within and between chromosome arms, increasing detection robustness and enabling CNV detection without the need for paired normal samples, significantly simplifying the clinical analysis process. Validated using low-depth WGS (sWGS) and glioma clinical samples, this method achieves over 95% accuracy in detecting key markers such as EGFR amplification, CDKN2A / B homozygous deletion, 1p / 19q co-deletion, and +7 / -10, maintaining robust performance even with tumor concentrations as low as 0.05. Simultaneously, it can detect CNVs from tissues, FFPE, cerebrospinal fluid, and peripheral blood sources, including genomic DNA (gDNA) and cfDNA, in multiple cancers such as brain tumors, bladder cancer, and lung adenocarcinoma brain metastases. Compared with existing technologies, this invention has outstanding substantive features and significant advancements in statistical modeling and event segmentation, demonstrating promising clinical application prospects. Attached Figure Description

[0039] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0040] Figure 1 This is a distribution diagram showing the accuracy of the detection method provided in Embodiment 2 of the present invention under different tumor contents;

[0041] Figure 2 This is a special case illustrating the inconsistency between the detection method provided in Embodiment 4 of the present invention and the FISH detection of CDKN2A / B;

[0042] Figure 3 This is a CNV diagram of a case in gDNA derived from FFPE in glioma provided in Example 6 of the present invention;

[0043] Figure 4 This is the CNV map of urinary gDNA from bladder cancer provided in Example 7 of the present invention;

[0044] Figure 5 This is a CNV map of cfDNA from cerebrospinal fluid and peripheral blood of lung adenocarcinoma provided in Example 8 of the present invention. Detailed Implementation

[0045] The present invention will be further described in detail below through specific embodiments. However, those skilled in the art will understand that the following embodiments are for illustrative purposes only and should not be considered as limiting the scope of the invention. Where specific techniques or conditions are not specified in the embodiments, they are performed in accordance with the techniques or conditions described in the literature in the field or according to the product instructions. Where the manufacturers of the reagents or instruments used are not specified, they are all conventional products that can be obtained commercially.

[0046] As used herein, the terms “comprising,” “including,” “having,” or any other variations thereof are intended to cover non-exclusive inclusion. For example, a process, method, article, or apparatus that includes listed elements is not necessarily limited to those elements, but may include other elements not expressly listed or inherent to such process, method, article, or apparatus. Unless the context clearly specifies otherwise, the singular forms “an” and “the” include a plural of objects under discussion.

[0047] Unless otherwise defined, all technical and scientific terms used in this invention have the same meaning as commonly understood by those skilled in the art. In addition to the specific methods, apparatus, and materials used in the embodiments, this invention can be implemented using any prior art methods, apparatus, and materials similar to or equivalent to those described in the embodiments of this invention, based on the knowledge of those skilled in the art and the description of this invention.

[0048] Unless otherwise stated, the experimental methods, detection methods, and preparation methods disclosed in this invention all employ conventional techniques in molecular biology, immunology, laboratory medicine, gene sequencing technology, bioinformatics technology, and related fields.

[0049] Example 1

[0050] This invention provides a method for detecting copy number variations, comprising:

[0051] S1. Obtain the sequencing data of the sample to be tested, and calculate the logR of each bin based on the sequencing data;

[0052] S2. Use a Gaussian mixture model to fit the logR distribution and estimate the baseline logR of normal diploid individuals;

[0053] S3. Based on the normal diploid baseline logR, the logR of each bin is corrected to obtain the corrected Z-score. The corrected Z-score values ​​are statistically analyzed by chromosome arm. The CNV status of each chromosome arm is determined by comparing the Z-score value within each chromosome arm with the threshold.

[0054] S4. Based on the noise of chromosome arms and the whole genome, the logR in S1 is corrected to obtain the noise-reduced Z-score of each bin. Significant windows are screened and adjacent significant windows are merged into candidate peaks. The logR of each bin (derived from S1) is clustered to obtain five components. Based on the five components, the posterior probability of each candidate peak belonging to each group is calculated, and the local CNV state of the candidate peak is determined according to the posterior probability.

[0055] S5. Output the CNV status of the chromosome arm and / or the local CNV status.

[0056] Specifically, it includes the following steps:

[0057] 1. After obtaining the sequencing data of the sample to be tested, preprocess it, including: data filtering, alignment and sorting, deduplication and indexing, coverage statistics and coverage correction.

[0058] 1.1 Data Filtering: Trimmomatic (v0.38) software was used to remove adapter sequences, low-quality sequences (Q30 < 80%), and sequences with N (base type undetermined) content greater than 10% from the original data to obtain high-quality clean sequences.

[0059] 1.2 Alignment and Sorting: The obtained clean sequence was aligned to the reference genome (hg19) using the mem module of BWA (0.7.12-r1039) software to obtain the original bam file after alignment; the obtained bam file was sorted using the SortSam command of picard software to obtain the sorted bam file.

[0060] 1.3 Deduplication and Indexing: Use the MarkDuplicates command of picard (https: / / broadinstitute.github.io / picard / ) to deduplicat the sorted BAM files, and then use the index command of samtools (v0.1.20) to build an index on the deduplicated BAM files.

[0061] 1.4 Coverage Statistics: The readCounter function of HMMcopy (Version 1.42.0) software is used to count the number of reads in each fixed-length window (the window size can be set to 50kb-500kb, and this embodiment uses 200kb with no overlap) in the deduplicated BAM file, which is the coverage. Each window is called a bin area.

[0062] 1.5 Coverage Correction: The wig files generated in 1.4 were corrected for GC content and mappability using ichorCNA software. The corrected coverage of each bin was obtained (denoted as logR, where R represents the ratio of the actual number of reads obtained from sequencing in each bin to the theoretical number of reads that fall within that bin), as well as the chromosome number, start position, end position, and gene name of each bin. Based on the centromere position of the reference genome (hg19), each chromosome was divided into a short arm (p arm) and a long arm (q arm), and gene annotation information was added.

[0063] 2. Chromosome arm level CNV status determination, including: estimation of normal baseline using Gaussian mixture model (GMM), genome-wide logR correction, and chromosome CNV status determination.

[0064] 2.1 Baseline estimation.

[0065] The normal diploid baseline logR (the set of logR values ​​corresponding to multiple bins) is obtained as follows: a Gaussian mixture model (GMM) is used to fit the logR distribution, with the number of components G∈{1,2,3,4,5} set, and the remaining parameters set to default. The optimal model is selected using the Bayesian information criterion. The mean (mean of logR values ​​of all bins in each component) and standard deviation (standard deviation of logR values ​​of all bins in each component) of each component are extracted from the fitting results.

[0066] 2.1.1 The component with the mean closest to 0 was selected as the baseline logR for normal diploid individuals, with a mean of μ. b The standard deviation is σ b ;

[0067] 2.1.2 If GMM fails (possibly due to polyploidy) or |μ b If |>0.5, then revert to robust estimation and calculate the normal diploid baseline based on the median logR of each component.

[0068] The robust estimate is calculated as follows:

[0069] ;

[0070] ;

[0071] Where median is the median value; for example... Let R be the logR of the first bin in component G. Component G has n bins. Calculate the mean and standard deviation of logR for each of the five components.

[0072] The introduction of GMM can adaptively identify the main peak, remain robust to aneuploid and highly heterogeneous samples, and avoid global misjudgment due to baseline errors.

[0073] 2.2 Genome-wide logR correction.

[0074] Based on the normal diploid baseline logR, the logR (derived from 1.5) of each bin is corrected to obtain the corrected Z-score for each bin:

[0075] ;

[0076] Z i Represents the Z-score of the i-th bin. Let R be the logR of the i-th bin. The mean of the normal diploid baseline logR described in section 2.1; The standard deviation of the normal diploid baseline logR described in 2.1.

[0077] For all bins, the same Gaussian mixture model fitting method as in step 2.1 is used to obtain the mean Z-score of the normal diploid baseline. and standard deviation .

[0078] Since tumor samples may contain large fragment amplifications / deletions, the actual Z-score distribution often deviates from the standard normal distribution. The quadratic GMM fitting directly estimates the true distribution parameters of the Z-score and obtains a data-adaptive threshold, rather than using 0 and 1.5. This significantly improves the matching degree between the threshold and the current sample noise level, and further enhances the robustness of abnormal signal detection.

[0079] 2.3 Chromosome arm CNV status determination.

[0080] Specifically, the CNV status of each chromosome arm is determined by comparing the corrected Z-score (from 2.2) values ​​of all bins within each chromosome arm with a threshold, including:

[0081] For each chromosome arm C, assuming it contains T bins, define the following statistics:

[0082] ;

[0083] ;

[0084] Where T C This represents the number of bins contained in chromosome arm C. Z-score of chromosome arm C The number of bins, Z-score of chromosome arm C The number of bins;

[0085] The method for determining the CNV state S of a chromosome arm is as follows:

[0086] ;

[0087] Gain represents amplification, Loss represents deletion, and None represents normal.

[0088] 3. Determination of local CNV status in the genome.

[0089] The aim is to identify focal copy number alterations in the genome, including local high-level copy number amplification (Amp), local medium-level copy number amplification (Gain), heterozygous deletion (HetDel), and homozygous deletion (HMD), including: screening significant bins within chromosome arms, k-means clustering and component definition, merging of significant bins and peak filtering, and empirical Bayesian identification of the status of local CNVs.

[0090] 3.1 Screening for significant bins within chromosome arms.

[0091] The noise levels of different chromosome arms vary significantly due to GC bias and alignment quality, and distinguishing between focal amplification and outliers is crucial for detection accuracy; simply using a global Z-score reduces sensitivity. This invention uses intra-arm MAD to resist outliers and incorporates global MSV to smooth the MAD influence of extremely short arms. Specifically, the MAD of each chromosome arm is calculated, the whole-genome background noise MSV is calculated, and the standard deviation is adjusted for each bin using the following formula:

[0092] ;

[0093] ;

[0094] ;

[0095] ;

[0096] Chromosome arm C contains a total of T bins. The logR of the first bin in chromosome arm C. The median absolute deviation of logR in chromosome arm C; for example, Let R be the variance of logR for all bins in chromosome arm 1; This represents the standard deviation among all bins after noise reduction of chromosome arm C. Normally, 23 chromosomes need to be counted, each divided into long and short arms. Seven chromosome arms—13p, 14p, 15p, 21p, 22p, and XY—are excluded from the calculation because they are too short, resulting in a total of 39 chromosome arms.

[0097] Using chromosome arms and genome-wide background noise, the denoised Z-score was obtained:

[0098] ;

[0099] in, Let R be the Z-score value of the i-th bin of chromosome arm C after noise reduction, logR i Let logR be the value of the i-th bin, and μ be the value of the i-th bin. b The mean of logR at the baseline of normal diploid individuals (from 2.1);

[0100] Calculate using the standard normal distribution The bilateral p-values ​​were used, and the p-values ​​of all bins on each chromosome arm were FDR corrected using the Benjamini-Hochberg method, denoted as q. i If q i If the value is less than 0.05, then the bin is marked as significant (significant). i =TRUE). This algorithm can detect weaker events in the low-noise arm, while avoiding false positives due to excessive discretization in the high-noise arm. The adaptive mechanism significantly improves the detection capability of focal events.

[0101] The algorithm described above is used to calculate the other chromosome arms. The algorithm is demonstrated using chromosome C as an example. The other 39 chromosomes are calculated using the same method.

[0102] 3.2 Merging of significant bins and peak filtering.

[0103] Adjacent significant bins are merged into candidate peaks, including:

[0104] Significant bins are merged according to chromosome, arm, and positional continuity. If adjacent significant bins have the same chromosome, the same arm, and the difference between the start position and the end position of the previous window is ≤1bp, they are classified as the same candidate peak.

[0105] Record the mean logR of all bins within each candidate peak interval. And width (number of bins), preferably, candidate peaks with widths between [3, 100] are retained to exclude false positive regions that are too narrow or too wide due to random noise.

[0106] The introduction of this method greatly reduces false positives caused by local outliers resulting from nonspecific amplification introduced by the NGS experimental process, while retaining the detection capabilities for typical ecDNA / local amplification and homozygous deletion.

[0107] 3.3 k-means clustering and component definition.

[0108] For the corrected Clustering was performed using k-means clustering (k=1~5), automatically identifying five components: normal (N), high-copy amplification (H), medium-copy amplification (G), heterozygous deletion (L), and homozygous deletion (M). After sorting the means from largest to smallest, the following was defined:

[0109] Component H: mean μ H Standard deviation σ H The proportion of bins contained in this component to the total number of bins (π) H .

[0110] Component G: mean μ G Standard deviation σ G The proportion of bins contained in this component to the total number of bins (π) G .

[0111] Component N: mean μ N Standard deviation σ N The proportion of bins contained in this component to the total number of bins (π) N .

[0112] Component L: Mean μ L Standard deviation σ L The proportion of bins contained in this component to the total number of bins (π) L .

[0113] Component M: Mean μ M Standard deviation σ M The proportion of bins contained in this component to the total number of bins (π) M .

[0114] Where μ H >μ G >μ N >μ L >μ M .

[0115] 3.4 Empirical Bayesian identification of the state of local CNVs in the genome.

[0116] The proportions of each component obtained in 3.3 are smoothed by Laplace plus 1 to obtain the prior probability. The Gaussian likelihood of each candidate peak belonging to each component is calculated based on the median logR of each candidate peak selected in 3.2. The likelihood is multiplied by the prior probability and normalized to obtain the posterior probability of each candidate peak belonging to each component.

[0117] Taking component H as an example:

[0118] ;

[0119] Where ∅ is the normal distribution density function.

[0120] Laplace smoothing with 1 is used as the prior probability:

[0121] ;

[0122] Then calculate the posterior probability:

[0123] ;

[0124] Similarly, calculate and .

[0125] Based on the comparison results between the posterior probability and the preset threshold θ, and the relationship between the mean logR of the candidate peak and the mean and standard deviation of each component, the candidate peak is determined to belong to which category among high copy amplification, medium copy amplification, heterozygous deletion, or homozygous deletion according to the following criteria.

[0126] like >θ and > If so, it is determined to be local high copy amplification (Amp);

[0127] like >θ and > and If so, it is determined to be local copy gain (Gain);

[0128] like >θ and If the result is negative, it is determined to be a local copy number homozygous deletion (HMD).

[0129] like >θ and > and If the result is negative, it is determined to be a local copy number heterozygous missing (HetDel).

[0130] In this embodiment, the natural selection of the loss minimization criterion (range 0-1) in Bayesian decision theory is used, taking into account sensitivity and specificity, and multi-indicator joint judgment logic (such as the mean of candidate peaks). (The relationship between the mean and standard deviation of each component) was investigated, and a posterior probability threshold of θ=0.5 was set.

[0131] This local detection module, based on low-depth whole-genome sequencing technology and utilizing statistical modeling and lightweight Bayesian methods, achieves reliable identification of focal copy number alterations in tumors without relying on external marker data. Furthermore, it classifies focal copy number variations into two categories, distinguishing between true high-copy-number variations (such as ecDNA) and copy number variations caused by chromosomal amplification, making it more suitable for large-scale cohort studies, clinical auxiliary analysis, and treatment guidance.

[0132] This invention also provides a device for detecting copy number variations, comprising: a calculation module for acquiring sequencing data of a sample to be tested and calculating the logR value of each bin based on the sequencing data; a CNV status discrimination module at the chromosome arm level for fitting the logR distribution using a Gaussian mixture model to estimate the normal diploid baseline logR; correcting the logR value of each bin based on the normal diploid baseline logR to obtain a corrected Z-score; statistically analyzing the corrected Z-score values ​​by chromosome arm; and determining the Z-score value within each chromosome arm relative to a threshold. The comparison determines the CNV status of each chromosome arm; the local CNV status discrimination module corrects the logR of each bin based on the noise of the chromosome arm and the whole genome to obtain the noise-reduced Z-score of each bin, filters significant windows, and merges adjacent significant windows into candidate peaks; the logR values ​​of each bin are clustered to obtain five components; the posterior probability of each candidate peak belonging to each group is calculated based on the five components, and the local CNV event type of the candidate peak is determined according to the posterior probability; the output module outputs the CNV status of the chromosome arm and / or the local CNV status.

[0133] This invention also provides a computer-readable medium containing a program capable of running the above-described copy number variation detection method.

[0134] Obviously, those skilled in the art will understand that the modules or steps of the present invention described above can be implemented using general-purpose computing devices. They can be centralized on a single computing device or distributed across a network of multiple computing devices. Optionally, they can be implemented using computer-executable program code, thereby storing them in a storage device for execution by a computing device, or fabricating them separately as individual integrated circuit modules, or fabricating multiple modules or steps as a single integrated circuit module. Thus, the present invention is not limited to any particular hardware and software combination.

[0135] Example 2: Comparison of chromosome arm CNV status with other algorithms.

[0136] sWGS sequencing was performed on 174 clinical samples of gliomas. The algorithm described in this invention was used to analyze CNV variations across all chromosome arms. A comparative method was used, referring to the algorithm published in invention CN121896357A. Compared to the comparative method, this method demonstrated extremely high accuracy (overall mean 98.8%), specificity (98.9%), and negative predictive value (99.6%) in CNV detection at the whole chromosome arm level. Sensitivity (95.5%) and positive predictive value (92.1%) were also high, and the detection reliability for deletion events was higher than that for amplification events. A very small number of arms (such as 22q amplification and 13q amplification) had wider confidence intervals due to sample size limitations or shorter chromosome arms, but they still have reference value. Specific comparison results are shown in Table 1.

[0137] Table 1. Comparison of the performance of CNV status and comparison methods in gliomas.

[0138] Tumor fraction (TF) is crucial for CNV detection. Therefore, we compared the consistency of two algorithms at the chromosome arm level under different tumor fractions. The results show that the CNV detection method described in this invention can obtain high accuracy in arm-level CNV determination under different tumor fractions. When TF ≥ 0.2, the accuracy of all chromosome arms exceeds 0.95, and some arms reach 1.0 (such as 10p and 10q). Even under the extreme condition of TF as low as 0.1, 95% of the arms still have an accuracy of 100%, and the accuracy of the remaining two arms, 9p and 18p, is 85.7%. This method has good stability and clinical applicability, and is particularly suitable for low-depth WGS CNV analysis of solid tumors such as gliomas. See [link to relevant documentation]. Figure 1 .

[0139] Example 3: Comparison of important CNV status in glioma with NGS.

[0140] To evaluate the consistency between the CNV detection method described in this invention (this method) and the published comparative method (reference method) in glioma core biomarkers, using the comparative method as the reference standard, a comparative analysis was performed on 174 glioma samples for 1p / 19q co-deletion, +7 / -10, CDKN2A / B homozygous deletion, and EGFR amplification status. Different gold standards were used for different detection indicators. 1p / 19q and +7 / -10 were statistically analyzed using the algorithm of patent CN121896357A, CDKN2A / B homozygous deletion was statistically analyzed using the algorithm of patent CN117524301B, and EGFR local high-copy amplification was statistically analyzed using the algorithm of patent CN121896357A. The performance indicators of each indicator are summarized in Table 2. The results show that the CNV detection method of this invention is highly consistent with the comparative method in determining the four core molecular markers of glioma (1p / 19q co-deletion, +7 / -10, CDKN2A / B homozygous deletion, and EGFR amplification). All performance indicators (accuracy, sensitivity, specificity, PPV, and NPV) reach or exceed 0.94, with most exceeding 0.98. This demonstrates that the method of this invention can stably replace the comparative method for detecting key CNV events in gliomas, and has the advantages of not requiring paired normal samples and being compatible with low-depth sequencing data, thus possessing good clinical application value.

[0141] Table 2. Performance comparison results of important CNV states and comparison methods in gliomas.

[0142] Example 4: Comparison of important CNV status in glioma with FISH.

[0143] To verify the clinical accuracy of the CNV detection method described in this invention in glioma core biomarkers, fluorescence in situ hybridization (FISH) was used as the gold standard to compare and analyze the 1p / 19q co-deletion (43 cases), CDKN2A / B homozygous deletion (53 cases), and EGFR amplification (33 cases) statuses of glioma samples. The performance of each indicator is summarized in Table 3. The results show that the method of this invention has extremely high detection accuracy (98%) and sensitivity of 100% for the 1p / 19q co-deletion status. Only one sample that was not identified as co-deleted by FISH was identified as co-deleted by the method of this invention (specificity 93%), which is highly consistent with the gold standard overall. The accuracy for detecting CDKN2A / B homozygous deletion was 91%, the sensitivity was 90%, and the specificity was 91%. The accuracy for detecting EGFR amplification was 91%, the specificity was as high as 95%, and the PPV was 92%, indicating that the positive results are highly reliable.

[0144] Table 3. Comparison of important CNV states and FISH performance in gliomas.

[0145] FISH has limitations in gene detection, being susceptible to factors such as probe design, hybridization efficiency, and the subjectivity of signal interpretation, thus carrying the risk of false positives and false negatives. For example, in this invention, there is an instance where the two detection methods are inconsistent; see [link to relevant documentation]. Figure 2 Clinical FISH testing showed that CDKN2A / B was a non-homozygous deletion, but the algorithm of this invention detected it as a homozygous deletion. The read coverage of the FISH probe region and gene region in sWGS showed that the read coverage of the left probe was higher and the read coverage of the right probe was lower. The final FISH result determined it to be a non-homozygous deletion, but the actual read coverage of the CDKN2A / B gene region was extremely low, indicating a homozygous deletion. This inconsistency is caused by the inherent limitations of FISH, indicating that this invention has more advantages.

[0146] The comparative verification of the three core molecular markers of gliomas mentioned above shows that the CNV detection method of the present invention is highly consistent with the FISH gold standard, and the sensitivity, specificity and predictive value of each parameter are all at a good level. This proves that the method of the present invention can accurately identify key driver CNVs of gliomas in low-depth whole-genome sequencing data without the need for paired normal samples, and has clinical translation potential.

[0147] Example 5: Comparison of local gene CNV detection status with comparative methods.

[0148] sWGS sequencing was performed on clinical glioma samples. The CNV variations of specific chromosome arms and local genes were analyzed using the algorithm in this invention. The results were compared with the copy number results output by the ichorCNA software. Generally, a copy number higher than 8 indicates high copy amplification.

[0149] The performance of each indicator is summarized in Table 4. The local high copy amplification of genes output by the algorithm of this invention is divided into two cases: one is that the chromosome CNV where the gene is located is normal, and the gene is locally high copy amplified (such as samples 13 and 12); the other is that the chromosome CNV where the gene is located is amplified, and the gene is high copy amplified (such as samples 1 and 2). The copy number of the gene in both cases is extremely high.

[0150] In this algorithm, the local copy gain of the gene output is divided into two cases: one is that the chromosome CNV where the gene is located is normal, and the gene has local copy gain (such as samples 9 and 6); the other is that the chromosome CNV where the gene is located is amplified, and the gene has local copy gain (such as samples 4 and 5). In both cases, the copy number of the gene fluctuates around 4.

[0151] All the above results indicate that the invention is reliable and has good stability and clinical applicability.

[0152] Table 4 Comparison of local gene CNV detection status and comparison methods

[0153] Example 6: Identification of CNVs in gDNA derived from FFPE in glioma.

[0154] One FFPE brain tumor sample (Case 1) was identified as having +7 / -10 and CDKN2A / B homozygous deletions using gold standard assays. Another FFPE brain tumor sample (Case 2) was identified as having a 1p / 19q co-deletion using gold standard assays. Here, the gold standard is the result of integrating multiple indicators including NGS, FISH, and clinicopathological studies. Genomic DNA (gDNA) was extracted, and sWGS library construction and sequencing were performed. The algorithm of this invention was used to analyze the above three indicators. The results showed that the detection results of this invention were consistent with those of the gold standard. Specifically, +7 / -10 and CDK2A / B homozygous deletions were detected in Case 1; and a 1p / 19q co-deletion was detected in Case 2. The detection results of this invention are consistent with known results, indicating that this invention can be used for CNV discrimination of gDNA from FFPE sources. See [link to relevant documentation]. Figure 3 .

[0155] Example 7: Identification of CNVs in urinary gDNA from bladder cancer.

[0156] One urine sample from a bladder cancer case (Case 3) was selected. Using a comparative method (invention patent CN121896357A), it was confirmed that the chromosome arms were normal and PTPRB showed high copy number amplification (23 copies). Genomic DNA was extracted, and sWGS library construction and sequencing were performed. The algorithm of this invention was used to analyze the above two indicators. The results showed that the detection results of this invention were consistent with the comparative method. Specifically, in Case 3, normal chromosome arms and high copy number amplification of PTPRB were detected. The detection results of this invention are consistent with known results, indicating that this invention can be used for CNV discrimination of gDNA from bladder cancer urine. See [link to relevant documentation]. Figure 4 .

[0157] Example 8: Identification of CNVs in cerebrospinal fluid and peripheral blood cfDNA from lung adenocarcinoma brain metastases.

[0158] Paired cerebrospinal fluid (Case 4 and follow-up Case 6) and peripheral blood samples (Case 5 and follow-up Case 7) from two cases of lung adenocarcinoma with brain metastases were selected. The CNV status of chromosome arms and genes was confirmed using a comparative method (analysis results from ichorCNA software). cfDNA was extracted, and sWGS library construction and sequencing were performed. The CNV status was analyzed using the algorithm of this invention. The results showed that the detection results of this invention were consistent with the comparative method. Specifically, in Case 4 and its subsequent follow-up cerebrospinal fluid sample (Case 6), amplifications of chromosome arms 1q, 5p, 5q, 16q, and 17q, deletions of chromosome arms 3p, 8p, 10q, 12q, 13q, 15q, 18q, 19q, and 21q, localized high-copy amplifications of EGFR, NRP1, ZEB1, and MYC genes, and homozygous deletions of STK11 and CDKN2A / B genes were detected. No CNV changes were detected in peripheral blood samples from Cases 5 and 7. The detection results of this invention are consistent with known results, indicating that this invention can be used for CNV identification of cfDNA from cerebrospinal fluid and peripheral blood in lung adenocarcinoma brain metastases. See [link to relevant documentation]. Figure 5 .

[0159] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for detecting copy number variation, characterized in that, include: Obtain sequencing data of the sample to be tested, and calculate the logR value of each bin based on the sequencing data; A Gaussian mixture model was used to fit the logR distribution to estimate the baseline logR of normal diploid individuals; Based on the normal diploid baseline logR, the logR value of each bin is corrected to obtain the corrected Z-score. The corrected Z-score values ​​are statistically analyzed by chromosome arm. The CNV status of each chromosome arm is determined by comparing the corrected Z-score value within each chromosome arm with the threshold. Based on the noise of chromosome arms and the whole genome, the logR of each bin is corrected to obtain the noise-reduced Z-score of each bin. Significant windows are screened and adjacent significant windows are merged into candidate peaks. The logR values ​​of each bin are clustered to obtain five components. Based on the five components, the posterior probability of each candidate peak belonging to each group is calculated, and the local CNV state of the candidate peak is determined according to the posterior probability. Output the CNV status of the chromosome arm and / or the local CNV status.

2. The method for detecting copy number variation according to claim 1, characterized in that, The normal diploid baseline is obtained as follows: a Gaussian mixture model is used to fit the logR distribution, the number of components G∈{1,2,3,4,5} is set, and the optimal model is selected by Bayesian information criterion; the mean and standard deviation of logR of each component are extracted from the fitting results, and the component with the mean closest to 0 is selected as the normal diploid baseline. If the Gaussian mixture model fails or the absolute value of the mean is >0.5, then a fallback to robust estimation is performed, and the normal diploid baseline is calculated based on the median logR of each component.

3. The method for detecting copy number variation according to claim 2, characterized in that, The robust estimate is calculated as follows: ； ； in, The mean of the normal diploid baseline logR; The standard deviation of the normal diploid baseline logR; Let R be the logR of the i-th bin in component G. Component G has n bins. Calculate the mean and standard deviation of logR for each of the five components.

4. The method for detecting copy number variation according to claim 1, characterized in that, Based on the normal diploid baseline logR, the logR values ​​of each bin are corrected to obtain the corrected Z-score values: ； Z i This represents the Z-score value of the i-th bin. Let logR be the value of the i-th bin. The mean of the normal diploid baseline logR; The standard deviation of the normal diploid baseline logR; The corrected Z-score values ​​for all bins were then fitted using a Gaussian mixture model to obtain the mean Z-score of the normal diploid baseline. and standard deviation .

5. The method for detecting copy number variation according to claim 4, characterized in that, The CNV status of each chromosome arm is determined by comparing the corrected Z-score value within each chromosome arm with a threshold, including: ； ； Where T C This represents the number of bins contained in chromosome arm C. Z-score of chromosome arm C The number of bins, Z-score of chromosome arm C The number of bins, The mean of the Z-score is... The standard deviation of the Z-score; The method for determining the CNV state S of a chromosome arm is as follows: ； In this context, Gain represents amplification, Loss represents deletion, and None represents normal.

6. The method for detecting copy number variation according to claim 1, characterized in that, Calculate the Z-score after denoising for each bin, and filter for significant windows, including: ； in: ； ； ； Suppose that chromosome arm C contains a total of T bins. Let logR be the value of the i-th bin in chromosome arm C. This represents the median absolute deviation of logR in chromosome arm C; Let be the variance of logR for all bins in chromosome arm i; The standard deviation among all bins after noise reduction of chromosome arm C; The corrected Z-score was recalculated using chromosome arms and genome-wide background noise. ； in, Let R be the Z-score value of the i-th bin of chromosome arm C after noise reduction, logR i μ is the coverage value of the i-th bin. b The mean of logR is the baseline value for normal diploid individuals. Calculate using the standard normal distribution The two-sided p-values ​​were obtained, and FDR correction was performed on the p-values ​​of all bins to obtain q. i ; If q i If the value is less than 0.05, then the bin is marked as significant; The algorithm described above is used to calculate the other chromosome arms.

7. The method for detecting copy number variation according to claim 1, characterized in that, Adjacent salient windows are merged into candidate peaks, including: Significant bins are merged according to chromosome, arm, and positional continuity. If adjacent significant bins have the same chromosome, the same arm, and the difference between the start position and the end position of the previous window is ≤1bp, they are classified as the same candidate peak.

8. The method for detecting copy number variation according to claim 6, characterized in that, For the corrected Clustering was performed using k-means clustering, which automatically identified five components: normal N, high-copy amplification H, medium-copy amplification G, heterozygous deletion L, and homozygous deletion M.

9. The method for detecting copy number variation according to claim 8, characterized in that, The calculation of the posterior probability includes: using the proportion of each component after Laplace smoothing as the prior probability; calculating the likelihood of the candidate peak belonging to each component based on the logR mean of the candidate peak; multiplying the likelihood by the prior probability and normalizing it to obtain the posterior probability of the candidate peak belonging to each component. Based on the comparison between the posterior probability and the preset threshold θ, and the relationship between the mean logR of the candidate peak and the mean and standard deviation of each component, the candidate peak is determined to be high-copy amplification, medium-copy amplification, heterozygous deletion, or homozygous deletion.

10. A device for detecting copy number variation, characterized in that, include: The calculation module acquires the sequencing data of the sample to be tested and calculates the logR value of each bin based on the sequencing data; The CNV status discrimination module at the chromosome arm level uses a Gaussian mixture model to fit the logR distribution and estimate the normal diploid baseline. Based on the normal diploid baseline, the logR value of each bin is corrected to obtain the corrected Z-score. The corrected Z-score values ​​are statistically analyzed by chromosome arm, and the CNV status of each chromosome arm is determined by comparing the Z-score value within each chromosome arm with the threshold. The local CNV state discrimination module corrects the logR of each bin based on the noise of chromosome arms and the whole genome to obtain the noise-reduced Z-score of each bin, filters significant windows, and merges adjacent significant windows into candidate peaks; it clusters the logR values ​​of each bin to obtain five components; it calculates the posterior probability of each candidate peak belonging to each group based on the five components, and determines the local CNV state of the candidate peak based on the posterior probability. The output module outputs the CNV status of the chromosome arm and / or the local CNV status.

Images