Method for detecting copy number of sequencing gene in target area based on next-generation sequencing
Through UMI technology, Fastp tools, BWA-MEM algorithm, Samtools and LOWESS regression method, NGS detection of unpaired samples was achieved, which solved the accuracy and accessibility issues of CNV detection, simplified the configuration process, and improved the stability and consistency of detection results.
Patent Information
- Application Number
- CN202510758480.6
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-06-09
- Publication Date
- 2025-09-23
AI Technical Summary
Existing NGS-based CNV detection methods rely on paired samples, have complex configurations, large numerical deviations, poor detection performance, and insufficient clinical accessibility and accuracy.
Using UMI technology, Fastp tools, BWA-MEM algorithm, Samtools, and LOWESS regression method, we eliminated the effects of sequencing depth differences and GC content through inter-sample depth normalization, baseline correction, and GC correction, thus achieving unpaired sample detection.
It improves the accuracy and accessibility of gene copy number detection in the target region, reduces dependence on external databases, simplifies the configuration process, and enhances the stability and consistency of detection results.
Smart Images
Figure CN120690290A_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the field of mNGS detection technology, and in particular to a method for detecting the copy number of sequenced genes in a target region based on second-generation sequencing, aiming to improve the accuracy, accessibility and simplicity of gene copy number detection. Background Art
[0002] Gene copy number variation (CNV) refers to the phenomenon of an increase or decrease in the number of DNA copies in certain gene regions in the genome. Changes in gene copy number will affect its expression level, thereby having an important impact on the function and activity of cells. In cancer research, CNV variation is considered an important biomarker that can be used for early diagnosis of tumors, prognosis assessment, and prediction of treatment response. For example, the amplification of the ERBB2 gene is closely related to the aggressiveness, prognosis, and effectiveness of targeted therapy of tumors such as breast cancer and gastric cancer; the amplification of the MET gene is associated with drug resistance and poor prognosis of non-small cell lung cancer; and the loss of the PTEN gene is related to the treatment response of advanced or metastatic breast cancer.
[0003] Currently, commonly used CNV detection technologies include array comparative genomic hybridization (aCGH), fluorescence in situ hybridization (FISH), quantitative PCR (qPCR), and next-generation sequencing (NGS). Compared to other methods, NGS has become one of the primary techniques for CNV detection due to its high throughput, high resolution, and wide range of applications. NGS can detect CNVs across the genome and provide high coverage and accurate detection in targeted regions.
[0004] According to the coverage of the target region, NGS sequencing can be divided into whole genome sequencing (WGS), whole exome sequencing (WES) and targeted region sequencing (Targeted Sequencing). Among them, targeted region sequencing is widely used in clinical research due to its low cost, high depth and ability to directly capture clinically significant variant sites. However, the existing NGS-based CNV detection methods are mainly for WGS and WES detection. There is limited research on methods for detecting gene copy number by NGS target region sequencing, and there is still much room for improvement. The existing methods that can be used for CNV detection by NGS have the following main shortcomings:
[0005] Reliance on paired samples and poor clinical accessibility: Many current CNV detection tools (such as control-freec and exomeCNV) rely on paired samples (e.g., normal control and tumor sample pairing) to correct for target region bias. However, in actual clinical applications, in many cases only a single tumor sample is available for analysis, and paired normal samples are unavailable, which limits the application of these methods.
[0006] Dependence on external databases and complex configuration: Existing methods require not only a reference genome for deep correction of target regions but also multiple external databases, such as genomic repetitive region files and genomic alignment accessibility files. These external dependencies complicate the data analysis process and increase the barrier to entry.
[0007] Large numerical deviations and poor detection performance: Existing tools often exhibit large numerical deviations when testing samples of different sequencing depths, different batches, or known copy numbers. The detection results are not stable, affecting accuracy and reliability.
[0008] In summary, it is necessary to develop and design a method for detecting the copy number of sequenced genes in the target region based on second-generation sequencing to solve the aforementioned technical problems. Summary of the Invention
[0009] The purpose of the present invention is to provide a method for detecting the copy number of sequenced genes in the target region based on second-generation sequencing, which does not rely on paired samples, significantly reduces the dependence on external databases, and can effectively improve the accuracy, accessibility and simplicity of detection.
[0010] The technical solution adopted by the present invention is: a method for detecting the copy number of sequenced genes in a target region based on second-generation sequencing comprises the following steps:
[0011] Step S1, preprocessing of raw FASTQ data to generate BAM files;
[0012] Step S2, calculation of the original depth of the target region: based on the aligned BAM file and the BED file of the target region, the original sequencing depth of each sample in the target region is calculated;
[0013] Step S3, target region depth normalization: normalization makes the sequencing depth of the target region between samples comparable, eliminating the difference in sequencing depth between samples;
[0014] Step S4, target interval depth baseline correction: by collecting baseline samples of white blood cells from multiple healthy individuals, using the same experimental and analytical procedures, a target area depth file of the baseline samples is obtained, and the target area depth of the sample to be tested is corrected to eliminate the interval depth preference caused by probe coverage deviation and repeated sequence factors;
[0015] Step S5, GC correction of target interval depth: performing local weighted smoothing on the target region copy number ratio obtained in the previous step to eliminate the influence of GC content on the depth of the target region of the sample, obtaining the corrected target region copy number ratio, and then calculating the copy number of each target region;
[0016] Step S6, gene interval copy number calculation: calculate the average of the corrected copy numbers of all target regions within the same gene as the final copy number of the gene.
[0017] Preferably, in step S1, the sample FASTQ data is first split using UMI technology, and then the Fastp tool is used to perform data quality control and remove adapters. Then, the BWA-MEM algorithm is used to align the cleaned reads to the human reference genome hg19, and finally Samtools is used for sorting processing to generate a BAM file.
[0018] Preferably, in step S2, the depth command of Samtools is used to calculate the original sequencing depth of each sample in the target region.
[0019] Preferably, the normalization calculation formula in step S3 is:
[0020]
[0021] Among them, ND i D is the sequencing depth of the sample after normalization in the target region. i is the original sequencing depth of the sample in the i-th target region.
[0022] Preferably, in step S4, baseline samples of leukocytes from 20 healthy individuals are collected, and the baseline depth correction formula used is:
[0023]
[0024] Among them, log2ratio i BD is the log2 transformation of the copy number ratio of the target interval i of the sample to be tested after baseline correction, i ND is the copy number ratio of the target interval i of the sample to be tested after baseline correction, i NDB is the sequencing depth of the sample to be tested after normalization in the i-th target region. ik represents the depth of the k-th baseline sample after normalization in the i-th target interval.
[0025] Preferably, in step S5, the copy number ratio log2ratio of the target region obtained in the previous step is calculated using the LOWESS regression method. iPerform local weighted smoothing to obtain the corrected target region copy number ratio corrected_log2ratio i , the calculation formula for calculating the copy number of each target region is,
[0026]
[0027] Among them, CN i is the copy number of the i-th target interval of the sample to be tested after baseline correction.
[0028] Preferably, the calculation formula for the final copy number in step S6 is,
[0029]
[0030] Among them, CN gene is the copy number of the target gene region in the sample to be tested, CN k is the copy number of the kth target region within the gene range of the sample to be tested.
[0031] The advantages and positive effects of the present invention are:
[0032] The present invention provides a method for detecting the copy number of sequenced genes in a target region based on second-generation sequencing. The method does not rely on paired samples, reduces the complex external database configuration process, and eliminates the deviations caused by different samples, different target regions, and different GC contents through effective depth normalization, baseline correction, and GC correction steps, thereby improving the accessibility and accuracy of target region gene copy number detection. First, through normalization, the problem of incomparable depth of samples with different sequencing depths is solved. Secondly, by using healthy human leukocyte samples processed by the same experimental and analytical process as a baseline, the depth deviation in different target regions caused by factors such as probe coverage and repetitive sequence content is solved, and the dependence on paired samples in copy number analysis is removed by the baseline; furthermore, the LOWESS regression method is used to correct the GC deviation, thereby improving the accuracy of gene copy number detection between different batches and libraries. BRIEF DESCRIPTION OF THE DRAWINGS
[0033] Figure 1 This is a flowchart for analyzing the copy number of sequenced genes in the target region of the second-generation sequencing test;
[0034] Figure 2 This is a comparison of interval depth distribution before and after normalization at different sequencing depths;
[0035] Figure 3 is the copy number distribution of the negative control before and after baseline correction;
[0036] Figure 4 This is the depth distribution diagram of negative and positive quality control products before and after GC calibration;
[0037] Figure 5 This is a panoramic view of the copy number distribution in the target area of the negative control product;
[0038] Figure 6 It is a panoramic view of the copy number distribution in the target area of the positive quality control product. DETAILED DESCRIPTION
[0039] In order to further understand the content, features and effects of the present invention, the following embodiments are given to illustrate in detail.
[0040] See Figure 1 The method of the present invention for detecting the copy number of sequenced genes in the target region based on second-generation sequencing includes the following steps: step S1, preprocessing the original FASTQ data to generate a BAM file; step S2, calculating the original depth of the target region; step S3, normalizing the depth of the target region; step S4, baseline correction of the depth of the target interval; step S5, GC correction of the depth of the target interval; step S6, calculating the copy number of the gene interval.
[0041] The following is a detailed description with reference to specific embodiments.
[0042] 1. Sample preparation
[0043] Baseline samples of leukocytes from 20 healthy individuals were collected. In addition, negative and positive reference samples with known gene copy numbers were collected. ERBB2 amplification reference samples were 3, 4, 4.5, and 6 copies. PTEN deletion reference samples were 0.5, 0.8, 1, and 1.5 copies. The negative reference samples contained 2 copies of PTEN and ERBB2, respectively, while the positive control samples contained 1 copy of PTEN and 2 copies of ERBB2.
[0044] 2. Experimental Process
[0045] The above collected samples were prepared experimentally to obtain data for 20 negative samples and 3 replicates of other reference samples. The following experimental steps were used:
[0046] DNA extraction from 20 negative samples: 200 μL of leukocytes were collected and DNA was extracted using the QIAamp DNA MiniKit (Qiagen, Cat: 51304), and quantified using Qubit (Vazyme, Cat: EQ121-02).
[0047] Reference DNA sources: gDNA from a cell line containing ERBB2 CNV amplification (GW-OGTM001) and a human genomic DNA reference (NA12878) were purchased, mixed, and diluted to prepare reference samples with different CNV amplification copy numbers. These were then accurately quantified using ultra-low-frequency digital PCR and high-depth sequencing methods.
[0048] Prelibrary construction: DNA is sheared, end-repaired, and A-added. DNA ligase is used to connect adapters with single-molecule barcode sequences to both ends of the DNA template, and the prelibrary is obtained by PCR amplification.
[0049] Final library construction: The prelibrary is liquid-phase hybridized with a biotin-labeled oligonucleotide probe, and the library bound to the probe is captured and enriched using streptavidin-coated magnetic beads. Finally, PCR amplification is performed using primers with a tag sequence (Index) and a polymerase to obtain the captured library.
[0050] Sequencing: Sequencing was performed using a DNBSEQ-T7 sequencer (BGI) in PE150 mode, with each sample sequenced at 1.6M.
[0051] 3. Data Analysis
[0052] The method of detecting the number of copies of the target region sequenced gene based on the second generation sequencing provided in the present invention is used, that is, according to Figure 1 The data were analyzed according to the process steps given in .
[0053] 1. Preprocessing of raw FASTQ data: First, use UMI technology to split the sample FASTQ data; then, use the Fastp tool to perform data quality control and remove adapters; then, use the BWA-MEM algorithm to align the cleaned reads to the human reference genome hg19; finally, use Samtools for sequencing and generate BAM files.
[0054] in:
[0055] 1) UMI technology: Unique Molecular Identifier (UMI), a molecular tagging technology primarily used in gene sequencing to improve detection accuracy, also has differentiated applications in materials science. Its core principle is to add unique identifiers to DNA molecules or material surfaces, thereby mitigating sequencing errors or enabling precise identification.
[0056] 2) FASTQ data: It is a standard text format for storing biological sequencing data, containing sequence information and base quality scores, and is widely used in second-generation sequencing (NGS) data analysis.
[0057] 3) Fastp tool: It is a tool for quickly processing NGS offline sequences, namely FASTQ data.
[0058] 4) BWA-MEM algorithm: This is an open source tool set developed by Dr. Heng Li, specifically for processing SAM, BAM, and CRAM format files in high-throughput sequencing data. It supports format conversion, sorting, indexing, statistics, and data filtering, and is widely used in basic medical omics and bioinformatics.
[0059] 5) BAM file: It is a binary compressed format of SAM (Sequence Alignment Map) file, mainly used to store the alignment results of high-throughput sequencing data. It is characterized by small file size, high read and write efficiency, and support for random access.
[0060] 2. Calculation of the original depth of the target region: Based on the aligned BAM file and the BED file of the target region, use the depth command of Samtools to calculate the original sequencing depth of each sample in the target region.
[0061] in:
[0062] 1) Samtools: An open-source toolset developed by Dr. Heng Li, specifically for processing SAM, BAM, and CRAM format files in high-throughput sequencing data. It supports format conversion, sorting, indexing, statistics, and data filtering, and is widely used in basic medical omics and bioinformatics.
[0063] 2) BED file: Browser Extensible Data, a text format file used to describe genomic features. It is mainly used to store genome coordinates and annotation information and is widely used in bioinformatics and genomics research.
[0064] 3. Target region depth normalization: To eliminate the differences in sequencing depth between samples, the sequencing depth of the target region is made comparable through the normalization step. The specific formula is as follows:
[0065]
[0066] Among them, ND i D is the sequencing depth of the sample after normalization in the target region. i is the original sequencing depth of the sample in the i-th target region.
[0067] 4. Target interval depth baseline correction: Baseline samples of leukocytes from 20 healthy individuals were collected and the same experimental and analytical procedures were used to obtain a baseline depth file for the target area. The baseline depth correction formula was then used to correct the target area depth of the sample to be tested to eliminate interval depth bias caused by factors such as probe coverage deviation and repeated sequences. The correction formula is as follows:
[0068]
[0069] Among them, log2ratio i BD is the log2 transformation of the copy number ratio of the target interval i of the sample to be tested after baseline correction, i ND is the copy number ratio of the target interval i of the sample to be tested after baseline correction, i NDB is the sequencing depth of the sample to be tested after normalization in the i-th target region. ik represents the depth of the k-th baseline sample after normalization in the i-th target interval.
[0070] 5. Target interval deep GC correction: Use the LOWESS regression method to calculate the copy number ratio log2ratio of the target region obtained in the previous step. i Perform local weighted smoothing to eliminate the influence of GC content on the depth of the target region of the sample and obtain the corrected target region copy number ratio corrected_log2ratio i The copy number of each target region is then calculated using the following formula:
[0071]
[0072] Among them, CN i is the copy number of the i-th target interval of the sample to be tested after baseline correction.
[0073] in:
[0074] 1) GC content: The percentage of guanine (G) and cytosine (C) in a DNA or RNA sequence.
[0075] 2) LOWESS regression: Locally weighted regression is a nonparametric regression method used to address nonlinear trends and fluctuations in data. Its core concept is to fit a low-dimensional polynomial to a subset of data points at each point in the data set and then use weighted least squares to estimate the value of the dependent variable for the data points of the independent variable near that point. The farther away from the point, the smaller the weight, resulting in the regression function value at that point.
[0076] 6. Gene interval copy number calculation: Calculate the average of the corrected copy numbers of all target regions within the same gene as the final copy number of the gene. The calculation formula is as follows:
[0077]
[0078] Among them, CN gene is the copy number of the target gene region in the sample to be tested, CN k is the copy number of the kth target region within the gene range of the sample to be tested.
[0079] In addition, to illustrate the accuracy of the present invention, a control method is introduced in the analysis process. This method removes the depth correction process in steps 4 and 5 above. The calculation formula is as follows:
[0080]
[0081] Then, the corresponding gene copy number in the control method was calculated according to the same calculation formula as in step 6 above.
[0082] 4. Experimental Results
[0083] Depend on Figure 2 As can be seen, after normalization of the raw depth differences at different sequencing depths, the depths of the same target region in negative controls sequenced at 500x and 1000x are consistent. This step effectively addresses the target region differences introduced by different sequencing depths, making the depths of the same target region comparable across samples of different sequencing depths.
[0084] from Figure 3 As can be seen, the copy number data of the negative control before baseline sample depth correction showed large fluctuations in the target region (coefficient of variation 50.75%), while the copy number fluctuation was significantly reduced after baseline correction (coefficient of variation 14.94%). Reducing data fluctuation has a significant positive effect on detecting amplification and deletion of weakly positive samples.
[0085] In addition, after baseline correction, we can see that the target interval depth is related to GC content. Ideally, the target interval depth should be independent of GC content. Therefore, we further corrected the GC content and depth preference of each library by using LOWESS regression technology. We can see that in negative and positive quality control products, the corrected target interval depth is independent of GC content. Figure 4 .
[0086] Figure 5 This is a panoramic view of the copy number distribution in the target area of the negative control product. Figure 6This is a panoramic view of the copy number distribution of the target region of the positive quality control. As can be seen, after normalization of sample sequencing depth, baseline depth correction, and GC depth correction, the whole-genome distribution of negative and positive quality controls is more stable. Furthermore, the known 7-copy amplification of ERBB2 and 1-copy deletion of PTEN in the positive quality control are also stably detected.
[0087] The detection results of all samples in ERBB2 amplification (see Table 1) and PTEN deletion (see Table 2) were further summarized.
[0088] Table 1 ERBB2 amplification reference gene copy number detection statistics
[0089]
[0090]
[0091] Table 2 PTEN deletion reference gene copy number detection statistics
[0092] Gene Sample No. Mutation Type Theoretical copy number Control method The present invention PTEN NC_R1 No amplification deletion 2 1.26 1.43 PTEN NC_R2 No amplification deletion 2 1.48 1.63 PTEN NC_R3 No amplification deletion 2 1.44 1.58 PTEN LOD_PTEN_1_5_R1 Missing 1.5 1.22 1.28 PTEN LOD_PTEN_1_5_R2 Missing 1.5 1.24 1.30 PTEN LOD_PTEN_1_5_R3 Missing 1.5 1.25 1.31 PTEN PC_R1 Missing 1 0.99 1.08 PTEN PC_R2 Missing 1 0.84 0.95 PTEN PC_R3 Missing 1 0.93 1.03 PTEN LOD_PTEN_0_8_R1 Missing 0.8 0.64 0.68 PTEN LOD_PTEN_0_8_R2 Missing 0.8 0.64 0.68 PTEN LOD_PTEN_0_8_R3 Missing 0.8 0.65 0.68 PTEN LOD_PTEN_0_5_R1 Missing 0.5 0.39 0.40 PTEN LOD_PTEN_0_5_R2 Missing 0.5 0.39 0.40 PTEN LOD_PTEN_0_5_R3 Missing 0.5 0.38 0.40
[0093] It can be seen that the detection value of ERBB2 amplification by the present invention is more stable and accurate than that of the control method.
Claims
1. A method for detecting the copy number of sequenced genes in a target region based on second-generation sequencing, characterized by: The following steps are included: Step S1, preprocessing of raw FASTQ data to generate BAM files; Step S2, calculation of the original depth of the target region: based on the aligned BAM file and the BED file of the target region, the original sequencing depth of each sample in the target region is calculated; Step S3, target region depth normalization: normalization makes the sequencing depth of the target region between samples comparable, eliminating the difference in sequencing depth between samples; Step S4, target interval depth baseline correction: by collecting baseline samples of white blood cells from multiple healthy individuals, using the same experimental and analytical procedures, a target area depth file of the baseline samples is obtained, and the target area depth of the sample to be tested is corrected to eliminate the interval depth preference caused by probe coverage deviation and repeated sequence factors; Step S5, GC correction of target interval depth: performing local weighted smoothing on the target region copy number ratio obtained in the previous step to eliminate the influence of GC content on the depth of the target region of the sample, obtaining the corrected target region copy number ratio, and then calculating the copy number of each target region; Step S6, gene interval copy number calculation: calculate the average of the corrected copy numbers of all target regions within the same gene as the final copy number of the gene.
2. The method for detecting the copy number of sequenced genes in a target region based on second-generation sequencing according to claim 1, wherein: In step S1, the sample FASTQ data was first split using UMI technology, and then the Fastp tool was used for data quality control and linker removal. The cleaned reads were then aligned to the human reference genome hg19 using the BWA-MEM algorithm. Finally, Samtools was used for sequencing to generate a BAM file.
3. The method for detecting the copy number of sequenced genes in a target region based on second-generation sequencing according to claim 1, characterized in that: In step S2, the depth command of Samtools is used to calculate the original sequencing depth of each sample in the target region.
4. The method for detecting the copy number of sequenced genes in a target region based on second-generation sequencing according to claim 1, wherein: The normalization calculation formula in step S3 is: Among them, ND i D is the sequencing depth of the sample after normalization in the target region. i is the original sequencing depth of the sample in the i-th target region.
5. The method for detecting the copy number of sequenced genes in a target region based on second-generation sequencing according to claim 1, wherein: In step S4, baseline samples of leukocytes from 20 healthy individuals are collected, and the baseline depth correction formula used is: Among them, log2ratio i BD is the log2 transformation of the copy number ratio of the target interval i of the sample to be tested after baseline correction, i ND is the copy number ratio of the target interval i of the sample to be tested after baseline correction, i NDB is the sequencing depth of the sample to be tested after normalization in the i-th target region. ik represents the depth of the k-th baseline sample after normalization in the i-th target interval.
6. The method for detecting the copy number of sequenced genes in a target region based on second-generation sequencing according to claim 1, wherein: In step S5, the LOWESS regression method is used to calculate the copy number ratio log2ratio of the target region obtained in the previous step. i Perform local weighted smoothing to obtain the corrected target region copy number ratio corrected_log2ratio i , the calculation formula for calculating the copy number of each target region is, Among them, CN i is the copy number of the i-th target interval of the sample to be tested after baseline correction.
7. The method for detecting the copy number of sequenced genes in a target region based on second-generation sequencing according to claim 1, wherein: The calculation formula for the final copy number in step S6 is, Among them, CN gene is the copy number of the target gene region in the sample to be tested, CN k is the copy number of the kth target region within the gene range of the sample to be tested.