Gene copy number variation detection method and device, equipment, storage medium

Abstract

The application belongs to the technical field of biological information detection, and discloses a gene copy number variation detection method, according to the confidence interval determination of copy abnormal sites, clustering, determining the copy abnormal sites belonging to the high copy variation sub-class and the RD value being greater than the corresponding confidence interval and the copy abnormal sites belonging to the low copy variation sub-class and the RD value being less than the corresponding confidence interval as the copy variation sites, merging the copy variation sites with adjacent positions and the same variation type to obtain more accurate CNV fragments, separating the interval size of the RD site and the interval size of the CNV fragment, so that the resolution of the RD site can be improved while the stability of the CNV fragment is ensured, and the unsupervised machine learning method is used for the clustering classification of the copy abnormal sites and the comparison with the confidence interval, so that the copy variation sites are double-recognized, the recognition accuracy of the copy variation sites is further improved, and the accuracy, sensitivity and stability of the CNV detection are improved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of biological information detection, and particularly relates to a gene copy number variation detection method and device, equipment and a storage medium. BACKGROUND

[0002] The macrogenomic clinical detection technology is a new type of clinical technology for identifying and diagnosing microbial infection from the perspective of genetic material by using second-generation high-throughput sequencing. The microbial sequences (i.e. non-human sequences) contained in high-throughput sequencing are relatively limited, and more than 90% of the sequence content is human sequence. However, in the current clinical application, only a small amount of microbial sequence is used for infectious identification, and a large amount of human sequence still lacks application and analysis. The macrogenomic technology uses free nucleic acid for microbial identification, and tumor cells are easily detected in free nucleic acid due to their high metabolic rate. Free nucleic acid has become an important experimental sample for tumor screening and examination. If macrogenomics is used for infectious detection while performing early tumor screening detection, not only can the current disease problem be solved, but also early intervention and treatment can be performed on the positive results of early tumor screening, which has great significance for human health. However, so far, there are few algorithms and software that combine macrogenomics and tumor gene copy number variation (CNV) detection. Therefore, the technology limitation prevents the use of macrogenomic infection detection methods for simultaneous tumor screening.

[0003] Currently, the mainstream theories for CNV detection include Read-Pair (RP) method, Split-read (SR) method, Read-Depth (RD) method and Assembly (AS) method. Among them, RP is the earliest algorithm, which uses the distribution of double-end sequencing insert fragment length to detect CNV, also known as pair end mapping (PEM) method. When the insert fragment length is too long or too short, it represents that the genome has undergone structural variation. The SR method uses reads that can be aligned on one end and cannot be aligned on the other end to identify CNV. The other end cannot be aligned, which may indicate CNV. By splitting the individual reads, they can be correctly aligned to the reference genome, and the splitting point is the CNV breakpoint. The RD method uses the correlation between copy number and corresponding region sequencing depth for analysis. The basic model is that the sequencing depth of the deletion region is relatively low, and the sequencing depth of the insertion region is relatively high. The AS method uses short sequences obtained by sequencing to assemble, compares the assembled contig with the reference genome, and determines the region where structural variation has occurred.

[0004] Because the macro gene detection is a 50-75 bp second generation single end sequencing technology, it does not conform to the current mainstream software, and the site difference introduced by the regional population site and different experimental conditions cannot be well corrected, so that its resolution and results are not good in the current software, and there is still no algorithm that can be used for macro gene CNV detection. According to the mainstream theory of CNV detection, the RP and SR methods depend on the double end sequencing technology, which is not suitable for macro gene data, and the relative algorithm is not accurate enough. The AS method depends on the flux and sequencing coverage technology, and the coverage of the macro gene is quite different from the genome assembly technology. This method cannot be run on the macro gene data. The three kinds of theoretical methods do not conform to the macro gene sequencing data, and only a small amount of software under the RD method can be applied.

[0005] However, the RD method depends on the flux and sequencing depth technology, and needs higher and more stable depth changes to identify CNV, so that in the data of macro gene with low sequencing depth, the method will introduce many false positive sites. And the analysis model of the traditional RD method uses the same RD resolution and CNV resolution, and too small resolution will lead to too strong data dispersion, too high false positive, and too large resolution will lead to CNV averaging, resulting in false negative results, and the edge position of CNV may also form a transition type due to the interval coverage in RD calculation, affecting the recognition and judgment of CNV, so that the accuracy and sensitivity of CNV detection are not enough. SUMMARY

[0006] The purpose of the present application is to provide a gene copy number variation detection method and device, equipment and storage medium, which can use different resolutions for RD sites and CNV fragments, improve the resolution of RD sites while ensuring the stability of CNV fragments, and improve the accuracy and sensitivity of CNV detection.

[0007] The first aspect of the present application discloses a gene copy number variation detection method, comprising:

[0008] The sequencing depth value of a plurality of specified sites is determined from the sequencing data of the sample to be tested; wherein each specified site corresponds to a confidence interval and an RD set mean value;

[0009] If the sequencing depth value of any specified site is not located in the confidence interval corresponding to the specified site, the specified site is determined as a copy abnormal site;

[0010] According to the RD set mean value and the sequencing depth value corresponding to each copy abnormal site, the copy number value of each copy abnormal site is calculated; wherein the copy number value and the sequencing depth value are in a positive correlation relationship;

[0011] copy number values of all the copy abnormal sites are spatially clustered and classified to obtain a normal copy class and a copy variation class, the copy variation class including two copy variation sub-classes, namely a high copy variation sub-class and a low copy variation sub-class;

[0012] copy abnormal sites belonging to the high copy variation sub-class and having sequencing depth values greater than the corresponding confidence intervals, and copy abnormal sites belonging to the low copy variation sub-class and having sequencing depth values less than the corresponding confidence intervals, are determined as copy variation sites respectively;

[0013] copy variation sites that are adjacent in position and belong to the same copy variation sub-class are merged to obtain copy variation fragments.

[0014] The second aspect of the present application discloses a gene copy number variation detection device, comprising:

[0015] a depth determination unit configured to determine sequencing depth values of a plurality of specified sites from sequencing data of a sample to be tested; wherein each of the specified sites corresponds to a confidence interval and an RD set mean value;

[0016] an abnormality detection unit configured to determine a specified site as a copy abnormal site when the sequencing depth value of the specified site is not located within the confidence interval corresponding to the specified site;

[0017] a calculation unit configured to calculate copy number values of each of the copy abnormal sites according to the RD set mean value and the sequencing depth value corresponding to each of the copy abnormal sites; wherein the copy number values and the sequencing depth values are in a positive correlation relationship;

[0018] a clustering unit configured to spatially cluster and classify copy number values of all the copy abnormal sites to obtain a normal copy class and a copy variation class, the copy variation class including two copy variation sub-classes, namely a high copy variation sub-class and a low copy variation sub-class;

[0019] a variation determination unit configured to determine copy abnormal sites belonging to the high copy variation sub-class and having sequencing depth values greater than the corresponding confidence intervals, and copy abnormal sites belonging to the low copy variation sub-class and having sequencing depth values less than the corresponding confidence intervals, as copy variation sites respectively;

[0020] a merging unit configured to merge copy variation sites that are adjacent in position and belong to the same copy variation sub-class to obtain copy variation fragments.

[0021] The third aspect of the present application discloses an electronic device comprising a memory storing executable program codes and a processor coupled with the memory; the processor invokes the executable program codes stored in the memory to execute the gene copy number variation detection method disclosed in the first aspect.

[0022] The fourth aspect of the present application discloses a computer readable storage medium storing a computer program, wherein the computer program enables a computer to execute the gene copy number variation detection method disclosed in the first aspect.

[0023] The present application has the beneficial effect that the provided gene copy number variation detection method and device, equipment, storage medium, the method mainly comprises: determining the copy abnormal site by using the confidence interval of each specified site, calculating the copy number value of the copy abnormal site for spatial clustering classification, obtaining the normal copy class and the copy variation class, wherein the copy variation class includes two copy variation subclasses, which are high copy variation subclass and low copy variation subclass, at this time the copy abnormal site can be well classified by spatial clustering, and then the automatic identification of high copy variation sites and low copy variation sites can be realized, and the copy abnormal sites belonging to the normal copy class are filtered out to reduce the misjudgment rate; then, the copy abnormal sites belonging to the high copy variation subclass and having a sequencing depth value greater than the corresponding confidence interval, and the copy abnormal sites belonging to the low copy variation subclass and having a sequencing depth value less than the corresponding confidence interval are determined as copy variation sites; then, the copy variation sites with adjacent positions and the same variation type are merged to obtain more accurate CNV fragments.

[0024] Therefore, the present application separates the interval size of the RD site and the interval size of the CNV fragment, uses different resolutions for the RD site and the CNV fragment, improves the stability of the RD site to improve the stability of the CNV fragment, so that the resolution of the RD site can be improved while the stability of the CNV fragment is ensured, and the data structure problems such as unstable large fragment depth and too large unit point dispersion are prevented, and the unsupervised machine learning method is used to cluster and classify the copy abnormal sites and compare them with the confidence interval to double-identify the copy variation sites, which can further improve the identification accuracy of the copy variation sites and improve the accuracy, sensitivity and stability of the CNV detection. BRIEF DESCRIPTION OF DRAWINGS

[0025] The drawings herein show specific examples of the technical solutions of the present application, and constitute part of the specification together with the specific embodiments, for explaining the technical solutions, principles and effects of the present application.

[0026] Unless specifically stated or defined otherwise, the same reference signs in different drawings represent the same or similar technical features, and different reference signs may also be used to represent the same or similar technical features.

[0027] Figure 1 is a flowchart of a gene copy number variation detection method disclosed by an embodiment of the present application;

[0028] Figure 2 is a basic structure diagram of metagenomic sequencing data disclosed by the embodiment of the present application;

[0029] Figure 3 is a detection effect comparison diagram of the embodiment of the present application and other CNV software at the same resolution;

[0030] Figure 4 is a detection effect comparison diagram of the embodiment of the present application and other CNV software at the best resolution;

[0031] Figure 5 is a detection effect comparison diagram of the embodiment of the present application and other CNV software for positive samples;

[0032] Figure 6 is a structural schematic diagram of a gene copy number variation detection device disclosed by the embodiment of the present application;

[0033] Figure 7 is a structural schematic diagram of an electronic device disclosed by the embodiment of the present application.

[0034] Legend:

[0035] 601, depth determination unit; 602, anomaly detection unit; 603, calculation unit; 604, clustering unit; 605, variation determination unit; 606, merging unit; 701, memory; 702, processor. DETAILED DESCRIPTION

[0036] Unless specifically stated or otherwise defined, all technical and scientific terms used herein are of the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. In the case of combining the technical solutions of the present application with the actual scene, all technical and scientific terms used herein can also have meanings corresponding to the purpose of implementing the technical solutions of the present application. The terms "first, second" used herein are only used for distinguishing names and do not represent specific quantities or sequences. The term "and / or" used herein includes any and all combinations of one or more related listed items.

[0037] Unless specifically stated or otherwise defined, "the", "this" used herein refers to the technical features or technical contents mentioned or described before the corresponding position, which can be the same as the technical features or technical contents mentioned, or can be similar to the technical features or technical contents mentioned.

[0038] Without doubt, the technical content or technical features contrary to the purpose of the present application or obviously contradictory should be excluded. In order to facilitate the understanding of the present application, the specific embodiments of the present application will be described in more detail below with reference to the drawings of the specification.

[0039] It should be noted that the gene copy number variation detection method and device disclosed by the application can be applied to low sequencing depth data such as metagenomic data, and can also be applied to normal sequencing high sequencing depth data, such as whole exon sequencing data, whole gene sequencing data, and chip sequencing data. In the embodiments of the application, the metagenomic sequencing data is taken as an example for illustration, and should not be considered as a limitation of the application.

[0040] (I) Baseline establishment process

[0041] Preferably, in the embodiments of the application, before the sequencing data of the to-be-tested sample is subjected to copy number variation detection, a large number of clinical sample data of the same batch or condition and source can be used in combination with a machine learning method and a variance index to filter unstable and unreasonable noise sites, a K-mean algorithm is used to identify centromere, mitochondrial and repetitive sequence sites in the sites, and a variance is used to remove fluctuating sites, and stable sites that can be used for CNV detection are retained, and then a baseline of each stable site is established.

[0042] Specifically, before the sequencing data of the to-be-tested sample is subjected to copy number variation detection, the following steps S01-S07 can be performed:

[0043] S01, obtaining sequencing data of N1 training samples; the sequencing data of each training sample includes sequencing depth values of M1 candidate sites.

[0044] In the application, the site actually refers to a site interval, which refers to a chromosome segment, and the sequencing depth value corresponds to the Read-Depth (RD) data of the chromosome segment. It should be noted that when early screening is performed on a tumor, most CNV fragments related to the tumor usually occur on autosomes, and the gene copy number variation detection method is mainly applied to the detection scene of CNV fragments on autosomes, but is not limited thereto, and in some possible embodiments, the method can also be applied to the detection scene of CNV fragments on sex chromosomes.

[0045] Since an unsupervised clustering method is performed, the obtained sequencing data can only include raw fastq format sequencing data, and does not require clinical verification and follow-up data, for example, the training sample is the result data of CNV detection positive or negative. Therefore, the application scenario of the application can be double-blind, which is more meaningful for early screening. After obtaining the raw fastq format sequencing data of the N1 training samples, further preprocessing work such as filtering, adapter removal, alignment, and deduplication is performed on the raw fastq format sequencing data according to a pre-designed comparison and quality control process.

[0046] Then, according to the preset block interval gradient parameter, the size of the site interval (i.e., the chromosome segment) is determined, and the chromosomes in the sequencing data are divided into M1 chromosome segments according to the size of the site interval, and each chromosome segment corresponds to a candidate site. For example, assuming that the preset block interval gradient parameter is 10 KB, the stable site interval size under this process is 10 KB, and based on this, the chromosomes are divided into M1 chromosome segments, and the size of each chromosome segment is 10 KB.

[0047] After obtaining the M1 candidate sites of each training sample, the sequencing depth values of the M1 candidate sites of the N1 training samples can be preferably subjected to data standardization processing, so that the sequencing depth values of the N1 training samples under the same candidate site that differ greatly in characteristics fall within the set range [0, 1], thereby eliminating the batch differences of different genomes and sample sequencing depths, and reducing the adverse effects of singular characteristic data with excessively large or small values in the sequencing data.

[0048] S02, unsupervised clustering classification is performed on the sequencing depth values of all candidate sites to obtain multiple classification categories.

[0049] In step S02, interval quality control is performed on the sequencing depth values of the N1 x M1 candidate sites, mainly by using the K-mean algorithm to identify a classification tree to unsupervisedly cluster and classify the sequencing depth values of the N1 x M1 candidate sites into multiple classification categories.

[0050] S03, the first RD mean value of each classification category is calculated, and the noise category is identified from the multiple classification categories according to the first RD mean value.

[0051] In step S03, the sequencing depth values of all candidate sites included in each classification category are averaged to obtain the first RD mean value of each classification category, and then based on the first RD mean value, the classification category with an abnormal mean value can be identified, so that the classification category with an abnormal mean value is determined as a noise category and is excluded. Therefore, background noise sites such as mitochondrial sites, repetitive sequence sites, and centromere sites can be identified by a machine learning clustering algorithm, and efficient and detailed noise reduction can be achieved.

[0052] Specifically, the classification category with a first RD mean value less than a first specified threshold value is determined as a centromere category; the classification category with a first RD mean value greater than a second specified threshold value is determined as a repetitive sequence category; the classification category with a first RD mean value greater than a third specified threshold value is determined as a mitochondrial category; and the centromere category, the repetitive sequence category, and the mitochondrial category are all regarded as noise categories.

[0053] The first specified threshold is smaller than the second specified threshold, and the second specified threshold is smaller than the third specified threshold. The first specified threshold, the second specified threshold and the third specified threshold can be preset by the developer according to actual needs, and in some preferred embodiments, the first specified threshold, the second specified threshold and the third specified threshold can also be determined according to the overall sequencing depth average value of the plurality of classification categories.

[0054] For example, after calculating the first RD average value of each classification category, the first RD average value of all classification categories can be further averaged to obtain the second RD average value of all classification categories, and then one tenth of the second RD average value is taken as the first specified threshold, and five times the second RD average value is taken as the second specified threshold, and twenty times the second RD average value is taken as the third specified threshold. That is, assuming that the second RD average value of all classification categories is C U , the first specified threshold is C U / 10, the second specified threshold is C U *5, and the third specified threshold is C U *20.

[0055] S04, the candidate sites included in the noise category are removed to obtain M2 set sites. M2 is less than M1.

[0056] In the clustering result, the candidate sites included in the centromere category are typical noise sites due to high sequencing difficulty, and generally have nothing to do with tumor CNV variation. The candidate sites included in the repeat sequence category have high copy number and unstable genomic location, and have no clear diagnostic significance, so they are also removed as noise sites; the candidate sites included in the mitochondrion category differ greatly in different cells and become important noise points in low sequencing depth data, therefore, after removing these candidate sites, the remaining M2 set sites, including but not limited to set sites located on the Y chromosome, set sites located on the X chromosome and set sites located on the autosome, are used as stable CNV detection targets.

[0057] S05, calculating the RD set average value and variance of each set site under N1 training samples.

[0058] Specifically, the average value of the sequencing depth value (preferably the standardized RD data) corresponding to each set site of N1 training samples is calculated as the RD set average value corresponding to the set site, and there is a difference between the RD set average values of each set site; and the variance of the sequencing depth value corresponding to each set site of N1 training samples is calculated.

[0059] S06, determining all or part of the set sites from the M2 set sites as specified sites.

[0060] After determining the M2 set positions, in some possible embodiments, all the set positions can be designated positions, i.e., stable positions used for reference in subsequent detection of the to-be-tested sample. In the embodiments of the present application, the M2 set positions are preferably filtered again, and a small number of set positions with large variances in the sample dimension are removed. Specifically, after calculating the variance of each set position under the N1 training samples, the set positions with variances less than a specified variance threshold can be determined as designated positions to obtain M3 designated positions, and the set positions with variances greater than or equal to the specified variance threshold are removed.

[0061] It should be noted that M3 is less than or equal to M2. When the variances of the M2 set positions are all less than the specified variance threshold, the M3 designated positions include all the M2 set positions, and when there are some set positions in the M2 set positions with variances greater than or equal to the specified variance threshold, the M3 designated positions include some set positions.

[0062] S07, set the mean and variance according to the RD, and calculate the confidence interval of each designated position.

[0063] After determining the M3 designated positions, the confidence interval of each designated position is calculated by the following formula:

[0064]

[0065] In the formula, μi represents the RD set mean corresponding to the ith designated position, S represents the variance corresponding to the ith designated position, since the data conforms to the normal distribution, i 2 is 1.96.

[0066] Therefore, the confidence interval expression of the ith designated position is

[0067] That is to say, after the calculation of the N1 training samples is completed, each designated position corresponds to a baseline, and the content of the baseline includes the confidence interval and the RD set mean. The confidence interval is used to define the normal value range of the sequencing depth value of the designated position, and the points outside the confidence interval will be identified as abnormal values. By setting ci=95% as the confidence interval of the designated position, a dedicated residual range is established for each designated position, and with the accumulation of sample quantity, the residual range is more accurate and does not depend on known outcome samples.

[0068] ​The steps S01-S07 are implemented, a baseline is established by using the identified relatively stable set points, the differences between different sites are removed by using the baseline, the noise points of uneven depth between sites introduced by the differences of GC content, population genome source, experimental method and quality control method are removed, the heterogeneity of background sites can be well processed, the noise brought by batch effect and local human characteristics can be better removed, and it is ensured that the identified copy number variation is a real change.

[0069] (II) Detection process of the to-be-detected sample / training sample

[0070] As shown in Figure 1 , the embodiment of the present application discloses a gene copy number variation detection method, which comprises the following steps S10-S60:

[0071] S10, determining the sequencing depth values of a plurality of specified sites from the sequencing data of the to-be-detected sample.

[0072] In the present application, the specified site refers to a pre-specified chromosome segment, which can be specified by the developer according to actual needs in some other possible embodiments, and in the present embodiment, it is determined by the M3 specified sites determined in the establishment process of the baseline. In actual sequencing, the to-be-detected sample may not measure all the specified M3 specified sites, so in step S10, the plurality of specified sites of the to-be-detected sample determined includes all or part of the M3 specified sites. Then the baseline content corresponding to each specified site, i.e. the corresponding confidence interval and RD set mean value, can be retrieved according to the determination of each specified site of the to-be-detected sample.

[0073] S20, if the sequencing depth value of any specified site is not located in the confidence interval corresponding to the specified site, determining the specified site as a copy abnormal site.

[0074] In which, the sequencing depth value of each specified site is compared with the confidence interval corresponding thereto, as long as the sequencing depth value is not located in the confidence interval corresponding to the specified site, the specified site is determined as a copy abnormal site, and after traversing a plurality of specified sites, a plurality of copy abnormal sites can be determined. Specifically, if the sequencing depth value is less than the corresponding confidence interval, the specified site is determined as a low copy abnormal site, so as to obtain a low copy abnormal site set; if the sequencing depth value is greater than the corresponding confidence interval, the specified site is determined as a high copy abnormal site, so as to obtain a high copy abnormal site set.

[0075] For example, assuming that the sequencing depth value of the first specified site of the to-be-detected sample is RD1, if , the first specified site is identified as a low copy abnormal site; if , the first specified site is identified as a high copy abnormal site.

[0076] The high and low copy abnormal sites mentioned above can be copy number variation sites or error points under ecological analysis (a = 0.05). A single abnormal site is not very clear. In order to further improve the accuracy, the copy number value of the copy abnormal site identified above can be calculated for further analysis, that is, step S30 is performed.

[0077] S30, according to the RD setting mean value and the sequencing depth value corresponding to each copy abnormal site, the copy number value of each copy abnormal site is calculated. Wherein, the copy number value and the sequencing depth value are positively correlated.

[0078] In this embodiment, the copy number value of each copy abnormal site is calculated by the following formula:

[0079]

[0080] In the formula, CP j represents the copy number value of the jth copy abnormal site, RD j represents the sequencing depth value of the jth copy abnormal site, represents the RD setting mean value corresponding to the jth copy abnormal site.

[0081] It should be noted that the copy number value of each copy abnormal site is not limited to the calculation method of taking the ratio between the sequencing depth value of the copy abnormal site and the RD setting mean value as the copy number value of the copy abnormal site in the above formula (2). In some other possible embodiments, the calculation method of the copy number value can also use other alternative formulas or obvious variant formulas equivalent to the above formula (2), which can also ensure that the copy number value and the sequencing depth value are positively correlated. For example and the like.

[0082] S40, the copy number values of all copy abnormal sites are spatially clustered and classified to obtain a normal copy class and a copy variation class, wherein the copy variation class includes two copy variation subclasses, which are a high copy variation subclass and a low copy variation subclass.

[0083] Specifically, in the embodiment, the sample conclusion is not pre-known, and the copy abnormal sites are identified by using unsupervised clustering. The preselected range K value is from 1 to 10, and the optimal K value is searched by using the cross-validation method. After a large number of sample verifications, the optimal K value is 4. Therefore, in the embodiment, K is set to 4, and the copy values of all copy abnormal sites are clustered by using the K nearest neighbor (KNN) method to obtain 4 classification categories. Then, the average values of the copy values (CP) of the classification categories are calculated, the 4 classification categories are sorted according to the average values of the CP, and the classification categories are identified as a high copy variation subcategory, a high copy discrete category, a low copy discrete category, and a low copy variation subcategory in the order from high to low according to the average values of the CP. The high copy discrete category and the low copy discrete category can be regarded as normal copy categories, and are used to filter out the copy abnormal sites belonging to the normal copy categories and reduce the misjudgment rate.

[0084] S50, the copy abnormal sites belonging to the high copy variation subcategory and having a sequencing depth value greater than the corresponding confidence interval and the copy abnormal sites belonging to the low copy variation subcategory and having a sequencing depth value less than the corresponding confidence interval are determined as copy variation sites, respectively.

[0085] Then, the copy abnormal sites included in the identified high copy variation subcategory are taken as the first intersection set with the high copy abnormal site set obtained by comparing with the confidence interval, and the copy abnormal sites in the first intersection set are determined as high copy variation sites. Similarly, the copy abnormal sites included in the identified low copy variation subcategory are taken as the second intersection set with the low copy abnormal site set obtained by comparing with the confidence interval, and the copy abnormal sites in the second intersection set are determined as low copy variation sites. The high copy variation sites and the low copy variation sites are collectively referred to as copy variation sites.

[0086] By using the above intersection taking method, that is, by using the K nearest neighbor data clustering method to cluster and classify the copy abnormal sites and the method of comparing with the baseline confidence interval, the copy variation sites are double-identified, which can further improve the identification accuracy of the copy variation sites. Meanwhile, the copy number variation sites can be well clustered by using the space clustering, and the automatic identification of the mutation types such as the increase and decrease of the copy number can be realized. The high copy variation sites and the low copy variation sites can be more accurately identified, and the accuracy, sensitivity, and stability of the CNV detection can be improved.

[0087] S60, the copy variation sites having adjacent positions and belonging to the same copy variation subcategory are merged to obtain copy variation fragments.

[0088] wherein high copy variation sites belonging to the same high copy variant class and being adjacent in position are merged to obtain a high copy variation fragment; and low copy variation sites belonging to the same low copy variant class and being adjacent in position are merged to obtain a low copy variation fragment. The high copy variation fragment and the low copy variation fragment are collectively referred to as a copy variation fragment.

[0089] Specifically, copy variation sites of the same variation type (i.e., belonging to the same high copy variant class or low copy variant class) are continuously merged within a chromosome, i.e., chromosome segments of the same variation type and being adjacent in position are merged. In this embodiment, a sliding window is used to scan the multiple copy variation sites. The length of the sliding window should cover at least 2 copy variation sites. In this embodiment, the interval of each site is 10 KB, and thus the length of the sliding window can be set to 20 KB. The moving step of the sliding window is 10 KB, i.e., the sliding window moves by a distance of 1 copy variation site each time. There are always 2 copy variation sites in the window. After each movement of the sliding window, it is determined whether the latter copy variation site in the sliding window belongs to the same copy variant class as the former copy variation site adjacent thereto. If yes, the latter copy variation site is merged with the former copy variation site adjacent thereto.

[0090] For example, assuming that the multiple copy variation sites are A1, A2, A3, …, A8, A9, …, A j , the sliding window first selects {A1, A2}. If A2 belongs to the same copy variant class as A1, the two are merged. Then the sliding window moves by a step of 10 KB to {A2, A3}. If A3 belongs to the same copy variant class as A2, the two are merged. This is repeated until all copy variation sites are traversed.

[0091] In this process, it is assumed that A1 is merged to A8, and then the sliding window moves to {A8, A9}. If it is determined that A9 belongs to a different copy variant class from A8, and the total step of the sliding window moving from A1 does not exceed 100 KB, A9 can be regarded as the first non-copy variation site appearing in the 100 KB step, and A9 is still regarded as a linked site and linked to A8. Then, if it is determined that A 10 belongs to a different copy variant class from A9, it is determined that the current merged copy variation fragment is interrupted, and A 10 is not merged. If it is determined that A 10 belongs to the same copy variant class as A9, A 10Linkage merge with A9. Based on this example, if continuously merged to A j-1 , and it is determined that A j and A j-1 belong to different copy variant classes, it is necessary to determine the distance between A j and the last-appeared non-copy variant site A9. If the distance between A j and A9 exceeds 100 KB, A j may still be considered as a linkage site and be merged with A j-1 ; otherwise, if the distance between A j and A9 does not exceed 100 KB, it is determined that the current-merged copy variant fragment is interrupted.

[0092] It can be seen that, when the positions of multiple copy variant sites of the same variant type are adjacent, the linkage algorithm is used to extend the sites, the copy variant sites of the same variant type and adjacent in position are merged, more accurate CNV fragments are obtained, the interval size of the RD site and the interval size of the CNV fragment are separated, different resolutions are used for the RD site and the CNV fragment, the stability of the RD site is improved to improve the stability of the CNV fragment, so that the resolution of the RD site is improved while the stability of the CNV fragment is ensured, data structure problems such as large fragment depth instability and unit point dispersion are prevented, the copy variant site is double-recognized by using the unsupervised machine learning method to cluster and classify the copy variant site and compare the copy variant site with the baseline confidence interval, and the recognition accuracy of the copy variant site is further improved, and the accuracy, sensitivity and stability of CNV detection are improved.

[0093] Further, all the calculated copy variant fragments can be screened by setting a copy variant fragment threshold, so as to further improve the accuracy. Specifically, in the embodiment of the present application, after step S60 is performed, the following steps S70 or S80 can also be performed:

[0094] S70, if the size of the copy variant fragment reaches a first fragment threshold, the copy variant fragment is determined as a target copy variant fragment.

[0095] S80, if the size of the copy variant fragment reaches a second fragment threshold and is less than the first fragment threshold, the copy variant fragment is determined as a suspected copy variant fragment.

[0096] For example, the first fragment threshold is set as 10 MB, and the second fragment threshold is set as 1 MB, that is, the copy variant fragment with a size of 1 MB and less than 10 MB is a suspected copy variant fragment, and the copy variant fragment with a size of 10 MB is a target copy variant fragment.

[0097] After step S70 is performed, the following steps S71-S72 can also be performed:

[0098] S71, obtain the copy number value of each target copy abnormal site included in the target copy variant fragment, and perform log conversion on the copy number value of the target copy abnormal site by the following formula:

[0099] CP' z = log(CP z + 0.001) (3)

[0100] In the formula, CP' z represents the copy parameter for plotting, CP z represents the copy number value of the zth target copy abnormal site.

[0101] S72, plot according to the order of chromosomal distribution according to the plotting reference copy parameter of each target copy abnormal site.

[0102] By performing steps S71-S72, the copy number value of the target copy abnormal site is log converted to obtain stable and high readability plotting reference copy parameters. The auxiliary line is added to the identified target copy variant fragment to prompt, and the log transformed site copy number is plotted and visualized with the help of the plotting tool.

[0103] After steps S01-S07 are performed, the above steps S10-S80 are detection steps for unknown single test samples, and in actual application, unknown batch test samples can also be implemented. The sample and sample quantity targeted in step S10 are different, and in the detection process of the test sample, the sequencing data of N2 test samples is first obtained; the sequencing depth values of a plurality of specified sites are determined from the sequencing data of each test sample, and then the above steps S20-S80 are implemented. The test function can be realized, and the baseline of each specified site determined in steps S01-S07 is optimized and updated.

[0104] For example, after N2 test samples are implemented according to step S40, the copy abnormal sites in the copy variant class obtained are all retained and used as new specified sites, and each new specified site uses historical sample (at least including N1 training samples) data and N2 test data of the batch to calculate a new baseline again, and stores the new baseline corresponding to each new specified site and the test data of the batch to the database. In order to detect unknown single test samples next time or optimize specified sites and their baselines for the next batch of test sample data.

[0105] The test samples mentioned above can be largely derived from clinical laboratory data, without the need for additional research training sets and clinical information. The advantages of large clinical laboratory sample size and data volume can be translated into the stability of the specified loci and their baselines. Through continuous testing and optimization of the specified loci and their baselines, better compatibility and stability can be achieved for the testing methods, testing populations, analysis methods, and sequencing methods of the testing institutions.

[0106] like Figure 2 As shown, Figure 2 The data shows the GC content and depth of the autosomal, X / Y sex chromosome, and mitochondrial genomes. Abnormalities with unstable GC content are visible on the chromosomes; these are anomalies that introduce false signals during the analysis. The depth data indicates that the overall genome sequencing depth is very low, between 1 and 2, which does not meet the requirements for CNV identification techniques such as assembly in current technologies.

[0107] Figure 3 For comparison of results at the same resolution, the upper half of the image shows other existing CNV identification software, while the lower half shows the results of this invention. It is evident that the upper half of the image shows severe site dispersion at the same resolution, failing to remove differences between different sites. Additionally, there are a few outliers at both the top and bottom, making it impossible to determine whether they are CNVs or background noise. The lower half of the image shows the detection results of this invention. It is clear that the chromosome markers are obvious, making the image more readable. Furthermore, all sites are normally distributed around the normal 2x copy line (the solid line in the middle). The two dashed lines above and below the solid line represent the 3x and 1x copy lines, respectively, indicating no significant CNV variation. This sample is a negative sample.

[0108] Figure 4 This image shows the resolution performance of other existing CNV identification software at optimal resolution. It's evident that the software produces significantly different results for the same sample at different resolutions. The upper half of the image shows numerous discrete points below the haplotype line, each representing a large 1M mutation segment, indicating a positive sample. However, manual inspection reveals that these abnormal sites are mostly low-depth anomalies originating from low-depth data or low-quality sites, leading to false positives from the software. The lower half of the image shows the detection results of this embodiment. It shows no obvious discrete points below the haplotype line, and all sites are stably distributed within the normal range (between the triplotype and haplotype lines), indicating a clearly normal negative sample.

[0109] Figure 5The resolution of the positive sample is shown. The existing other CNV identification software uses the optimal parameters for analysis, and it can be seen that there is a copy number doubling variation in chromosome 13, but there are differences in different sites, especially the sites at the edge of the variation region present discrete excessive points, so that the position of the whole variation interval cannot be confirmed, and the variation multiple is not clear, and only a few points of chromosome 14 present deletion mutation signals, but the number is small, and it is not yet possible to judge whether it is a variation. The detection result of the embodiment of the present application is shown in the following figure. It can be seen that the mutation result of chromosome 13 is very clear, and a large segment of chromosome presents a copy number increase phenomenon, and the copy number center is between the three-fold copy line and the two-fold copy line, that is, the intercellular variation heterogeneity is identified, which indicates that there may be CNV signals caused by tumors rather than genetic diseases, and the deletion mutation of chromosome 14 is obviously abnormal. The machine learning tool of the present application also identifies the mutation of the site and assists in marking.

[0110] As can be seen from the above, the present application has better specificity and stability than other limited methods. The present application can filter low-depth and repetitive sequence sites, and can use a large number of site baselines to remove site depth differences caused by batch effects, experimental methods, quality control means, regional population characteristics, etc. Through linkage calculation, the stability of CNV detection can be improved from the genomic position dimension and plotted, so it can be well applied to CNV identification under the low-depth sequencing data structure of metagenomic sequencing, and can be used for exploration and early screening of tumors and genetic diseases in infected patients, which has great clinical significance.

[0111] Moreover, by separating the interval size of the RD site and the interval size of the CNV fragment, using different resolutions for the RD site and the CNV fragment, and using the improved stability of the RD site to improve the stability of the CNV fragment, the resolution of the RD site can be improved while the stability of the CNV fragment is ensured, preventing problems such as unstable large fragment depth and excessive unit point dispersion, and using an unsupervised machine learning method to cluster and classify the copy abnormal sites and compare them with the confidence interval to double-identify the copy variation sites, which can further improve the identification accuracy of the copy variation sites and further improve the accuracy, sensitivity and stability of CNV detection.

[0112] As shown in Figure 6 The embodiment of the present application discloses a gene copy number variation detection device, which comprises a depth determination unit 601, an abnormality detection unit 602, a calculation unit 603, a clustering unit 604, a variation determination unit 605 and a merging unit 606, wherein,

[0113] The depth determination unit 601 is used for determining the sequencing depth values of a plurality of specified sites from the sequencing data of the sample to be tested; wherein each specified site corresponds to a confidence interval and an RD set mean value;

[0114] The abnormality detection unit 602 is configured to determine that a specified site is a copy abnormal site when a sequencing depth value of any specified site is not located in a confidence interval corresponding to the specified site.

[0115] The calculation unit 603 is configured to calculate a copy number value of each copy abnormal site according to a RD set mean value and the sequencing depth value corresponding to the copy abnormal site, wherein the copy number value and the sequencing depth value are in a positive correlation.

[0116] The clustering unit 604 is configured to perform spatial clustering classification on the copy number values of all copy abnormal sites to obtain a normal copy class and a copy variation class, and the copy variation class includes two copy variation sub-classes, i.e., a high copy variation sub-class and a low copy variation sub-class.

[0117] The variation determination unit 605 is configured to determine the copy abnormal sites belonging to the high copy variation sub-class and having a sequencing depth value greater than a corresponding confidence interval and the copy abnormal sites belonging to the low copy variation sub-class and having a sequencing depth value less than a corresponding confidence interval as copy variation sites, respectively.

[0118] The merging unit 606 is configured to merge the copy variation sites that are adjacent in position and belong to the same copy variation sub-class to obtain a copy variation fragment.

[0119] Optionally, Figure 6 The gene copy number variation detection device shown can further include the following units not shown in the figure:

[0120] The data acquisition unit is configured to acquire sequencing data of a plurality of training samples; the sequencing data includes sequencing depth values of M1 candidate sites.

[0121] The classification unit is configured to perform unsupervised clustering classification on the sequencing depth values of all candidate sites to obtain a plurality of classification categories.

[0122] The mean calculation unit is configured to calculate a first RD mean value of each classification category.

[0123] The noise identification unit is configured to identify a noise category from the plurality of classification categories according to the first RD mean value.

[0124] The noise elimination unit is configured to eliminate the candidate sites included in the noise category to obtain M2 set sites.

[0125] The baseline calculation unit is configured to calculate a RD set mean value and a variance of each set site under a plurality of training samples, determine all or part of the set sites from the M2 set sites as specified sites, and calculate a confidence interval of each specified site according to the RD set mean value and the variance.

[0126] As Figure 7As shown, the embodiment of the present application discloses an electronic device, comprising a memory 701 storing executable program code and a processor 702 coupled with the memory 701;

[0127] The processor 702 calls the executable program code stored in the memory 701 to execute the gene copy number variation detection method described in each of the above embodiments.

[0128] The embodiment of the present application also discloses a computer readable storage medium storing a computer program, wherein the computer program enables the computer to execute the gene copy number variation detection method described in each of the above embodiments.

[0129] The purpose of the above embodiments is to exemplarily reproduce and deduce the technical solutions of the present application, and to completely describe the technical solutions, purposes and effects of the present application, so as to make the public understand the disclosure of the present application more thoroughly and comprehensively, and not to limit the protection scope of the present application.

[0130] The above embodiments are not based on the exhaustive enumeration of the present application, and there can be many other unlisted embodiments. Any replacement and improvement made without violating the concept of the present application is within the protection scope of the present application.

Claims

1. A method for detecting gene copy number variations, characterized in that, include: Sequencing depth values ​​for multiple specified sites are determined from the sequencing data of the sample to be tested; wherein each specified site corresponds to a confidence interval and a set mean value of RD. If the sequencing depth value of any specified site is not within the confidence interval corresponding to that specified site, the specified site is determined to be a copy aberration site; Based on the mean RD setting and sequencing depth value corresponding to each of the aforementioned copy aberration sites, the copy value of each of the aforementioned copy aberration sites is calculated; wherein, the copy value is positively correlated with the sequencing depth value; Spatial clustering is performed on the copy values ​​of all the copy aberration sites to obtain normal copy class and copy mutation class. The copy mutation class includes two copy mutation subclasses, namely high copy mutation subclass and low copy mutation subclass. Copy aberration sites belonging to the high copy variant subclass and having a sequencing depth value greater than the corresponding confidence interval, and copy aberration sites belonging to the low copy variant subclass and having a sequencing depth value less than the corresponding confidence interval, are respectively identified as copy variant sites. Copy mutation sites that are adjacent in location and belong to the same copy mutation subclass are merged to obtain copy mutation fragments.

2. The gene copy number variation detection method as described in claim 1, characterized in that, The confidence interval and the RD set mean corresponding to the specified site are calculated through the following steps: Sequencing data from multiple training samples is obtained; the sequencing data includes sequencing depth values ​​for M1 candidate sites. Unsupervised clustering was performed on the sequencing depth values ​​of all candidate sites to obtain multiple classification categories; Calculate the first RD mean for each of the stated classification categories; The noise category is identified from multiple classification categories based on the first RD mean; Candidate sites included in the noise category are removed to obtain M2 set sites; Calculate the mean and variance of the RD setting for each of the specified sites across multiple training samples; All or some of the M2 designated sites are identified as designated sites; Based on the mean and variance of the RD, calculate the confidence interval for each specified site.

3. The gene copy number variation detection method as described in claim 2, characterized in that, The step of identifying the noise category from multiple classification categories based on the first RD mean includes: The classification category whose first RD mean is less than the first specified threshold is determined as the centromere category; The classification category whose first RD mean is greater than the second specified threshold is determined as the repeat sequence category; The classification category whose first RD mean is greater than the third specified threshold is determined as the mitochondrial category; The centromere category, the repetitive sequence category, and the mitochondrial category are used as noise categories.

4. The gene copy number variation detection method as described in claim 3, characterized in that, After calculating the first RD mean for each of the classification categories, the method further includes: Calculate the second RD mean of multiple classification categories based on the first RD mean of each classification category; Use one-tenth of the mean of the second RD as the first specified threshold; Five times the mean of the second RD is used as the second specified threshold; Twenty times the mean of the second RD is used as the third specified threshold.

5. The gene copy number variation detection method according to any one of claims 2 to 4, characterized in that, The step of calculating the confidence interval for each specified locus based on the mean and variance set by the RD includes: The confidence interval for each specified site is expressed by the following formula: In the formula, S represents the set mean of RD corresponding to the i-th specified site. i 2 This represents the variance corresponding to the i-th specified site.

6. The gene copy number variation detection method according to any one of claims 1 to 4, characterized in that, The step of calculating the copy number of each of the copy aberration sites based on the mean RD setting and sequencing depth value corresponding to each of the copy aberration sites includes: The copy number of each of the aforementioned copy aberration sites is calculated using the following formula: In the formula, CP j RD represents the copy value of the j-th copy aberration site. j The sequencing depth value representing the j-th copy aberration site. This represents the mean value of the RD setting corresponding to the j-th copy aberration site.

7. The gene copy number variation detection method according to any one of claims 1 to 4, characterized in that, After merging copy mutation sites that are adjacent in position and belong to the same copy mutation subclass to obtain a copy mutation fragment, the method further includes: If the size of the copy mutation fragment reaches the first fragment threshold, the copy mutation fragment is determined to be the target copy mutation fragment; If the size of the copy mutation fragment reaches the second fragment threshold but is less than the first fragment threshold, the copy mutation fragment is determined to be a suspected copy mutation fragment.

8. A gene copy number variation detection device, characterized in that, include: A depth determination unit is used to determine the sequencing depth values ​​of multiple specified sites from the sequencing data of the sample to be tested; wherein each specified site corresponds to a confidence interval and an RD set mean value; An anomaly detection unit is used to determine that a specified site is a copy anomaly site when the sequencing depth value of any specified site is not within the confidence interval corresponding to that specified site. The calculation unit is used to calculate the copy number of each of the copy aberration sites based on the mean RD setting and the sequencing depth value corresponding to each of the copy aberration sites; wherein the copy number is positively correlated with the sequencing depth value; Clustering unit is used to spatially cluster and classify the copy values ​​of all the copy aberration sites to obtain normal copy class and copy mutation class. The copy mutation class includes two copy mutation subclasses, namely high copy mutation subclass and low copy mutation subclass. The mutation determination unit is used to determine copy mutation sites as belonging to the high copy mutation subclass and having a sequencing depth value greater than the corresponding confidence interval, and copy mutation sites belonging to the low copy mutation subclass and having a sequencing depth value less than the corresponding confidence interval, respectively. The merging unit is used to merge copy mutation sites that are adjacent in location and belong to the same copy mutation subclass to obtain copy mutation fragments.

9. An electronic device, characterized in that, It includes a memory storing executable program code and a processor coupled to the memory; the processor calls the executable program code stored in the memory to execute the gene copy number variation detection method according to any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program, wherein the computer program causes a computer to perform the gene copy number variation detection method according to any one of claims 1 to 7.

Images