Copy number variation detection process method based on sequencing depth

Abstract

The invention relates to the field of biological detection and analysis, and discloses a copy number variation detection process method based on sequencing depth, which comprises the following steps: step S1: obtaining and processing sequencing data of a negative sample and a to-be-detected sample to obtain reference sequencing data and to-be-detected sequencing data; s2, comparing with a reference genome to obtain a reference comparison file and a comparison file to be detected; s3, dividing a main region and a sub-region, and performing GC content correction; s4, performing filtering screening and normalization processing on the sample main region according to the parameter Mappability to obtain parameter information; s5, judging copy number variation of the to-be-detected sample according to the parameter information, and calculating a copy number value; and S6, outputting a result. According to the method, GC content correction is performed by using the non-overlapping main region and the plurality of overlapping sub-regions, so that the problem of tiny deviation of window sequencing depth caused by window division can be effectively avoided, and copy number variation can be detected 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 copy number variation detection process method based on sequencing depth. Background Art

[0002] Copy Number Variation (CNV) refers to a mutation caused by the deletion or amplification of base fragments larger than 1 kbp. It is an important category of genomic structural variation and is closely related to the occurrence of cancer.

[0003] At present, the methods for analyzing copy number variation using sequencing technology mainly include: (1) Paired-end mapping (PEM) based detection of copy number variation: by comparing the distance between the end read pairs on the reference genome and the known insertion fragments to determine the amplification and deletion of the copy number; (2) Split read (SR) based detection of copy number variation: by evaluating the gap fragments in the reference genome to identify the copy number variation; (3) Read depth (RD) based detection of copy number variation: by searching for regions with large changes in read depth in the genome to detect the regions of copy number variation; (4) De novo assemble based (AS) based detection of copy number variation: first, the read sequence of the sample is aligned with the reference genome and assembled to obtain the contigs sequence, and the copy number variation region is identified by comparing the differences between the contigs and the reference genome. Among them, the RD-based method can detect the insertion and deletion of large fragments in complex genomic regions and can accurately detect the copy number, so it has become the main strategy for detecting copy number variation.

[0004] The calculation of sequencing depth is particularly important in RD-based copy number detection methods, and amplification and sequencing during the sequencing process may lead to GC bias, which may lead to deviations in the sequencing data of the genomic region, resulting in deviations in the calculation of sequencing depth. Traditional correction methods usually divide the genome into continuous small windows and use normalized GC content or local weighted regression scatter point smoothing (LOWESS) to correct each small window. However, the division of different windows and the correlation between windows and adjacent windows may lead to slight deviations in GC correction. Summary of the invention

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

[0006] The specific technical scheme of the present invention is:

[0007] A copy number variation detection process method based on sequencing depth, the specific steps are:

[0008] 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;

[0009] 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;

[0010] Step S3: dividing the main region and sub-regions in turn according to the reference comparison file and the comparison file to be tested, and performing GC content correction on each region;

[0011] Step S4: filtering and screening the main area of ​​each sample according to the parameter Mappability, and performing normalization processing to obtain parameter information;

[0012] 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;

[0013] Step S6: Output the result.

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

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

[0016] Preferably, the specific steps of 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 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;

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

[0020] Step S304: perform GC content correction on the RD value of each regional window in turn, and the correction formula is:

[0021]

[0022] in, is the depth value of window i in the region after correction by GC content, 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;

[0023] Step S305: recalculate the corrected depth value of each site according to each region, and the calculation formula is:

[0024]

[0025] 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.

[0026] Step S306: recalculate the RD value of each window in the main area according to the corrected depth value of each site.

[0027] Preferably, the filtering and screening in step S4 is to calculate the parameter Mappability of each window in the main area, and when Mappability ≤ filtering threshold 1 or Mappability ≥ filtering threshold 2, the window is filtered, wherein the calculation formula of the parameter Mappability is as follows:

[0028]

[0029] 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.

[0030] Preferably, the normalization process is as follows: according to the amount of sequencing data of each sample (negative sample and sample to be tested), a normalization factor N is obtained. s ; Secondly, according to the normalization factor N sEach sample was normalized to obtain the RD value;

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

[0032]

[0033] 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;

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

[0035]

[0036] 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.

[0037] Preferably, the parameter information of step S4 is obtained based on multiple samples, including: the RD value of the main area and the depth variance of the main area. The calculation formula of the depth variance of the main area 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 depth mean of the same window in multiple samples.

[0040] Preferably, the determination steps of step S5 are as follows:

[0041] 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;

[0042] 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;

[0043] 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;

[0044] 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;

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

[0046] Step S533: If only one of the depth value or the depth variance of the negative sample does not have a significant difference, the ploidy number 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;

[0047] 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.

[0048] The calculation formula of the copy value is as follows:

[0049]

[0050] 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.

[0051] Compared with the prior art, the present invention has the following beneficial effects:

[0052] The present invention provides a copy number variation detection process method based on sequencing depth. (1) The present invention uses a combination of a main region and a sub-region to jointly perform GC content correction on the depth value of each window of the main region. Compared with only performing GC content correction on 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 test sample and multiple negative samples to jointly evaluate whether the test sample has copy number variation, and can more accurately and precisely detect the copy number variation of the test sample. Compared with other copy number variation detection methods, this method can establish a baseline data by using negative samples in the early stage, and perform multiple sample copy number variation detection on each single sample, which can better and more accurately detect copy number variation. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 Flow chart of a copy number variation detection process method based on sequencing depth in an embodiment of the present invention.

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

[0055] Figure 34 is a flow chart of copy number variation determination in an embodiment of the present invention. DETAILED DESCRIPTION

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

[0057] Embodiment 1:

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

[0059] Step S1: Acquire sequencing data of multiple negative samples and test samples that do not contain copy number variations and perform quality control to obtain processed sequencing data, which are benchmark sequencing data and test sequencing data.

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

[0061] 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;

[0062] In this example, the sequencing data were aligned with the reference genome hg19 using software Bowtie2 and software BWA, respectively, and the obtained results were processed using software Picard to filter out duplicate sequences, thereby obtaining Bowtie2.bam files and bwa.bam files, respectively.

[0063] Step S3: dividing the main region and sub-regions in turn according to the reference comparison file and the comparison file to be tested, and performing GC content correction on each region;

[0064] Both the main region and the sub-region are continuous, non-overlapping base segments with a length of M bp. There is a repeated base segment between the main region and the sub-region, and the repeat length is m bp; wherein, the total number of main regions and sub-regions is M / m, M can be specifically selected as 1000bp to 50000bp, and m can be specifically selected as 100bp to 5000bp. In the present embodiment, M is selected as 1000bp, m is selected as 100bp, and the total number of main regions and sub-regions is 10.

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

[0066] Step S31: dividing a continuous, non-overlapping main region with a length of 1000 bp;

[0067] Step S32: Divide the main region into 9 continuous, non-overlapping sub-regions with a length of 1000 bp, wherein each region partially overlaps with an overlap length of 100 bp;

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

[0069] Step S34: perform GC content correction on the RD value of each regional window in turn, and the correction formula is:

[0070]

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

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

[0073]

[0074] In this embodiment, is the corrected depth of site t, is the corrected depth value of the ith window in the mth region, RD (m,i) is the original depth value of the ith window in the mth region, SD t is the depth value of site t.

[0075] Step S36, recalculating the RD value of each window in the main area according to the corrected depth value of each site.

[0076] The RD value refers to the depth value of each window, which is calculated by counting the sequencing depth of each base in the window and dividing the sum of the sequencing depths in the window by the length of the window; in this embodiment, the window length is 1000bp, and the RD value of the window is the sum of the base sequencing depths of the window divided by 1000.

[0077] In this embodiment, depth correction is performed on each sequencing sample in turn according to the bwa.bam file according to step S3 to obtain the window RD value of each sequencing sample.

[0078] Step S4: filtering and screening the samples to be tested according to the parameter Mappability based on the main area, and normalizing each window of the main area of ​​the samples to be tested according to the reference sample to obtain parameter information;

[0079] Mappability refers to sequence alignment. Genome alignment may result in errors due to homology regions. Therefore, by calculating Mappability, the validity of the data in the window in the main region can be screened.

[0080] The filtering rule is: when Mappability ≤ filtering threshold 1 or Mappability ≥ filtering threshold 2, the window is filtered. 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 of the parameter Mappability is as follows:

[0082]

[0083] Among them, r bwa is the number of sequences aligned using the software bwa in window i, r bowtie2 The number of sequences aligned using Bowtie2 for window i.

[0084] Since different sequencing samples will produce different sequencing amounts during the sequencing process, in order to more accurately detect copy number variation, the present invention normalizes the sequencing depth, that is, the sequencing depth of each sample is placed at the same data amount level for research. The normalization step includes: first, obtaining the normalization factor N of the sequencing data amount of each sample; s , and then each sample (negative sample and test sample) is normalized according to the normalization factor RD.

[0085] The calculation formula of 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 sum 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 RD normalization is as follows:

[0089]

[0090] 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.

[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 the depth variance of the copy number variation may be significantly different from the window depth value and the depth variance of the normal copy number, the window depth value and the window depth variance are an indicator for measuring whether there is a copy number variation in the window; wherein the window depth value is calculated by the normalized RD value, and the calculation formula of each window depth variance is as follows:

[0092]

[0093] Wherein, 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, is the depth mean of the same window in multiple samples.

[0094] Step S5: judging the copy number variation in the sample to be tested according to the negative sample and the sample to be tested, and calculating the copy value;

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

[0096] The calculation formula of the copy value is as follows:

[0097]

[0098] Among them, CN iis the ploidy value of the sample window i 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 The mean sequencing depth of this window in all samples;

[0099] Among them, the value of a is: if it is on the autosome, the value of a is 2; if it is on the sex chromosome, when the Y chromosome is detected, the value of a is 1, and if the Y chromosome is not detected, the value of a is 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 is <1.9 or the copy number value of the sample to be tested is >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 is <0.9 or the copy number value of the sample to be tested is >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 variation are merged and the output results include: sample number, copy number value, whether variation occurs, copy number variation type, and variation region.

[0103] In this example, a sample to be tested was used to perform copy number variation detection 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 No. Whether mutation occurs Copy value Mutation Type Variation region S1 yes 2.98 Amplification Chr17:37843997-37886810

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

Claims

1. A copy number variation detection process method based on sequencing depth, It is 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: dividing the main region and sub-regions in turn according to the reference comparison file and the comparison file to be tested, and performing GC content correction on each region; 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. According to the copy number variation detection process method based on sequencing depth according to claim 1, It is characterized in that The specific steps of GC content correction in step S3 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: recalculate the RD value of each window in the main area according to the corrected depth value of each site.

3. According to the copy number variation detection process method based on sequencing depth of claim 1, It is characterized in that The filtering rule in step S4 is: when the parameter Mappability ≤ filtering threshold 1 or Mappability ≥ filtering threshold 2, the window is filtered, wherein 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.

4. According to the copy number variation detection process method based on sequencing depth according to claim 1, It is characterized in that The normalization processing steps described in step S4 are: first, obtaining a normalization factor Ns according to the sequencing data volume 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.

5. According to the copy number variation detection process method based on sequencing depth according to claim 1, It is 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.

6. According to the copy number variation detection process method based on sequencing depth according to claim 1, It is 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.

7. According to the copy number variation detection process method based on sequencing depth according to claim 1, It is 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