A method for detecting copy number variations based on sequencing depth

Abstract

The present invention relates to the field of biological detection and analysis, and discloses a method for detecting copy number variation based on sequencing depth, comprising the following steps: Step S1: Obtain the sequencing data of a negative sample and a sample to be tested and perform processing to obtain reference sequencing data and sequencing data to be tested; Step S2: Align them with a reference genome respectively to obtain a reference alignment file and an alignment file to be tested; Step S3: Divide the main region and sub-regions, and perform GC content correction; Step S4: Filter, screen and normalize the main region of the sample according to the parameter Mappability to obtain parameter information; Step S5: Judge the copy number variation of the sample to be tested according to the parameter information and calculate the copy number value; Step S6: Output the result. By using a non-overlapping main region and multiple overlapping sub-regions for GC content correction, the present invention can effectively avoid the problem of slight deviation in the window sequencing depth caused by window division problems, and can detect copy number variation more accurately.

Description

Technical Field

[0001] The present invention relates to the field of biological detection and analysis. More specifically, the present invention relates to a method for detecting copy number variation based on sequencing depth. Background Art

[0002] Copy Number Variation (CNV) refers to mutations caused by deletions or amplifications of base fragments greater than 1 kbp, which is an important category of genomic structural variations and is closely related to the occurrence of cancer.

[0003] Currently, the methods for analyzing copy number variation using sequencing technology mainly include: (1) detecting copy number variation based on the paired-end mapping (PEM) method: judging the amplification and deletion of copy number by comparing the distance between the end reads aligned to the reference genome and the known insert fragments; (2) detecting copy number variation based on the split-read (SR) method: identifying copy number variation by evaluating the gaps fragments aligned to the reference genome; (3) detecting copy number variation by the read depth (RD) method: detecting the regions of copy number variation by searching for regions with large changes in read depth in the genome; (4) detecting copy number variation based on the de novo assembly method (De novo assemble based, AS): first aligning the read sequences of the sample with the reference genome and performing assembly to obtain contigs sequences, and identifying copy number variation regions by comparing the differences between the contigs and the reference genome. Among them, the RD-based method can detect large insertions and deletions in complex genomic regions and can accurately detect copy number, so it has become the main strategy for detecting copy number variation.

[0004] In the RD-based copy number detection method, the calculation of sequencing depth is particularly important. However, amplification and sequencing during the sequencing process may cause GC bias, which may in turn cause bias in the sequencing data of genomic regions, resulting in bias in the calculation of sequencing depth. Traditional correction methods usually divide the genome into continuous small windows and perform correction on each small window using the GC content based on normalization or using the locally weighted regression scatterplot smoothing method (LOWESS). However, different window divisions and the correlation between windows and adjacent windows may cause slight bias in GC correction. Summary of the Invention

[0005] Based on the above description, the present invention provides a method for detecting copy number variations based on sequencing depth. By using a non-overlapping main region and multiple overlapping sub-regions for GC content correction, it can effectively avoid the problem of slight deviations in sequencing depth between adjacent windows caused by inconsistent window division, and can more accurately detect copy number variations.

[0006] The specific technical solution of the present invention is as follows:

[0007] A method for detecting copy number variations based on sequencing depth, the specific steps are as follows:

[0008] Step S1: Obtain the sequencing data of multiple negative samples without copy number variations and the sample to be tested, and perform quality control to obtain the processed sequencing data, which are the reference sequencing data and the sequencing data to be tested respectively;

[0009] Step S2: Align the processed reference sequencing data and the sequencing data to be tested with the reference genome respectively to obtain a reference alignment file and a to-be-tested alignment file;

[0010] Step S3: Sequentially divide the main region and sub-regions according to the reference alignment file and the to-be-tested alignment file, and perform GC content correction on each region;

[0011] Step S4: Filter and screen the main region of each sample according to the parameter Mappability, and perform normalization processing to obtain parameter information;

[0012] Step S5: Judge the copy number variations in the sample to be tested according to the parameter information of the negative sample and the sample to be tested, and calculate the copy number value;

[0013] Step S6: Output the result.

[0014] Preferably, in step S2, the software Bowtie2 and the software BWA are used to align the reference sequencing data and the sequencing data to be tested with the reference genome respectively.

[0015] Preferably, the reference genome described in step S2 is the human genome reference sequence. In one embodiment, the hg19 genome is obtained from the NCBI database as the reference genome for alignment.

[0016] Preferably, the specific steps of the GC content correction in step S3 include:

[0017] Step S301: Divide a continuous and non-overlapping main region with a length of M bp;

[0018] Step S302: Divide n continuous and non-overlapping sub-regions with a length of M bp according to the main region, where there is partial overlap between each region, and the overlap length is m bp;

[0019] Step S303: Calculate the RD value of each window in each region;

[0020] Step S304: Successively perform GC content correction on the RD values of each regional window. The correction formula is:

[0021]

[0022] where, is the depth value of window i in this region after GC content correction, and RD i is the original depth value of window i, and RD avg is the average depth of all windows in this region, and RD gc is the average depth of windows with the same GC content in this region;

[0023] Step S305: Recalculate the corrected depth value of each locus according to each region. The calculation formula is:

[0024]

[0025] where, is the depth of locus t after correction, n is the total number of regions, and RD (m,i) is the original depth value of the i-th window in the m-th region, is the depth value of the i-th window in the m-th region after correction, and SD t is the depth value of locus t.

[0026] Step S306: Recalculate the RD value of each window in the main region according to the corrected depth value of each locus.

[0027] Preferably, the filtering and screening described in step S4 is to calculate the parameter Mappability of each window in the main region. When Mappability ≤ filtering threshold 1 or Mappability ≥ filtering threshold 2, then filter this window. Among them, the calculation formula of the parameter Mappability is as follows:

[0028]

[0029] where, r bwa is the number of sequences aligned to window i using software BWA, and r bowtie2 is the number of sequences aligned to window i using software Bowtie2.

[0030] Preferably, the normalization process is as follows: Obtain the normalization factor N according to the sequencing data volume of each sample (negative sample and test sample) s ; Secondly, according to the normalization factor N sNormalize each sample to obtain the RD value;

[0031] The calculation formula for the normalization factor is as follows:

[0032]

[0033] N s is the normalization factor for sample s; n is the number of samples, including the sum of negative samples and samples to be tested; R i and R s are the sequencing amounts of sample i and sample s respectively;

[0034] The calculation formula for the normalized RD value is as follows:

[0035]

[0036] is the depth value after normalization for the t-th window of sample s; N s is the normalization factor for sample s; is the depth value after correction for window t.

[0037] Preferably, the parameter information in step S4 is obtained based on multiple samples, including: the RD value of the main region and the depth variance of the main region. The calculation formula for the depth variance of the main region is:

[0038]

[0039] where n is the number of samples, including the sum of negative samples and samples to be tested; x i is the depth value of the same window in each sample, is the average depth of the same window in multiple samples.

[0040] Preferably, the judgment steps in step S5 are as follows:

[0041] Step S51: Perform a t-test on the depth values and depth variances of the samples to be tested and negative samples respectively to determine whether there are significant differences;

[0042] Step S52: If there are no significant differences in the depth values and depth variances of the samples to be tested and negative samples, it is considered that there is no copy number variation in the sample to be tested;

[0043] Step S53: If there are significant differences in both the depth values and depth variances of the samples to be tested and negative samples, perform a t-test on the depth values and depth variances of the negative samples;

[0044] Step S531: If there are significant differences in both the depth values and depth variances of the negative samples, it is considered that there is no copy number variation in the sample to be tested;

[0045] Step S532: If only one of the depth value or depth variance of the negative sample shows no significant difference, calculate the ploidy number of the sample to be tested. If there is a significant difference, it is considered that the sample to be tested has a copy number variation.

[0046] Step S533: If only one of the depth value or depth variance of the negative sample shows no significant difference, calculate the ploidy number of the sample to be tested. If there is no significant difference, it is considered that the sample to be tested has no copy number variation.

[0047] Step S534: If both the depth value and depth variance of the negative sample show no significant difference, it is considered that the sample to be tested has no copy number variation.

[0048] Among them, the copy number value calculation formula is as follows:

[0049]

[0050] Among them, CN i is the copy number value of window i of the sample to be tested; a is the copy number of the normal genome, and RD i is the sequencing depth of this window i in the sample to be tested; R mean is the average sequencing depth of this window among all samples.

[0051] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0052] A method for detecting copy number variation based on sequencing depth according to the present invention. (1) The present invention uses a combination of a main region and a sub-region to jointly correct the depth value of each window in the main region for GC content. Compared with only correcting the GC content of the main region, this method takes into account the depth correlation of adjacent positions and can more effectively correct the depth of each window. (2) The present invention uses a combination of a sample to be tested and multiple negative samples to jointly evaluate whether the sample to be tested has a copy number variation, and can more accurately and precisely detect the copy number variation of the sample to be tested. Compared with other copy number variation detection methods, this method can establish a reference data using negative samples in the early stage and perform copy number variation detection on each single sample with multiple samples, and can better and accurately detect copy number variations. Brief Description of the Drawings

[0053] Figure 1 is a flowchart of the method for detecting copy number variation based on sequencing depth in an embodiment of the present invention.

[0054] Figure 2 is a flowchart of correcting the depth of the GC content window based on the main region and the sub-region in an embodiment of the present invention.

[0055] Figure 3This is the flowchart for copy number variation judgment in the embodiments of the present invention. Specific embodiments

[0056] The present invention will be further described below in conjunction with embodiments.

[0057] Embodiment 1:

[0058] A method for detecting copy number variation based on sequencing depth, the process is as Figure 1 shown, and the specific steps include:

[0059] Step S1: Obtain the sequencing data of multiple negative samples without copy number variation and the sample to be tested, and perform quality control to obtain the processed sequencing data, namely the reference sequencing data and the sequencing data to be tested.

[0060] Among them, high-throughput sequencing and quality control are performed on 50 negative samples without copy number variation and 1 sample to be tested. The quality control steps include: removing low-quality (Q30 < 85%), adapter sequences, and sequences containing N.

[0061] Step S2: Align the processed reference sequencing data and the sequencing data to be tested with the reference genome respectively to obtain a reference alignment file and a test alignment file;

[0062] In this embodiment, the software Bowtie2 and the software BWA are used to align the sequencing data with the reference genome hg19 respectively, and the software picard is used to process the obtained results to filter out duplicate sequences, and the Bowtie2.bam file and the bwa.bam file are obtained respectively.

[0063] Step S3: Divide the main region and the sub-region in sequence according to the reference alignment file and the test alignment file, and perform GC content correction on each region;

[0064] Both the main region and the sub-region are continuous and non-overlapping base fragments with a length of M bp. There are repetitive base fragments between the main region and the sub-region, and the repetitive length is m bp; among them, the total number of the main region and the sub-region is M / m. M can be specifically selected as 1000bp - 50000bp, and m can be specifically selected as 100bp - 5000bp. In this embodiment, M is selected as 1000bp, m is selected as 100bp, and the total number of the main region and the sub-region is 10.

[0065] In this embodiment, the process of GC content correction is as Figure 2 shown, and the specific steps are as follows:

[0066] Step S31: Divide 1 continuous and non-overlapping main region with a length of 1000bp;

[0067] Step S32: Divide the main region into 9 consecutive and non-overlapping sub-regions each with a length of 1000 bp, where there is partial overlap between each region, and the overlap length is 100 bp;

[0068] Step S33: Calculate the RD value of each window in each region;

[0069] Step S34: Sequentially perform GC content correction on the RD values of each region window, and the correction formula is:

[0070]

[0071] where, is the depth value of window i in this region after GC content correction, RD i is the original depth value of window i, RD avg is the average depth of all windows in this region, RD gc is the average depth of windows with the same GC content in this region.

[0072] Step S35: Recalculate the corrected depth value of each locus according to each region, and the calculation formula is:

[0073]

[0074] In this embodiment, is the depth after correction of locus t, is the corrected depth value of the i-th window in the m-th region, RD (m,i) is the original depth value of the i-th window in the m-th region, SD t is the depth value of locus t.

[0075] Step S36: According to the corrected depth value of each locus, recalculate the RD value of each window in the main region.

[0076] The RD value refers to the depth value of each window, and the calculation method is to count the sequencing depth of each base in this window and divide the total sequencing depth in this window by the length of this window; in this embodiment, the window length is 1000 bp, and the RD value of the window is the total base sequencing depth of this window divided by 1000.

[0077] In this embodiment, sequentially perform depth correction on each sequencing sample according to the bwa.bam file according to Step S3 to obtain the window RD value of each sequencing sample.

[0078] Step S4: Filter and screen the sample to be tested according to the main region based on the parameter Mappability, and perform normalization processing on each window in the main region of the sample to be tested based on the reference sample to obtain parameter information;

[0079] Mappability refers to sequence alignability. Since genomic regions with homology may cause errors in alignment results, calculating Mappability can be used to screen the data validity of windows in the main region.

[0080] The filtering rules are as follows: when Mappability ≤ filtering threshold 1 or Mappability ≥ filtering threshold 2, then filter the window. Filtering threshold 1 can be set to 0.8 - 0.95, and filtering threshold 2 can be set to 1.05 - 1.2. In this embodiment, filtering threshold 1 is set to 0.85, and filtering threshold 2 can be set to 1.15.

[0081] The calculation formula for the parameter Mappability is as follows:

[0082]

[0083] where r bwa is the number of sequences aligned to window i using the software bwa, and r bowtie2 is the number of sequences aligned to window i using the software Bowtie2.

[0084] Since different sequencing samples may generate different amounts of sequencing data during the sequencing process, in order to more accurately detect copy number variations, the present invention normalizes the sequencing depth, that is, places the sequencing depth of each sample at the same data volume level for research. The normalization steps include: first, obtaining the normalization factor N s of the sequencing data volume of each sample, and secondly, performing normalization processing RD on each sample (negative sample and sample to be tested) according to the normalization factor.

[0085] The calculation formula for the normalization factor is as follows:

[0086]

[0087] N s is the normalization factor of sample s; n is the number of samples, including the total of negative samples and samples to be tested, which is 51 in this embodiment; R i and R s are the sequencing amounts of sample i and sample s respectively.

[0088] The calculation formula for the RD normalization process is as follows:

[0089]

[0090] is the depth value of the t-th window of sample s after normalization processing; N s is the normalization factor of sample S; is the depth value of window t after calibration.

[0091] In this embodiment, the parameter information includes the main window depth value and the main window depth variance of each sample. Since the window depth value and depth variance of copy number variation may have significant differences from those of normal copy number, the window depth value and window depth variance are an index to measure whether there is copy number variation in the window. Among them, the window depth value is calculated from the normalized RD value, and the calculation formula for the depth variance of each window is as follows:

[0092]

[0093] where n is the number of samples, including the sum of negative samples and samples to be tested, which is 51 in this embodiment; x i is the depth value of the same window in each sample, and is the average depth of the same window in multiple samples.

[0094] Step S5: Based on the negative samples and samples to be tested, judge the copy number variation in the samples to be tested and calculate the copy value;

[0095] In this embodiment, the main window depth value and main window depth variance of the sample are used to judge the copy number variation. The t-tests are respectively performed on the depth value and depth variance. If the P-value ≤ 0.05, it indicates that there is a significant difference in this window. If the P-value > 0.05, it indicates that there is no significant difference in this window. The copy number variation judgment process is as Figure 3 shown. The specific steps are as follows: If there is a significant difference in the depth values of all samples of the same window and there is no significant difference in the depth value of the negative sample, it is considered that the sample to be tested in this window has copy number variation under the parameter depth value; if there is a significant difference in the depth variances of all samples of the same window and there is no significant difference in the depth variance of the negative sample, it is considered that the sample to be tested in this window has copy number variation under the parameter depth variance; if both the parameter depth value and depth variance are judged to have copy number variation, it is considered that there is copy number variation in this window and the copy value is calculated; if one of the parameter depth value and depth variance is judged to have copy number variation and the other is not, the copy value is calculated. If the copy value has a significant increase or decrease, it is judged that there is copy number variation in this window; if there is no copy number variation in both the parameter depth value and depth variance, it is judged that there is no copy number variation in this window.

[0096] Among them, the calculation formula for the copy value is as follows:

[0097]

[0098] where CN iis the ploidy value of window i of the sample to be tested; a is the copy number of the normal genome; RD i is the depth value of window i of the sample to be tested; is R mean is the average sequencing depth of this window in all samples;

[0099] Among them, the value of a is: if on the autosome, a takes the value of 2; if on the sex chromosome, when the Y chromosome is detected, a takes the value of 1, and when the Y chromosome is not detected, a takes the value of 2.

[0100] In this embodiment, the region where the copy number value significantly increases or decreases is within [-0.1, 0.1], that is, if the normal copy number value is 2, then when the copy number value of the sample to be tested < 1.9 or the copy number value of the sample to be tested > 2.1, it is considered that there is a copy number variation; if the normal copy number value is 1, then when the copy number value of the sample to be tested < 0.9 or the copy number value of the sample to be tested > 1.1, it is considered that there is a copy number variation.

[0101] Step S6: Output the result;

[0102] In this embodiment, all regions with copy number variations are merged and the output result includes: sample number, copy number value, whether there is a variation, copy number variation type, and variation region.

[0103] In this embodiment, 1 sample to be tested is used to detect copy number variations according to the above steps, and the detection results are shown in Table 1.

[0104] Table 1, Copy number variation detection results of 1 sample to be tested.

[0105] Sample number Whether mutation occurred Copy number Mutation type Mutation region S1 Yes 2.98 Amplification Chr17:37843997-37886810

[0106] According to the specific steps of the present invention, copy number amplification is detected in the sample to be tested S1, and the amplification region is Chr17:37843997 - 37886810.

Claims

1. A copy number variation detection process method based on sequencing depth, characterized in that: The specific steps are as follows: Step S1: obtaining sequencing data of multiple negative samples and test samples without copy number variation and performing quality control to obtain processed sequencing data, which are respectively benchmark sequencing data and test sequencing data; Step S2: aligning the processed benchmark sequencing data and the sequencing data to be tested with the reference genome respectively to obtain a benchmark alignment file and a alignment file to be tested; Step S3: Divide the main region and sub-regions in turn according to the reference comparison file and the comparison file to be tested, and perform GC content correction on each region. The specific steps include: Step S31: Divide a continuous and non-overlapping main region with a length of M bp; Step S32: Divide the main region into n continuous and non-overlapping sub-regions with a length of M bp, wherein each region partially overlaps with an overlap length of m bp; Step S33: Calculate the RD value of each window in each area; Step S34: perform GC content correction on the RD value of each regional window in turn, and the correction formula is: in, is the depth value of window i in the region after GC content correction, RD i is the original depth value of window i, RD avg is the mean depth of all windows in the area, RD gc is the mean depth of windows with the same GC content in the region; Step S35: recalculate the corrected depth value of each site according to each region, and the calculation formula is: in, is the corrected depth of site t, n is the total number of regions, RD (m,i) is the original depth value of the i-th window in the m-th region, is the corrected depth value of the i-th window in the m-th region, SD t is the depth value of site t; Step S36: recalculating the RD value of each window in the main area according to the corrected depth value of each site; Step S4: filtering and screening the main area of ​​each sample according to the parameter Mappability, and performing normalization processing to obtain parameter information; Step S5: determining the copy number variation in the sample to be tested according to the parameter information of the negative sample and the sample to be tested, and calculating the copy number value; Step S6: Output the result.

2. The copy number variation detection process method based on sequencing depth according to claim 1, characterized in that: The filtering in step S4 is to calculate the parameter Mappability of each window in the main area. When the parameter Mappability is less than or equal to the filtering threshold 1 or the parameter Mappability is greater than or equal to the filtering threshold 2, the window is filtered. The calculation formula of the parameter Mappability is as follows: Among them, r bwa The number of sequences aligned using the BWA software for window i, r bowtie2 The number of sequences aligned using Bowtie2 for window i.

3. The copy number variation detection process method based on sequencing depth according to claim 1, characterized in that: The normalization processing steps described in step S4 are: first, obtaining a normalization factor Ns according to the sequencing data amount of the negative sample and the sample to be tested; second, performing normalization processing on all samples according to the normalization factor to obtain an RD value; The calculation formula of the normalization factor is as follows: N s is the normalization factor of sample s; n is the number of samples, including the sum of negative samples and samples to be tested; R i and R s are the sequencing amounts of sample i and sample s respectively; The calculation formula of normalized RD value is as follows: is the normalized depth value of the t-th window of sample s; N s is the normalization factor of sample s; is the corrected depth value of window t.

4. The copy number variation detection process method based on sequencing depth according to claim 1, characterized in that: The parameter information in step S4 includes the RD value of the main area and the depth variance of the main area. The depth variance calculation formula of the main area is: Where n is the number of samples, including the sum of negative samples and samples to be tested; x i is the depth value of the same window in each sample; is the depth mean of the same window in multiple samples.

5. The copy number variation detection process method based on sequencing depth according to claim 1, characterized in that: The step of determining the copy number variation in step S5 is: Step S51: Perform a t-test on the depth values ​​and depth variances of the tested sample and the negative sample to determine whether there are significant differences; Step S52: If there is no significant difference in the depth value and depth variance between the test sample and the negative sample, it is considered that the test sample does not have copy number variation; Step S53: If there are significant differences in the depth value and depth variance between the sample to be tested and the negative sample, a t-test is performed on the depth value and depth variance of the negative sample; Step S531: If there are significant differences in both the depth value and the depth variance of the negative sample, it is considered that there is no copy number variation in the sample to be tested; Step S532: If only one of the depth value or the depth variance of the negative sample does not have a significant difference, the copy number value of the sample to be tested is calculated. If there is a significant difference, it is considered that the sample to be tested has a copy number variation; Step S533: If only one of the depth value or the depth variance of the negative sample does not have a significant difference, the copy number value of the sample to be tested is calculated. If there is no significant difference, it is considered that the sample to be tested does not have a copy number variation; Step S534: If there is no significant difference in the depth value and depth variance of the negative sample, it is considered that there is no copy number variation in the sample to be tested.

6. The copy number variation detection process method based on sequencing depth according to claim 1, characterized in that: The calculation formula of the copy value in step S5 is as follows: Among them, CN i is the copy number of the sample window i to be tested; a is the copy number of the normal genome; RD i is the sequencing depth of the window i in the sample to be tested; R mean is the mean sequencing depth of this window in all samples.

Images