A copy number variation detection method, device, equipment and computer readable medium

Abstract

The application discloses a copy number variation detection method, device and equipment and a computer readable medium, and belongs to the technical field of genetic engineering. The method comprises the following steps: dividing a genome into genome bins, generating an information profile of the genome, wherein the information profile comprises read depth signals and alignment qualities of each genome bin; performing global segmentation on the genome according to the information profile, and performing local segmentation on at least part of the genome after the global segmentation to obtain gene fragments and read depth signals and alignment qualities of the gene fragments; taking the read depth signals and the alignment qualities of the gene fragments as classification features, calculating abnormal scores of the gene fragments, and identifying copy number variation regions of the genome. The detection method disclosed in the embodiment of the application can improve the sensitivity of copy number variation detection and is effective and reliable in detecting low-amplitude copy number variations.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of genetic engineering, and particularly relates to a copy number variation detection method, device, equipment and computer readable medium. BACKGROUND

[0002] Copy number variations (CNVs) have a significant impact on human genome diversity and the occurrence of many complex diseases. Detection and identification of copy number variations are of great significance in biology and biomedical fields. Next-generation sequencing (NGS) technology provides rich data for detection of copy number variations, and many copy number variation detection methods based on NGS data have been proposed. However, the sensitivity of these methods is not reliable when detecting low-amplitude copy number variations, especially when the length of the copy number variation is very small. SUMMARY

[0003] In order to solve at least one technical problem existing in the prior art, the embodiments of the present application provide a copy number variation detection method, device, equipment and computer readable medium. The technical solution is as follows:

[0004] In a first aspect, a copy number variation detection method is provided, and the method comprises:

[0005] dividing a genome into genome bins, and generating an information profile of the genome, wherein the information profile comprises read depth signals and alignment qualities of each of the genome bins;

[0006] performing global segmentation on the genome according to the information profile, and performing local segmentation on at least part of the genome after the global segmentation, to obtain gene segments and read depth signals and alignment qualities of the gene segments;

[0007] using the read depth signals and the alignment qualities of the gene segments as classification features, calculating abnormal scores of the gene segments, and identifying a copy number variation region of the genome.

[0008] Further, the dividing the genome into genome bins and generating the information profile of the genome comprises:

[0009] obtaining a test sample and a reference sample of the genome;

[0010] aligning the test sample and the reference sample to obtain an alignment result;

[0011] dividing the genome into the genome bins according to the alignment result;

[0012] calculating read depth signals and alignment qualities in the genome bins to generate the information profile.

[0013] Further, the calculating the read depth signal and alignment quality in the genomic bin, generating the information profile, comprises:

[0014] calculating the original read depth signal in the genomic bin;

[0015] normalizing the original read depth signal for correction.

[0016] Further, the globally segmenting the genome according to the information profile, comprises:

[0017] determining a set of the genomic bins with continuous read depth signals;

[0018] comparing the average value of the read depth signal of the genomic bin with the average value of the read depth signal of the remaining genomic bins according to the maximum statistical quantity;

[0019] if the comparison result meets the variation threshold condition, determining that the continuous genomic bins have the read depth signal corresponding to the variation, and dividing the continuous genomic bins into a gene segment.

[0020] Further, the locally segmenting at least part of the genome after the global segmentation, comprises:

[0021] obtaining a preset segmentation length;

[0022] dividing part of the gene segments into a plurality of continuous and non-overlapping gene fragments according to the segmentation length;

[0023] calculating the read depth signal and alignment quality of the gene fragments.

[0024] Further, after obtaining the gene fragments, the method further comprises:

[0025] de-noising the read depth signal in the gene fragments.

[0026] Further, the taking the read depth signal and alignment quality of the gene fragments as classification features, calculating the abnormal value score of the gene fragments, and identifying the copy number variation region of the genome, comprises:

[0027] expressing the read depth signal and alignment quality of all the gene fragments as a standardization matrix, the read depth signal and alignment quality of one gene fragment in the standardization matrix as a data sample;

[0028] calculating a covariance matrix according to the standardization matrix;

[0029] calculating the eigenvalue and eigenvector of the covariance matrix;

[0030] calculate a projection distance of each of the data samples on the eigenvector as an anomaly score;

[0031] determine an abnormal sample in the data samples according to the anomaly score and a set threshold value;

[0032] determine a baseline according to a read depth signal of the gene fragment corresponding to the abnormal sample, and declare the copy number variation region.

[0033] In a second aspect, a copy number variation detection device is provided, and the device comprises:

[0034] a file generation module configured to divide a genome into genome bins, and generate an information configuration file of the genome, wherein the information configuration file comprises a read depth signal and an alignment quality of each of the genome bins;

[0035] a segmentation module configured to globally segment the genome according to the information configuration file, and locally segment at least part of the globally segmented genome, to obtain a gene fragment and a read depth signal and an alignment quality of the gene fragment;

[0036] a detection module configured to take the read depth signal and the alignment quality of the gene fragment as a classification feature, calculate an anomaly score of the gene fragment, and identify a copy number variation region of the genome.

[0037] In a third aspect, an electronic device is provided, and the device comprises:

[0038] one or more processors; and

[0039] a memory associated with the one or more processors, the memory being configured to store program instructions, which, when executed by the one or more processors, perform the method according to any one of the first aspect.

[0040] In a fourth aspect, a computer readable medium is provided, and the medium stores a computer program, wherein the program, when executed by a processor, implements the method according to any one of the first aspect.

[0041] The technical scheme provided by the embodiments of the present application has the following beneficial effects:

[0042] (1) The detection method disclosed by the embodiments of the present application can improve the sensitivity of copy number variation detection, and is effective and reliable in detecting low-amplitude copy number variations;

[0043] (2) The detection method disclosed by the embodiments of the present application adopts global and local segmentation, which effectively avoids the problem that low-amplitude copy number variations and small-length copy number variations are smoothed;

[0044] (3) The detection method disclosed by the embodiment of the application can make low comparison quality signals present higher abnormal value scores, thereby reducing the interference of mapping errors. BRIEF DESCRIPTION OF DRAWINGS

[0045] In order to more clearly illustrate the technical solutions in the embodiments of the application, the drawings needed to be used in the embodiments will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the application, and other drawings can be obtained by those skilled in the art without creative effort on the basis of these drawings.

[0046] Figure 1 is a flow chart of the copy number variation detection method provided by the embodiment of the application;

[0047] Figure 2 is a structural schematic diagram of the copy number variation detection device provided by the embodiment of the application;

[0048] Figure 3 is an evaluation result diagram of the precision, sensitivity and F1 score of each method in the evaluation experiment;

[0049] Figure 4 is a WLS histogram of each method in the evaluation experiment;

[0050] Figure 5 is a distribution chord diagram of the copy number variation detected by the five methods in the method effectiveness experiment;

[0051] Figure 6 is the number of the copy number variation detected by the five methods in the method effectiveness experiment;

[0052] Figure 7 is a structural schematic diagram of the electronic device provided by the embodiment of the application. DETAILED DESCRIPTION

[0053] In order to make the objects, technical solutions and advantages of the application clearer, the technical solutions in the embodiments of the application will be described clearly and completely with reference to the drawings in the embodiments of the application. Obviously, the described embodiments are only some embodiments of the application, but not all the embodiments. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative effort belong to the protection scope of the application.

[0054] In the past, the detection of copy number variations largely relied on microarray technology. But microarray technology is limited by the number of probes, and it can only detect copy number variations that exist in the reference components designed for the probes. In recent years, next-generation sequencing (NGS) technology has developed rapidly and become the mainstream sequencing method. Subsequently, sequence-based structural variation detection method strategies have also emerged. Among them, the strategy based on read depth is widely used for the detection of copy number variations. The basic idea is that compared with normal regions, regions with copy number amplification will obtain higher read depth, while regions with deletion will have lower read depth. There are many methods based on this strategy, such as CNVnator, FREEC, ReadDepth, GROM-RD and the recently released iCopyDAV, CNV-LOF and CNV_IFTV.

[0055] The first step of the read depth-based method is to align the reads in the genomic coordinates, and then obtain the read depth signal by calculating the average read count in the genomic bins. But the read depth signal is biased in regions with high or low GC content (GC-bias), so it needs to be normalized according to the GC content in the genomic bins, which is widely used in single sample cases. The basic assumption of the read depth-based method is that the read depth signal is proportional to the number of copies in the region. Limited by sequence coverage and disturbed by mapping errors, the signal intensity of low-amplitude copy number variations (gaining no more than 2 copies or losing 1 copy) changes less. At the same time, the lower tumor purity further weakens the signal intensity, leading to low-amplitude copy number variations being easily ignored in the detection process.

[0056] Segmentation is performed after removing GC bias, aiming to cluster adjacent genomic bins with similar read depth signals into the same segment. Segmentation determines the length and location of copy number variations, and popular segmentation algorithms include circular binary segmentation (CBS), mean shift, hidden Markov model (HMM) and logistic regression. For example, CNVnator uses the mean shift algorithm for segmentation. It calculates the mean shift vector in each bin according to the read depth signal in the adjacent genomic bin, and determines the segmentation breakpoint according to the direction of the vector. This method has high sensitivity and positioning accuracy. The segmentation of FREEC is completed by logistic regression. Then, the amplification and deletion of the genome are predicted by selecting the allele content corresponding to the maximum log-likelihood. It can estimate the tumor purity of the sequencing sample, and can estimate the absolute copy number (CN) of the predicted copy number variation. iCopyDAV combines CBS and total variation minimization (TVM) algorithm for segmentation, which makes up for the shortcomings of CBS in segmenting low coverage sequences, enabling it to detect a wider range of copy number variations with high sensitivity and accuracy. However, the above segmentation process is performed on the whole genome (global segmentation), and does not take into account the variability of local read counts. This leads to some local copy number variations with weak signals being ignored, especially for small copy number variations (<6kb). To avoid this problem, CNV-LOF segments from a local perspective. It first divides the target genome into multiple consecutive and non-overlapping regions of the same length, and then segments each sub-region using the CBS algorithm. Finally, an anomaly factor is assigned to each genomic segment to identify copy number variation regions. This method shows high sensitivity to low-amplitude copy number variations and performs well on low-tumor purity data. However, focusing only on local regions may limit its performance on sequencing data with high tumor purity.

[0057] In order to solve the problems in the prior art, the embodiments of the present application provide a copy number variation detection method, device, equipment and computer readable medium, and the specific technical solutions are as follows:

[0058] As shown in Figure 1 A copy number variation detection method, comprising:

[0059] S1, divide the genome into genomic bins, and generate an information configuration file of the genome, the information configuration file comprising: read depth signals and alignment quality of each genomic bin.

[0060] The above-mentioned dividing gene bins can be divided according to a preset division rule, for example, first counting the genome, and dividing the genome bins according to the counting order. The read depth and the alignment quality are calculated in a non-overlapping genome bin of appropriate size, which can reduce the random fluctuations of the read depth caused by noise signals. The read depth signal reflects the read count in the genome bin. The alignment quality reflects the average mapping quality level of the reads contained in the genome bin.

[0061] In one embodiment, step S1 comprises:

[0062] Obtaining a test sample and a reference sample of a genome;

[0063] Aligning the test sample and the reference sample to obtain an alignment result;

[0064] Dividing the genome into genome bins according to the alignment result;

[0065] Calculating the read depth signal and the alignment quality in the genome bins to generate an information profile.

[0066] The above-mentioned alignment of the sequencing sample (Fastq format) and the reference sequence (such as hg38) can generate a comparison file, for example, the alignment can be completed by the BWA-MEM method, and then sorted by the SAMTools software.

[0067] The above-mentioned dividing the genome into genome bins can be divided according to a fixed size, for example, using the b i (i = 1, 2, 3,..., m) to represent the i-th genome bin, and m represents the total number of genome bins. The read depth signal of each genome bin can be calculated by formula (1).

[0068]

[0069] wherein, rd i represents the value of the read depth signal of b i , rc j represents the read count of the j-th position in the genome bin, size_b i represents the size of b i , which is set to 1kb.

[0070] The region with mapping error presents a lower value. In particular, when the read cannot be uniquely mapped to a certain position, the related mapping quality is zero. Therefore, a higher mapping quality value indicates a more reliable alignment. The mapping quality of a genome bin, also called the alignment quality, can be calculated by formula (2):

[0071]

[0072] wherein, mq i represents the value of the read depth signal of b ithe mapping quality value of the jth position in the genomic bin, mapg j represents the mapping quality of the jth position in the genomic bin.

[0073] In one embodiment, the read depth signal and the mapping quality in the genomic bin are calculated, and an information profile is generated, comprising:

[0074] The raw read depth signal in the genomic bin is calculated.

[0075] The raw read depth signal is normalized to correct.

[0076] As mentioned above, GC bias (gas chromatography bias) is one of the main reasons for the inconsistency of read depth signal and sequence coverage. The read depth signal value will be biased in areas with low or high GC content. In order to obtain a representative and accurate read depth signal, a commonly used method is used for correction, and the formula is as follows.

[0077]

[0078] wherein and rd i respectively represent the corrected read depth signal value and the raw read depth signal value of the ith genomic bin b av g represents the average read depth signal value of all genomic bins, r gc b i b

[0079] S2, according to the information profile, globally segmenting the genome, and locally segmenting at least part of the globally segmented genome to obtain gene fragments and read depth signals and mapping qualities of the gene fragments.

[0080] In one embodiment, the global segmentation comprises:

[0081] determining a set of genomic bins with continuous read depth signals;

[0082] comparing the average value of the read depth signal of the genomic bin with the average value of the read depth signal of the remaining genomic bins according to the maximum statistic;

[0083] if the comparison result meets the change threshold condition, it is determined that there is a change in the read depth signal corresponding to the genomic bin in the continuous genomic bin, and the continuous genomic bin is divided into a gene segment.

[0084] As mentioned above, the global segmentation takes the CBS algorithm as an example, and the entire genome is segmented, and b1,…,b m are divided into many segments. In each step, it determines a set of continuous genomic bins b i ,b i+1 ,...,b j(1≤i < j≤m). Then the maximum t-statistic is utilized to compare the average of the read depth signal values from b i to b j with the average of the remaining genomic bins. If the p-value is less than a threshold (typically 0.01), b i and b j (if j < m) can maximize the test statistic and are considered as the location of the change point. In other words, the region from b i to b j is divided into a segment. This process is applied recursively to the entire genome and divides it into multiple segments.

[0085] In one embodiment, the local segmentation comprises:

[0086] obtaining a preset segmentation length;

[0087] dividing the partial gene segment into a plurality of continuous and non-overlapping gene segments according to the segmentation length;

[0088] calculating the read depth signal and alignment quality of the gene segments.

[0089] In the above, after the global segmentation is completed, the generated gene segments are further subjected to local segmentation. In one case, the length of the finally obtained gene segment can be preset, and if the gene segment obtained after the global segmentation is greater than the preset gene segment length, the corresponding gene segment needs to be subjected to local segmentation. This process can effectively identify the copy number variations that are smoothed in large segments, such as low amplitude and small copy number variations. First, the length of the sub-segment (Lrs) is specified, which is an integer multiple of the size of the genomic bin (1 kb). Then the segment with a length greater than Lrs is divided into a plurality of continuous and non-overlapping sub-segments. Each sub-segment has the same length (Lrs), and the last one can be smaller than Lrs. The size of Lrs is related to the resolution of the copy number variation. Generally, a smaller Lrs will provide higher detection resolution and sensitivity, but will result in a large number of false positive events. A larger Lrs will provide higher accuracy, but false negatives are difficult to control. Users can set the size of Lrs according to actual needs. In our study, the size of Lrs is set to 3 kb. After the local segmentation is completed, all segments (gene segments generated by local segmentation and segments without local segmentation) are arranged in order and represented by equation (4).

[0090] RS = {rs1, rs2, rs3,..., rs n} (4)

[0091] where rs i represents the i-th segment, and n represents the total number of segments.

[0092] In one embodiment, step S2 further comprises:

[0093] The read-depth signal in the gene segment is denoised.

[0094] After segmentation, the read-depth signal in the segment needs to be smoothed to remove noise. This is because noise data generated in the sorting and segmentation process can cause new errors. The TV algorithm implements the smoothing process, and the read-depth signal containing noise shows a higher total variance. The TV algorithm restores the original signal by reducing the total variance between adjacent segments, while retaining the edge information well. The smoothing formula of the read-depth signal is as follows:

[0095]

[0096] wherein and respectively represent the original read-depth signal value and the smoothed read-depth signal value of the i-th segment; n represents the number of segments; the former term of the formula represents the fitting error of the original RD value and the smoothed read-depth signal value, and the latter term is the L1 norm of the total variance. λ is the penalty parameter of this term, used to adjust the constraint size of the total variance. The larger the value of λ, the stronger the penalty. When it tends to infinity, all read-depth signal values converge to the same value. When λ is 0, the original signal is retained. The user can specify the value of λ.

[0097] S3, taking the read-depth signal and the alignment quality of the gene segment as classification features, calculating the anomaly score of the gene segment, and identifying the copy number variation region of the genome.

[0098] The above, the airport score of each gene segment can be calculated by using the principal component classifier algorithm (PCC), and the read-depth signal and the alignment quality are taken as two features of the classification algorithm.

[0099] In one embodiment, step S3 comprises:

[0100] S31, the read-depth signal and the alignment quality of the gene segment are represented as a standardized matrix, and the read-depth signal and the alignment quality of one gene segment in the standardized matrix are taken as a data sample;

[0101] S32, calculating the covariance matrix according to the standardized matrix;

[0102] S33, calculating the eigenvalue and eigenvector of the covariance matrix;

[0103] S34, calculating the projection distance of each data sample on the eigenvector as the anomaly score;

[0104] S35, determining the abnormal sample in the data sample according to the anomaly score and the set threshold;

[0105] S36. Determine the baseline based on the read depth signal of the gene fragment corresponding to the abnormal sample, and declare the copy number variation region.

[0106] The above examples illustrate this point:

[0107] In step S31, the read depth signal and the alignment quality can be represented as vectors r and m, i.e., r = [r1, r2, ..., rm]. n ] and m = [m1, m2, ..., m n ], where n represents the number of segments, r i and m i Let represent the read depth signal value and alignment quality value of the i-th segment, respectively; they are the average values ​​of the corresponding signals in that segment. These two features can be represented by a matrix N, where each column vector (r) i ,m i ) T It is represented as a sample.

[0108]

[0109] The standardized matrix N is represented by matrix X.

[0110] In steps S32 and S33, the covariance matrix is ​​calculated:

[0111] Find the eigenvalue-eigenvector pairs of C: (λ1,e1) and (λ2,e2), where λ1≥λ2;

[0112] In step S34, calculate x for each data sample. i The projected distance d on e1 is used as the outlier score:

[0113]

[0114] In step S35, the threshold t is set using the OTSU algorithm. When score(x) i When t > 0, it is considered an abnormal sample.

[0115] In step S36, the baseline is determined based on the read depth signal and the copy number variation region is declared.

[0116] The principal component classifier (PCC) is built on top of principal component analysis (PCA), which is an algorithm commonly used for dimensionality reduction of high-dimensional data. The main principle of PCA is to project the original high-dimensional data onto a lower-dimensional space through linear transformation, and make its variance as large as possible, so as to maximize the effective information of the data. PCA has been applied to copy number variation detection problems as a data correction technique, rather than as the main method for identifying copy number variation regions. The main goal of steps S31-S36 is to project the two-dimensional matrix N onto a one-dimensional vector V, and find abnormal samples according to the projection distance.

[0117] In step S31, two features r and m are normalized to the same scale. This is because the comparison quality value is generally larger than the value of the read signal, and when projected onto a low-dimensional space, the comparison quality will obtain a larger weight in the principal component. The following formula is used to standardize the two features:

[0118]

[0119]

[0120] where r' represents the normalized read depth signal, and r sd respectively represent the mean and standard deviation of the read depth signal. The normalization process of the comparison quality is the same as that of the read depth signal, as shown in equation (8). After standardization, the mean of each feature becomes 0 and the standard deviation becomes 1. This ensures that the two features have the same impact on the principal component variable.

[0121] In step S33, the covariance matrix C can be decomposed into orthogonal vectors associated with eigenvalues, called eigenvectors. Eigenvectors reflect different directions of variance changes in sample data, and eigenvalues represent the variance size of data in the corresponding direction. Eigenvectors e1 with high eigenvalues capture most of the variance in the data and serve as principal component vectors.

[0122] In step S34, the read depth signal is the main feature for identifying copy number variations, so only the projection distance of the sample to e1 needs to be calculated. The outlier score is the weighted Euclidean distance between each sample and the eigenvector e1. Samples with larger outliers indicate potential copy number variation or mapping error regions.

[0123] In step S35, a threshold is set to determine abnormal samples. The distance on e1 is quite different for samples covered by different sequences. To accommodate data with different sequence coverage, we use the OTSU algorithm to calculate the threshold. OTSU is a global binary segmentation algorithm, mainly used for the segmentation of grayscale images. The optimal threshold obtained can maximize the separability of the resulting grayscale levels. It dynamically obtains a threshold by traversing all scores in an interval to maximize the variance between the two classes. In this step, we first convert the outlier score to a floating-point number with two decimal places. Then we traverse the outlier scores between the 35th percentile and the 85th percentile, finding the optimal threshold t each time with an increment of 0.01. Samples with scores higher than t (score(xi)≥t) are considered abnormal samples.

[0124] In step S36, the baseline is defined as the average read depth signal value of the remaining samples after removing the abnormal samples (allowing a 15% error). Abnormal samples with read depth signal values higher than this baseline are identified as copy number increases, and those lower than the baseline are considered losses.

[0125] Based on the copy number variation detection method disclosed in the above embodiments of the present application, as shown in the figure, the embodiments of the present application also provide a copy number variation detection device, comprising: Figure 2

[0126] The file generation module 201 is configured to divide the genome into genome bins and generate an information configuration file of the genome, wherein the information configuration file comprises read depth signals and alignment qualities of each genome bin.

[0127] The segmentation module 202 is configured to globally segment the genome according to the information configuration file, and locally segment at least part of the globally segmented genome to obtain gene segments and read depth signals and alignment qualities of the gene segments.

[0128] The detection module 203 is configured to take the read depth signals and alignment qualities of the gene segments as classification features, calculate abnormal scores of the gene segments, and identify copy number variation regions of the genome.

[0129] Further, the file generation module 201 comprises:

[0130] The input module is configured to obtain test samples and reference samples of the genome.

[0131] The alignment module is configured to align the test samples and the reference samples to obtain alignment results.

[0132] The binning module is configured to divide the genome into genome bins according to the alignment results.

[0133] ​The computing module is configured to calculate the read depth signal and the alignment quality in the genomic bin, and generate an information profile.

[0134] Further, the computing module is specifically configured to:

[0135] calculate the original read depth signal in the genomic bin;

[0136] perform normalization processing correction on the original read depth signal.

[0137] Further, the segmenting module 202 comprises:

[0138] The global segmentation module is configured to:

[0139] determine a set of continuous genomic bins of the read depth signal;

[0140] compare the average value of the read depth signal of the genomic bin with the average value of the read depth signal of the remaining genomic bins according to the maximum statistical quantity;

[0141] if the comparison result meets the variation threshold condition, it is determined that the continuous genomic bins have the read depth signal corresponding to the variation, and the continuous genomic bins are divided into one gene segment.

[0142] Further, the segmenting module 202 comprises:

[0143] The local segmentation module is configured to:

[0144] obtain a preset segmentation length;

[0145] divide the partial gene segment into a plurality of continuous and non-overlapping gene segments according to the segmentation length;

[0146] calculate the read depth signal and the alignment quality of the gene segment.

[0147] Further, the segmenting module 202 further comprises:

[0148] The denoising module is configured to perform denoising processing on the read depth signal in the gene segment.

[0149] Further, the detecting module 203 is specifically configured to:

[0150] express the read depth signal and the alignment quality of the gene segment as a standardized matrix, and the read depth signal and the alignment quality of one gene segment in the standardized matrix as a data sample;

[0151] calculate a covariance matrix according to the standardized matrix;

[0152] calculate the eigenvalue and eigenvector of the covariance matrix;

[0153] calculate the projection distance of each data sample on the eigenvector as an anomaly score;

[0154] determining an abnormal sample in the data sample according to the abnormal score and a set threshold value;

[0155] determining a baseline according to the read depth signal of the gene fragment corresponding to the abnormal sample, and declaring a copy number variation region.

[0156] In order to further illustrate the beneficial effects of the technical solutions disclosed in the present application, the performance of the technical solutions for copy number variation detection disclosed in the present application is evaluated as follows:

[0157] First, a comparative experiment is established on simulated data. The basic facts possessed by the simulated data ensure the reliability of the evaluation. The copy number variation detection method (CNV-PCC) disclosed in the present application is compared with five popular methods (CNVnator, FREEC, GROM-RD, CNV-LOF and CNV_IFTV) in terms of precision, sensitivity and F1 score. In order to ensure the fairness of the experiment, the size of the genome bin of certain methods is adjusted so that they can detect small copy number variations. For example, the size of the genome bin of CNVnator is set to the recommended value (250bp for 30x coverage, 130bp for 20x coverage and 90bp for 10x coverage), and the size of the genome bin of FREEC is set to 1kb. The remaining methods use their default parameters. Subsequently, real samples are used to verify the effectiveness of CNV-PCC, which includes three samples from the 1000 Genomes Project.

[0158] The above detection methods are simulated by using sample data. The comprehensive simulation software SinC and the sequence processing tool seqtk are used to generate the simulation data set. All the simulated data are generated based on chromosome 21 in the reference genome hg38. The coverage is set to 10X, 20X and 30X, and the tumor purity is set to 0.4, 0.5 and 0.6, and 30 repeated samples are simulated for each configuration. The length of the copy number variation is limited to 1kb to 6kb, because most methods perform well in detecting large copy number variations. At the same time, high-amplitude copy number variations (homozygous deletion and copy gain >4) are not considered, because they are also easy to detect. A total of 24 copy number variations are generated for each simulation replication, including 15 amplifications and 9 deletions. The copy number of the amplification is 3 and 4, and all the deletions are hemizygous deletions. CNV-PCC is compared with the five methods on the generated simulation data set. If the declared copy number variation covers 50% of the area of the true copy number variation, it is considered as a true positive event. Precision, sensitivity and F1 score are used as indicators in the evaluation, and the results are shown in Figure 3 Figure 3 ​CNV-PCC, CNV_IFTV, FREEC and CNV-LOF achieved better detection results, while CNVnator and GROM-RD performed poorly in most simulated data. CNV-PCC achieved the highest sensitivity in each dataset, which indicates that it detected the most number of copy number variants. It also had the highest F1 score except for the first dataset. CNV-LOF had the second highest sensitivity on 10X data when tumor purity was below 0.6 and had the largest F1 score when tumor purity was 0.4. This reflects its effectiveness on low purity data. The rest in order were CNV_IFTV, FREEC, GROM-RD and CNVnator. The sensitivity and F1 score of CNV_IFTV surpassed CNV-LOF as coverage increased, ranking second. It performed well on high coverage data. The sensitivity and F1 score of CNVnator and GROM-RD were consistently low, even on high coverage data. This indicates that these two methods are not suitable for detection on low purity data. The sensitivity and F1 score of FREEC increased significantly when tumor purity rose to 0.6. However, the sensitivity and F1 score of CNV-LOF decreased on 20X and 30X data. This fact indicates that relying solely on local segmentation can limit its performance on data with high tumor purity. In terms of precision, FREEC had the largest value on 10X data, which produced reliable results. As sequence coverage increased, the precision of CNV-LOF, CNV_IFTV, CNV-PCC surpassed FREEC on 30X data. This is because the strength of read depth signal increases as sequence coverage increases, which reduces the production of false positive results. The precision of CNVnator and GROM-RD was consistently low, as both methods produced a large number of false positive results.

[0159] In summary, CNV-PCC performed well on all simulated data and outperformed the other five methods. CNV_IFTV also showed good performance. CNV-LOF was suitable for CNV detection with low tumor purity, and the results were stable, but the effect decreased when tumor purity was high. In contrast, FREEC was suitable for data with high tumor purity. CNVnator required high coverage in addition to high tumor purity. GROM-RD performed poorly on all simulated data, indicating that it was not suitable for detecting such copy number variants.

[0160] To further evaluate the effectiveness of the six methods, Equation (9) was designed to calculate the weighted length score (WLS) of each method. The WLS histogram of each method is shown in Figure 4

[0161] WLS = W x Lv (9) ​

[0162] where W represents the weight value, which is the ratio of the length of the true copy number variations identified by a method to the length of all copy number variations detected by the method. Lv represents the length of the true copy number variations detected by the method. WLS can indirectly reflect the accuracy level of the identified copy number variations. For example, when the tumor purity is 0.6, the sensitivity of FREEC and CNV_IFTV on 10X data is almost the same, but Figure 4 The WLS of FREEC is higher than that of CNV_IFTV. This is likely because the breakpoint bias of IFTV is greater than that of FREEC. Comparing the WLS between PCC and IFTV on 30X data with a tumor purity of 0.4 can also find the same situation. When the coverage is 30X, the WLS of CNVnator ranks second on data with a tumor purity of 0.6. CNVnator exceeds CNV_IFTV and FREEC, which have higher sensitivity than it. This indicates that CNVnator has a smaller breakpoint bias. CNV-PCC maintains the highest WLS on all simulated data, showing strong performance.

[0163] In terms of real data, the sequencing samples of the trio of the Yoruba family from the 1000 Genomes Project were selected to test the effectiveness of CNV-PCC. It includes NA19238 (mother), NA19239 (father) and NA19240 (daughter). CNV-PCC is applied to the whole genome of each sample, and compared with five existing methods (CNVnator, FREEC, GROM-RD, CNV-LOF, CNV_IFTV). The chord diagram Figure 5 ) shows the distribution of copy number variations detected by the five methods (in kb). The upper half of the circle has 22 sectors, representing chromosomes 1 to 22. The width of each sector represents the number of CNVs detected on that chromosome. The lower half of the circle is divided into five sectors, each representing a method: CNV-LOF (green), CNVnator (red), CNV-PCC (purple), FREEC (black), CNV_IFTV (blue) and GROM-RD (pink). The width of the sector represents the total number of copy number variations detected by the method. It is observed that CNV-LOF detects the most copy number variations, which detects large-scale copy number variation regions on many chromosomes. The rest is in turn CNVnator, CNV-PCC, FREEC, CNV_IFTV and GROM-RD. For ease of viewing the detection results, the present application Figure 5 Color display is used, and gray processing will greatly affect the display effect.

[0164] In order to further analyze the distribution of the six methods on real samples, the number of CNVs detected by each method on each chromosome is counted in Figure 6In the middle, the number of copy number variations detected by each method on three samples is shown, as well as the number of copy number variations overlapping between two methods. It is found that CNV-PCC has the most average overlapping copy number variations, which indicates that CNV-PCC has a high consistency with other methods. CNVnator and FREEC follow, and the number of overlapping copy number variations between the two methods is high. CNV-LOF detects the most copy number variations in each sample, but the number of overlapping copy number variations with other methods is less. GROM-RD detects the least copy number variations.

[0165] Since real samples lack complete ground truth like simulated data, in order to make a reasonable evaluation, the overlapping density score (ODS) of each method is further calculated as a reliable measure. The formula for calculating ODS is as follows:

[0166]

[0167] Where L m represents the average overlapping length of a method with other methods, L sum represents the total length of copy number variations detected by the method. The overlapping region between different methods is considered as a true positive result. The first term of the formula can be regarded as sensitivity, and the second term can be regarded as accuracy. The method with higher ODS has better performance. The ODS values of the five methods are shown in Table 1.

[0168] Table 1

[0169] CNVnator FREEC GROM-RD CNV-LOF CNV_IFTV CNV-PCC NA19238 15772 15469 9141 3538 10468 16044 NA19238 17291 15649 11018 2557 10955 16391 NA19238 12658 16878 9125 3541 13996 17919

[0170] As can be seen from Table 1, CNV-PCC has the largest ODS in NA19238 and NA19240, and CNVnator has the largest ODS in NA19239. For the average ODS of the three samples, CNV-PCC has the largest value, followed by FREEC, CNVnator, CNV_IFTV, GROM-RD, and CNV-LOF. Therefore, it can be concluded that the proposed method is effective and reliable in the application of real data.

[0171] In addition, the embodiment of the present application also provides an electronic device, comprising:

[0172] One or more processors; and

[0173] A memory associated with the one or more processors, the memory being configured to store program instructions, the program instructions, when read and executed by the one or more processors, performing the optical information encoding method disclosed in the above embodiment.

[0174] Wherein, as Figure 7As shown, the computer device 12 is in the form of a general- purpose computer device. The components of the computer device 12 can include, but are not limited to, one or more processors or processing units 16, a system memory 28, and a bus 18 that couples various system components, including the system memory 28 to the processing unit 16. The bus 18 represents one or more of any of several bus structures, including a memory bus or memory controller, a peripheral bus, a graphics accelerator bus, and a local bus using any of a variety of bus architectures. By way of example, these architectures include Industry Standard Architecture (ISA) bus, Micro Channel Architecture (MCA) bus, Enhanced ISA bus, Video Electronics Standards Association (VESA) local bus, and Peripheral Component Interconnect (PCI) bus.

[0175] The computer device 12 typically includes a variety of computer system readable media. Such media can be any available media that is located either internally or externally to the computer device 12, including both volatile and nonvolatile media, removable and non-removable media.

[0176] The system memory 28 can include computer system readable media in the form of volatile memory, such as random access memory (RAM) 30 and / or cache memory 32. The computer device 12 can further include other removable / non-removable, volatile / non-volatile computer system storage media. By way of example only, a storage system 34 can be provided for reading from and writing to non-removable, non-volatile magnetic media (not shown and typically called a "hard drive"). Although not specifically shown, a magnetic disk drive can also be used for reading from and writing to a removable, non-volatile magnetic disk (e.g., a "floppy disk"), and an optical disk drive can be used for reading from or writing to a removable, non-volatile optical disk (such as a CD-ROM, DVD-ROM or other optical media). In these instances, each can be connected to the bus 18 by one or more data media interfaces. The memory 28 can include at least one program product having a set (e.g., at least one) of program modules that are configured to carry out the functions of embodiments of the application.

[0177] The program / utility 40, having a set (at least one) of program modules 42, can be stored in memory 28 by way of example, and not limitation, as well as an operating system, one or more application programs, other program modules, and program data, each or some combination thereof, can include implementation of a networking environment. The program modules 42 generally carry out the functions and / or methodologies of embodiments of the application as described herein.

[0178] Computer device 12 can also communicate with one or more external devices 14 such as a keyboard, a pointing device, a display 24, etc.; one or more devices that enable a user to interact with computer device 12; and / or one or more devices that enable computer device 12 to communicate with one or more other computing devices. Such communication can occur via input / output (I / O) interface(s) 22. Still yet, computer device 12 in this example can communicate with one or more networks, such as a local area network (LAN), a wide area network (WAN), and / or the Internet, as described above, via network adapter 20. As depicted, network adapter 20 communicates with the other components of computer device 12 via bus 18. It should be appreciated that although not shown, other hardware and / or software modules could be used in conjunction with computer device 12. Examples, include, but are not limited to, microcode, device drivers, redundant processing units, external disk drive arrays, RAID systems, tape drives, and data archival storage systems, etc.

[0179] Processing unit 16 can execute instructions for various functions and data processing by running programs stored in system memory 28.

[0180] Each of the embodiments described in the specification is described in a progressive manner, and the same or similar parts between the embodiments can be referred to each other. Each embodiment focuses on the difference from other embodiments. In particular, for the system or system embodiments, since they are basically similar to the method embodiments, they are described more simply, and the relevant parts can be referred to the part of the method embodiments. The above-described system and system embodiments are merely illustrative, and the units described as separate components can be or can not be physically separated, and the components displayed as units can be or can not be physical units, i.e., they can be located in one place or distributed to multiple network units. Part or all of the modules can be selected to achieve the purpose of the embodiment according to actual needs. Those skilled in the art can understand and implement it without creative labor.

[0181] The above describes the technical solutions provided by the present application in detail, and the principles and implementation modes of the present application are described by applying specific examples. The above description of the embodiments is only to help understand the method and its core idea of the present application; at the same time, for those skilled in the art, according to the idea of the present application, the specific implementation mode and application range will be changed. In summary, the content of the specification should not be understood as a limitation of the present application.

[0182] All the optional technical solutions described above can be combined to form optional embodiments of the present application, which will not be described one by one here.

[0183] The above only describes the preferred embodiments of the present application and is not intended to limit the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A method of detecting copy number variations, characterized by, The method comprises the following steps: dividing a genome into genome bins, generating an information profile of the genome, wherein the information profile comprises read depth signals and alignment qualities of each of the genome bins; globally segmenting the genome according to the information profile, and locally segmenting at least part of the genome after the global segmentation, to obtain gene segments and read depth signals and alignment qualities of the gene segments; the global segmentation comprises: determining a set of genome bins with continuous read depth signals; comparing the average value of the read depth signals of a genome bin with the average value of the read depth signals of the remaining genome bins according to a maximum statistic; if the comparison result meets a variation threshold condition, it is determined that there is a genome bin with a read depth signal corresponding to a variation in the continuous genome bins, and the continuous genome bins are divided into a gene segment; the local segmentation comprises: obtaining a preset segmentation length; dividing part of the gene segments into a plurality of continuous and non-overlapping gene segments according to the segmentation length; and calculating the read depth signals and alignment qualities of the gene segments; using the read depth signals and alignment qualities of the gene segments as classification features, calculating an anomaly score of the gene segments, and identifying a copy number variation region of the genome, comprising: expressing the read depth signals and alignment qualities of the gene segments as a standardized matrix, wherein the read depth signals and alignment qualities of one of the gene segments in the standardized matrix are taken as a data sample; calculating a covariance matrix according to the standardized matrix; calculating eigenvalues and eigenvectors of the covariance matrix; calculating a projection distance of each of the data samples on the eigenvectors as an anomaly score; determining an abnormal sample in the data sample according to the anomaly score and a set threshold value; determining a baseline according to the read depth signals of the gene segment corresponding to the abnormal sample, and declaring the copy number variation region.

2. The method of claim 1, wherein, The method further comprises the following steps: obtaining a test sample and a reference sample of the genome; aligning the test sample and the reference sample to obtain an alignment result; dividing the genome into the genome bins according to the alignment result; calculating the read depth signals and alignment qualities in the genome bins to generate the information profile.

3. The method of claim 2, wherein, The method further comprises the following steps: calculating original read depth signals in the genome bins; performing normalization processing on the original read depth signals.

4. The method of claim 1, wherein, The method further comprises the following steps: performing denoising processing on the read depth signals in the gene segments.

5. A copy number variation detection device, comprising: a file generation module configured to divide a genome into genome bins, and generate an information profile of the genome, wherein the information profile comprises read depth signals and alignment qualities of each of the genome bins; a segmentation module configured to globally segment the genome according to the information profile, and locally segment at least part of the genome after the global segmentation, to obtain gene segments and read depth signals and alignment qualities of the gene segments; The global segmentation comprises: determining a set of continuous genomic bins of read depth signals; comparing the average of read depth signals of the genomic bins with the average of read depth signals of the remaining genomic bins according to the maximum statistic; If the comparison result meets the variation threshold condition, it is determined that there is a variation in the read depth signals corresponding to the genomic bins in the continuous genomic bins, and the continuous genomic bins are divided into a gene segment; The local segmentation comprises: obtaining a preset segmentation length; dividing part of the gene segment into a plurality of continuous and non-overlapping gene segments according to the segmentation length; and calculating the read depth signals and alignment quality of the gene segments; The detection module is configured to calculate an anomaly score of the gene segment by taking the read depth signals and alignment quality of the gene segment as classification features, and identify a copy number variation region of the genome, comprising: The read depth signals and alignment quality of the gene segment are represented as a standardized matrix, and the read depth signals and alignment quality of one gene segment in the standardized matrix are taken as a data sample; a covariance matrix is calculated according to the standardized matrix; eigenvalues and eigenvectors of the covariance matrix are calculated; a projection distance of each data sample on the eigenvectors is calculated as an anomaly score; an abnormal sample in the data sample is determined according to the anomaly score and a set threshold; a baseline is determined according to the read depth signals of the gene segment corresponding to the abnormal sample, and the copy number variation region is declared.

6. An electronic device, comprising: Comprise: One or more processors; And A memory associated with the one or more processors, the memory is used to store program instructions, the program instructions are read and executed by the one or more processors to perform the method in any one of claims 1-4.

7. A computer readable medium having stored thereon a computer program, wherein, The program is executed by the processor to realize the method in any one of claims 1-4.

Images