Methods, apparatuses, and media for detecting copy number variations based on amplicon sequencing data

Abstract

This invention relates to a method for detecting copy number variations (CNVs) based on amplicon sequencing data, comprising: aligning preprocessed sequencing data with sequencing data of a reference genome to obtain alignment results; filtering out sequencing data that meets predetermined filtering conditions from the alignment results; obtaining the overall uniformity of the amplicon sequencing regions based on the filtered alignment results to determine the stable amplified target analysis regions; constructing a control set baseline based on the sequencing depth of the stable amplified target analysis regions; determining the breakpoint locations of CNVs in the test sample using a prediction model; calculating the ratio of the sequencing depth of each bed region of the test sample to the mean sequencing depth of the constructed control set baseline bed regions; and comparing the calculated ratio with a predetermined ratio threshold to determine the ploidy of each bed region. This invention enables accurate copy number variation detection results based on amplicon sequencing data.

Description

[0001] This application is a divisional application of patent application No. 202310625489.0, filed on May 30, 2023, entitled "Method, apparatus and medium for detecting copy number variation". Technical Field

[0002] This invention relates generally to bioinformatics processing, and more specifically, to methods, computing devices, and computer storage media for detecting copy number variations based on amplicon sequencing data. Background Technology

[0003] Copy number variation (CNV) refers to the deletion or amplification of a DNA fragment that is at least 1 kbp in length compared to a reference genome. Traditional methods for detecting CNV based on high-throughput sequencing data mainly include: split read (SR), paired-end mapping (PEM), de novo assembly (DA), read depth (RD), and combinations of these four methods.

[0004] Among the traditional methods for detecting copy number variations (CNVs), most are based on read depth strategies. The principle is to determine the CNV based on the significant difference between the number of reads in abnormal regions and the number of reads in normal regions. The amplicon sequencing technology used typically employs multiplex PCR. However, the homogeneity of the products in multiplex PCR library preparation is affected by multiple factors, including template quality, primer concentration and quality, reaction system and conditions, and enzymes. This can lead to significant differences in read depth between different batches of samples, thus affecting the stability of copy number detection.

[0005] In summary, the traditional method for detecting copy number variations based on amplicon sequencing data has the following drawback: the stability of copy number variation detection results for amplicon sequencing kits is not good. Summary of the Invention

[0006] This invention provides a method, computing device, and computer storage medium for detecting copy number variations based on amplicon sequencing data, which can significantly improve the stability of copy number detection results.

[0007] According to a first aspect of the present invention, a method for detecting copy number variations based on amplicon sequencing data is provided. The method includes: aligning preprocessed sequencing data with sequencing data of a reference genome to obtain alignment result data, wherein the sequencing data is obtained via amplicon sequencing technology; filtering out sequencing data that meets predetermined filtering conditions from the alignment result data; obtaining the overall uniformity of the amplicon sequencing region based on the filtered alignment result data to determine a stable amplified target analysis region; constructing a control set baseline based on the sequencing depth of the stable amplified target analysis region; and generating a detection result for the copy number variation of the test sample based on the breakpoint location of the copy number variation in the test sample determined by a prediction model and the constructed control set baseline.

[0008] According to a second aspect of the invention, a computing device is also provided, the device comprising: a memory configured to store one or more computer programs; and a processor coupled to the memory and configured to execute one or more programs to cause the device to perform the method of the first aspect of the invention.

[0009] According to a third aspect of the invention, a non-transient computer-readable storage medium is also provided. This non-transient computer-readable storage medium stores machine-executable instructions that, when executed, cause a machine to perform the method of the first aspect of the invention.

[0010] In some embodiments, determining the target analysis region for stable amplification includes filtering the remaining alignment data based on the obtained overall uniformity to obtain alignment data filtered by uniformity.

[0011] In some embodiments, determining the target analysis region for stable amplification further includes: using the start and end positions of the sequence alignment to the reference genome as a bed region, counting the sequence coverage within each bed region; correcting the sequence coverage within each bed region using the number of aligned sequences of the test sample to obtain the corrected sequence coverage within each bed region; and based on the alignment results data filtered for uniformity, extracting regions in the negative control set where all samples have sequence coverage and the difference in corrected sequence coverage is less than a predetermined correction threshold, to determine the target analysis region for stable amplification.

[0012] In some embodiments, generating detection results for copy number variations in a test sample includes: determining the breakpoint location of copy number variations in the test sample via a prediction model; calculating the ratio of the sequencing depth of each bed region of the test sample to the mean sequencing depth of the baseline bed regions of the constructed control set; comparing the calculated ratio with a predetermined ratio threshold to determine the ploidy of each bed region; statistically analyzing the proportion of each ploidy bed region within the breakpoint region; determining whether the statistically analyzed proportion of each ploidy bed region is greater than or equal to a predetermined proportion threshold; and, in response to determining that the proportion of the bed region with the current ploidy is greater than or equal to the predetermined proportion threshold, determining the ploidy of the stably amplified targeted analysis region as the current ploidy.

[0013] In some embodiments, obtaining alignment result data filtered by uniformity includes: statistically analyzing the T20%X coverage and Fold 80 value of the sample to be tested in order to obtain the overall uniformity of the amplicon sequencing region of the sample to be tested; and filtering the remaining alignment result data based on the comparison between the obtained overall uniformity and a predetermined uniformity threshold to obtain alignment result data filtered by uniformity.

[0014] In some embodiments, sequencing data that meets predetermined filtering conditions includes: sequence data with multiple alignments and sequence data with alignment quality values ​​below a predetermined quality threshold.

[0015] In some embodiments, constructing a control set baseline based on the sequencing depth of a stably amplified target analysis region includes: performing sequencing depth correction on the identified stably amplified target analysis region; and calculating the mean and standard deviation of each bed region in the control set sample based on the sequencing depth of the corrected target region in order to construct the control set baseline.

[0016] In some embodiments, sequencing depth correction for the identified stable amplified target analysis region includes: using alignment result data filtered for uniformity, statistically analyzing the sequence coverage depth of each bed region in the stable amplified target analysis region of the test sample and the control set samples; correcting based on the number of aligned sequences, the number of bed regions, and the average sequencing depth of the test sample and the control set samples; and correcting GC preference and amplicon length preference using a local multinomial regression equation to obtain the standardized depth of the corrected stable amplified target analysis region.

[0017] The summary section is provided to present the chosen concepts in a simplified form, which will be further described in the detailed description below. The summary section is not intended to identify key or principal features of the invention, nor is it intended to limit the scope of the invention. Attached Figure Description

[0018] Figure 1 A schematic diagram of a system for implementing a method for detecting copy number variations based on amplicon sequencing data, according to an embodiment of the present invention, is shown.

[0019] Figure 2 A flowchart of a method for detecting copy number variations based on amplicon sequencing data according to an embodiment of the present invention is shown.

[0020] Figure 3 A flowchart of a method for determining a target analysis region for stable amplification according to an embodiment of the present invention is shown.

[0021] Figure 4 A flowchart of a method for sequencing depth correction for a identified stable amplification target analysis region is shown according to an embodiment of the present invention.

[0022] Figure 5 A flowchart of a method for detecting copy number variations based on amplicon sequencing data according to an embodiment of the present invention is shown.

[0023] Figure 6 A block diagram schematically illustrates an electronic device suitable for implementing embodiments of the present invention.

[0024] In the various figures, the same or corresponding reference numerals indicate the same or corresponding parts. Detailed Implementation

[0025] Preferred embodiments of the invention will now be described in more detail with reference to the accompanying drawings. While preferred embodiments of the invention are shown in the drawings, it should be understood that the invention can be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that the invention will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art.

[0026] The term "comprising" and its variations as used herein signify open inclusion, i.e., "including but not limited to". Unless otherwise stated, the term "or" means "and / or". The term "based on" means "at least partially based on". The terms "one example embodiment" and "one embodiment" mean "at least one example embodiment". The term "another embodiment" means "at least one additional embodiment". The terms "first", "second", etc., may refer to different or the same objects.

[0027] As described above, amplicon sequencing data-based detection of single nucleotide variants (SNPs) and insertions / deletions (InDels) has high accuracy. However, in traditional methods for copy number variation detection based on amplicon sequencing data, the amplicon sequencing technology typically employs multiplex PCR. However, the homogeneity of the products in multiplex PCR library preparation is affected by multiple factors, including template quality, primer concentration and quality, reaction system and conditions, and enzymes. This results in significant differences in read depth between different batches of samples, thus affecting the stability of copy number variation detection.

[0028] To at least partially address one or more of the aforementioned problems and other potential issues, exemplary embodiments of the present invention propose a method for detecting copy number variations (CNVs) based on amplicon sequencing data. This method filters alignment results between preprocessed sequencing data obtained via amplicon sequencing technology and sequencing data from a reference genome to remove sequencing data that meets predetermined filtering criteria. Based on the filtered alignment results, the overall uniformity of the amplicon sequencing regions is statistically analyzed to determine stable amplified target regions. A control set baseline is constructed based on the sequencing depth of the stable amplified target regions. Finally, based on the breakpoint locations of CNVs in the test sample determined by a prediction model and the constructed control set baseline, detection results for CNVs in the test sample are generated. This invention enables accurate CNV detection results based on amplicon sequencing data.

[0029] Figure 1 A schematic diagram of a system for implementing a method 100 for detecting copy number variations based on amplicon sequencing data, according to an embodiment of the present invention, is shown. Figure 1 As shown, system 100 includes a computing device 110 and a sequencing device 130. In some embodiments, the computing device 110 and the sequencing device 130 interact with each other directly or via a network (not shown).

[0030] Regarding sequencing device 130, it is used, for example, to provide a sequencing platform to generate sequencing data about the sample to be tested. The sequencing platform is, for example, but not limited to, sequencing platforms such as Illumina, Life, and BGI. Sequencing read lengths can be, for example, 100 bp, 150 bp, etc. Sequencing device 130 obtains sequencing data about the sample to be tested based on amplicon sequencing technology.

[0031] Regarding the computing device 110, it is used, for example, to detect copy number variations based on amplicon sequencing data. Specifically, the computing device 110 is used to align preprocessed sequencing data with sequencing data from a reference genome to obtain alignment result data; filter out sequencing data that meets predetermined filtering conditions; and obtain the overall uniformity of the amplicon sequencing regions to determine the stable amplified target analysis regions. The computing device 110 is also used to construct a control set baseline based on the sequencing depth of the stable amplified target analysis regions; and to generate detection results for copy number variations in the test sample based on the breakpoint locations of copy number variations in the test sample determined by a prediction model and the constructed control set baseline.

[0032] In some embodiments, the computing device 110 may have one or more processing units, including dedicated processing units such as GPUs, FPGAs, and ASICs, and general-purpose processing units such as CPUs. Additionally, one or more virtual machines may run on each computing device 110. The computing device 110 includes, for example, an alignment result data acquisition unit 112, a sequencing data filtering unit 114, a stable amplification targeted analysis region determination unit 116, a control set baseline construction unit 118, and a copy number variation detection result generation unit 120. The alignment result data acquisition unit 112, sequencing data filtering unit 114, stable amplification targeted analysis region determination unit 116, control set baseline construction unit 118, and copy number variation detection result generation unit 120 may be configured on one or more computing devices 110.

[0033] The alignment result data acquisition unit 112 is used to align the preprocessed sequencing data with the sequencing data of the reference genome in order to obtain alignment result data.

[0034] Regarding the sequencing data filtering unit 114, it is used to filter out sequencing data that meets predetermined filtering conditions from the alignment result data.

[0035] The target analysis region determination unit 116 for stable amplification is used to obtain the overall uniformity of the amplicon sequencing region based on the filtered alignment result data in order to determine the target analysis region for stable amplification.

[0036] Regarding the control set baseline building unit 118, it is used to construct the control set baseline based on the sequencing depth of the target analysis region of stable amplification.

[0037] The copy number variation detection result generation unit 120 is used to generate the copy number variation detection result of the test sample based on the breakpoint position of the copy number variation of the test sample determined by the prediction model and the constructed control set baseline.

[0038] Figure 2 A flowchart of a method 200 for detecting copy number variations based on amplicon sequencing data according to an embodiment of the present invention is shown. It should be understood that method 200 can be used, for example, in... Figure 6 The described electronic device is executed at point 600. It can also be used in... Figure 1 The described computing device 110 performs the operation. It should be understood that method 200 may also include additional actions not shown and / or the actions shown may be omitted, and the scope of the invention is not limited in this respect.

[0039] In step 202, computing device 110 compares the preprocessed sequencing data with the sequencing data of a reference genome to obtain alignment result data, wherein the sequencing data is obtained via amplicon sequencing technology.

[0040] Regarding sequencing data, for example, computing device 110 obtains sequencing data from a sequencing platform (such as, but not limited to, Illumina, Life, BGI, etc.) provided by sequencing device 130, using amplicon sequencing technology. The sequencing read length can be, for example, 100 bp, 150 bp, etc.

[0041] For example, the computing device 110 uses the MEM algorithm of the BWA software to compare and analyze the high-quality usable data generated in step 504 with the HG38 standard reference genome sequence to obtain the alignment position information of each sequence on the reference genome sequence and obtain the alignment result data.

[0042] In step 204, the computing device 110 filters out sequencing data that meets predetermined filtering conditions from the alignment result data.

[0043] Sequencing data that meets predetermined filtering criteria may include, for example, sequence data with multiple alignments and sequence data with alignment quality values ​​below a predetermined quality threshold. For example, sequencing data with multiple alignments and sequences with alignment quality values ​​below a predetermined quality threshold (e.g., 20) may be filtered out from the alignment results.

[0044] In step 206, the computing device 110 obtains the overall uniformity of the amplicon sequencing region based on the filtered alignment result data in order to determine the target analysis region for stable amplification.

[0045] A method for determining stable amplification target analysis regions includes, for example, the following: A computing device 110 filters the remaining alignment result data based on the acquired overall homogeneity to obtain homogeneity-filtered alignment result data; using the start and end positions of sequence alignment to the reference genome as a bed region, the sequence coverage number within each bed region is counted; the sequence coverage number of each bed region is corrected using the number of aligned sequences in the test sample to obtain a corrected sequence coverage number within each bed region; and based on the homogeneity-filtered alignment result data, regions in the negative control set where all samples have sequence coverage and the difference in corrected sequence coverage number is less than a predetermined correction threshold are extracted to determine stable amplification target analysis regions. The following will combine... Figure 3 The specific methods for determining the target analysis region for stable amplification will not be elaborated here.

[0046] Methods for obtaining alignment results data filtered by homogeneity may include, for example,: statistically analyzing the T20%X coverage and Fold 80 value of the test sample to obtain the overall homogeneity of the amplicon sequencing region of the test sample; and filtering the remaining alignment results data based on the comparison between the obtained overall homogeneity and a predetermined homogeneity threshold to obtain alignment results data filtered by homogeneity. It should be understood that T20%X coverage and Fold 80 value can be used as indicators to measure overall homogeneity; test samples with unsatisfactory overall homogeneity will not be included in subsequent analysis of the target region. For example, if the Fold 80 value of a test sample is >3.3, the test sample will be filtered out.

[0047] For example, Table 1 below illustrates the basic data of seven test samples. Among them, the Fold 80 value of test sample CN890 is 3.5, which is greater than the threshold of 3.3, so test sample CN890 needs to be filtered out.

[0048] Table 1

[0049]

[0050] It should be understood that by employing the above methods, the present invention can obtain the overall homogeneity of the sample to be tested, and exclude samples with poor overall homogeneity from subsequent analysis.

[0051] The following uses formulas (1) and (2) to explain the algorithms for calculating T 20%X coverage rate and Fold 80 value.

[0052] T 20%X coverage rate (%) = (number of bases in the amplicon target region greater than 20% of the average depth) / total number of bases in the amplicon target region × 100 (1)

[0053] Fold 80 = Average sequencing depth of amplicon sequencing regions / Depth of amplicon sequencing regions covered by more than 80% of the regions (2)

[0054] In formula (1) above, T 20% x coverage rate (%) represents the coverage percentage of the bases in the amplicon target region that is greater than 20% of the average depth. Generally, the lower the T 20% x coverage rate (%), the better.

[0055] In formula (2) above, the Fold 80 value represents the proportion of the average sequencing depth of the amplicon sequencing region. Generally speaking, the higher the Fold 80 value, the better.

[0056] In step 208, computing device 110 constructs a control set baseline based on the sequencing depth of the stably amplified targeted analysis region.

[0057] Methods for constructing a control set baseline include, for example, performing sequencing depth correction on the identified stable amplified target analysis region; and calculating the mean and standard deviation of each bed region in the control set sample based on the sequencing depth of the corrected target region in order to construct the control set baseline.

[0058] Methods for sequencing depth correction for the identified stable amplification target region include, for example: using alignment data filtered for uniformity, statistically analyzing the sequence coverage depth of each bed region in the target region of the stable amplification in both the test sample and the control set; correction based on the number of aligned sequences, the number of bed regions, and the average sequencing depth of the test sample and the control set; and using local multinomial regression equations to correct GC bias and amplicon length bias to obtain the standardized depth of the corrected stable amplification target region. The following will combine... Figure 4 The specific methods for sequencing depth correction for the identified stable amplified target analysis regions will not be elaborated here.

[0059] In step 210, the computing device 110 generates detection results for the copy number variation of the test sample based on the breakpoint location of the copy number variation of the test sample determined by the prediction model and the constructed baseline of the control set.

[0060] The predictive model is, for example, but not limited to, a Hidden Markov Model or a CBS Model. The predictive model is used to predict the breakpoint locations of copy number variations in the test sample.

[0061] A method for generating detection results for copy number variations in a test sample includes, for example, the following steps: a computing device 110 determines the breakpoint location of the copy number variation in the test sample via a prediction model; calculates the ratio of the sequencing depth of each bed region of the test sample to the mean sequencing depth of the baseline bed regions of a constructed control set; compares the calculated ratio with a predetermined ratio threshold to determine the ploidy of each bed region; counts the proportion of each ploidy bed region within the breakpoint region; determines whether the counted proportion of each ploidy bed region is greater than or equal to a predetermined proportion threshold; and, in response to determining that the proportion of the bed region with the current ploidy is greater than or equal to the predetermined proportion threshold, determines the ploidy of the stably amplified targeted analysis region as the current ploidy.

[0062] It should be understood that traditional methods for determining ploidy directly based on the mean copy number are prone to errors due to copy number errors, potentially leading to incorrect ploidy determinations in certain regions, such as classifying heterozygous individuals as diploid. The method of this invention corrects the ploidy of bed regions by using the ratio of the mean sequencing depth of each bed region in the test sample to the mean of the baseline bed regions in the constructed control set. This results in more accurate ploidy data for the entire region, avoiding the impact of local copy number errors on regional detection results. For example, if 6 out of 10 bed regions are classified as diploid, the traditional method would classify the entire region as diploid; however, if it is actually heterozygous (monozygous), this invention can correct it to an accurate detection result, such as "heterozygous," while the traditional method would result in a false negative. Therefore, this invention improves the overall accuracy of copy number variation detection results.

[0063] Specifically, in this invention, the method for determining the ploidy of each bed region includes, for example: if the ratio of the average sequencing depth of the current bed region to the baseline bed region of the constructed control set is between 0.2 and 0.7 (i.e., greater than or equal to 0.2 and less than 0.7), then the ploidy of the current bed region is heterozygous deletion; if the ratio is less than 0.2, then the ploidy of the current bed region is homozygous deletion; if the ratio is between 7 and 1.3 (i.e., greater than or equal to 0.7 and less than 1.3), then the ploidy of the current bed region is diploid; and if the ratio is between 1.3 and 1.7 (i.e., greater than or equal to 1.3 and less than 1.7), then the ploidy of the current bed region is triploid.

[0064] Regarding the predetermined percentage threshold, it is, for example, but not limited to, 50% of the number of beds in the breakpoint region. For example, if the percentage of a bed region of a certain ploidy is higher than 50% of the number of beds in the breakpoint region, then that region is of that ploidy. For example, if the percentage of a triploid bed region is higher than 50% of the number of beds in the breakpoint region, then the ploidy of the target region is triploid.

[0065] Table 2 below exemplifies the detection results using a conventional method for detecting copy number variations based on amplicon sequencing data. Table 3 below exemplifies the detection results using the method for detecting copy number variations based on amplicon sequencing data of the present invention.

[0066] Table 2

[0067]

[0068] Table 3

[0069]

[0070]

[0071] In Tables 2 and 3 above, the column information indicated by the item "Thalassemia deletion type identified by qPCR method" indicates the actual thalassemia deletion type identified for each sample number. The column information corresponding to the item "Copy Number" indicates the "heterozygous / homozygous status". "1" represents heterozygous, and "0" represents homozygous. Taking the detection result of "Copy Number Variation Type" of blood sample CN990 as an example, according to Table 2, the "Copy Number Variation Type" detected by the traditional method is "SEA heterozygous". This detection result is inconsistent with the actual identified thalassemia deletion type "SEA complex heterozygous deletion 4.2". Table 2 shows that the actual identified thalassemia deletion types show that the copy number variation type of blood sample CN990 is a complex type, with both homozygous and heterozygous components. In the traditional detection method, the detection data for the breakpoint region is accurate, but homozygous deletion regions are not detected. According to Table 3, the "copy number variation type" detected by the method of the present invention is "α4.2 / SEA complex". This detection result is consistent with the actual identified thalassemia deletion type "SEA complex heterozygous 4.2 deletion", accurately detecting homozygous deletion regions and heterozygous deletion regions.

[0072] In the above scheme, by using filtered alignment results data, the overall uniformity of the amplicon sequencing regions of each sample is statistically analyzed to determine the stable amplified target analysis regions; based on the sequencing depth of the stable amplified target analysis regions, a control set baseline is constructed; and based on the breakpoint location of copy number variants in the test sample determined by the prediction model and the constructed control set baseline, the detection results of copy number variants in the test sample are generated. This invention can obtain accurate copy number variant detection results based on amplicon sequencing data. Experimental data show that the accuracy of this invention is improved by 50% compared to traditional methods (e.g., directly using CNVkit software amplicon mode for copy number variant detection).

[0073] Figure 3 A flowchart of a method 300 for determining a stable amplified targeted analysis region according to an embodiment of the present invention is shown. It should be understood that method 300 can, for example, be implemented in... Figure 6 The described electronic device is executed at point 600. It can also be used in... Figure 1 The described computing device 110 performs the operation. It should be understood that method 300 may also include additional actions not shown and / or the actions shown may be omitted, and the scope of the invention is not limited in this respect.

[0074] At step 302, the computing device 110 uses the start and end positions of the sequence alignment to the reference genome as a bed region and counts the sequence coverage within each bed region. For example, the sequence coverage of a bed region is N.

[0075] At step 304, the computing device 110 corrects the sequence coverage number for each bed region using the number of aligned sequences of the test sample to obtain the corrected sequence coverage number for each bed region. For example, the computing device 110 corrects the sequence coverage number N of the bed region using the number of aligned sequences of each sample to remove the influence of data volume differences on the interval sequence coverage number.

[0076] The following example, using formula (3), illustrates the algorithm for calculating the corrected sequence coverage within each bed region.

[0077] Rn=N / M (3)

[0078] In formula (3) above, Rn represents the number of corrected sequence coverages in each bed region. N represents the number of sequence coverages in each bed region. M represents the number of aligned sequences.

[0079] In step 306, the computing device 110 extracts regions from the negative control set where all samples have sequence coverage and the difference in the number of corrected sequence coverages is less than a predetermined correction threshold, based on the alignment result data filtered for uniformity, to identify them as targeted analysis regions for stable amplification.

[0080] Regarding the predetermined correction threshold, it is, for example, but not limited to, 40%.

[0081] For example, computing device 110 performs statistical analysis on each sample in the negative control sample set (which is, for example, formed by summarizing multiple batches of data to reflect the diversity of the experiment as much as possible), screening for bed regions with read coverage and a difference in corrected Rn values ​​(CV) of <40%, in order to identify the screened bed regions as stable amplification target analysis regions. For example, the screened bed regions are labeled as CNV target beds for subsequent copy number variation analysis based on CNV target beds. In some embodiments, if the above requirements are not met for important regions of copy number variation, it is recommended to optimize the stability of local regions of the amplicon kit.

[0082] It should be understood that CNV detection itself involves detecting the copy number difference between the test sample and the control set. In traditional detection methods, the instability of the amplicon itself can easily lead to large fluctuations in the copy number difference between the test sample and the control set. These fluctuations stem from both systemic and probe design issues. By employing the aforementioned methods, the region screened by this invention becomes stable in negative samples.

[0083] Table 4 below exemplarily illustrates regions in the 213 samples extracted by the present invention that are covered and have a difference of less than 40% in the corrected sequence coverage, serving as targeted analysis regions for stable amplification.

[0084] Table 4

[0085]

[0086]

[0087]

[0088]

[0089]

[0090] In the above scheme, the present invention can significantly reduce the impact of the instability of the amplicon itself on the copy number variation detection results by determining the target analysis region of stable amplification, thereby making the copy number variation detection results more accurate.

[0091] Figure 4 A flowchart of a method 400 for sequencing depth correction targeting a determined, stably amplified region, according to an embodiment of the present invention, is shown. It should be understood that method 400 can, for example, be used in... Figure 6 The described electronic device is executed at point 600. It can also be used in... Figure 1 The described computing device 110 performs the operation. It should be understood that method 400 may also include additional actions not shown and / or the actions shown may be omitted, and the scope of the invention is not limited in this respect.

[0092] At step 402, computing device 110 uses the alignment result data filtered for uniformity to statistically analyze the sequence coverage depth of each bed region in the stable amplified targeted analysis region of the test sample and the control set sample.

[0093] The alignment result data filtered for uniformity is obtained, for example, through the following steps: the computing device 110 statistically analyzes the T 20%X coverage and Fold 80 value of the sample to be tested in order to obtain the overall uniformity of the amplicon sequencing region of the sample to be tested; and based on the comparison result of the obtained overall uniformity with a predetermined uniformity threshold, the remaining alignment result data is filtered to obtain the alignment result data filtered for uniformity.

[0094] At step 404, computing device 110 performs correction based on the number of aligned sequences, the number of bed regions, and the average sequencing depth of the test sample and the control set samples.

[0095] The following example, using formulas (4) and (5), illustrates the algorithm for calculating the corrected average sequencing depth within each bed region.

[0096] Ri=Ni / (Q / n) (4)

[0097] Mi=Ri / (median(Ri...Rn)) (5)

[0098] In formulas (4) and (5) above, i represents the bed region number. Q represents the number of aligned sequences in the sample to be tested. n represents the number of bed regions (i.e., the total number of bed regions). Mi represents the average sequencing depth of the i-th bed region after correction. Ri represents the average sequencing depth of the i-th bed region. Ni represents the sequence coverage of the i-th bed region.

[0099] At step 406, computing device 110 uses a local multinomial regression equation to correct GC preference and amplicon length preference to obtain the standardized depth of the targeted analysis region of the stable amplification after correction. In some embodiments, the multinomial regression equation is, for example, locally weighted regression (LOESS) based on the R language, i.e., R LOESS.

[0100] Figure 5 A flowchart illustrating a method 500 for detecting copy number variations based on amplicon sequencing data according to an embodiment of the present invention is shown. It should be understood that method 500 can, for example, be implemented in... Figure 6 The described electronic device is executed at point 600. It can also be used in... Figure 1 The described computing device 110 performs the operation. It should be understood that method 500 may also include additional actions not shown and / or the actions shown may be omitted, and the scope of the invention is not limited in this respect.

[0101] In step 502, DNA is extracted from different sample types of the sample to be tested to construct sequencing libraries using an amplicon sequencing kit and sequencing platform in order to generate sequencing data for the sample to be tested.

[0102] Methods for generating sequencing data for test samples include, for example: extracting and quantifying sample DNA; amplifying the target region using multiplex PCR primers (i.e., first round of PCR); purifying the amplification product using 0.9 times the volume of magnetic beads according to experimental requirements (i.e., magnetic bead purification); adding sequencing adapters to the purified product according to the requirements of the Illumina library preparation kit (i.e., second round of PCR); quantifying the library concentration using a Qubit or similar instrument, and determining the library fragment length using an Agilent DNA 1000 kit or other equivalent reagents; and performing high-throughput sequencing using the Illumina Nextseq CN500 next-generation sequencing platform to obtain the raw sequencing sequence, i.e., amplicon-based sequencing data.

[0103] In step 504, quality control is performed on the generated sequencing data to generate preprocessed sequencing data.

[0104] Methods for quality control of generated sequencing data include, for example, filtering out sequence data with multiple alignments and sequence data with alignment quality values ​​below a predetermined quality value. For example, for preprocessed sequencing data, preliminary quality control is performed using data quality control software (such as, but not limited to, FastP, FastQC, Trim-Adaptor, etc.). The purpose is to filter out sequences containing N, low-quality sequences, and sequences containing adapters from the raw sequencing data to obtain high-quality sequencing sequences for subsequent analysis. For example, using FastP software to remove sequencing adapters, low-quality bases (base quality <10), and filters sequences shorter than 40 bp from the raw sequencing data, thereby obtaining high-quality usable data.

[0105] In step 506, the preprocessed sequencing data is compared with the sequencing data of a reference genome to obtain alignment result data. This alignment result data at least indicates the positional information of each sequence on the reference genome.

[0106] For example, using the MEM algorithm of BWA software, the high-quality usable data generated in step 504 is compared and analyzed with the HG38 standard reference genome sequence to obtain the alignment position information of each sequence on the reference genome sequence and obtain the original alignment information file.

[0107] In step 508, sequencing data that meets predetermined filtering conditions is filtered out from the alignment result data.

[0108] For example, the samtools software can be used to filter out unaligned sequences, sequences aligned at multiple positions, and sequences with an alignment quality value lower than 20, in order to obtain the alignment results data that are filtered out.

[0109] In step 510, the overall uniformity of the amplicon sequencing region is obtained based on the filtered alignment result data, so as to filter the filtered alignment result data based on the obtained overall uniformity to obtain alignment result data filtered by uniformity.

[0110] For example, using the bamdst software to perform in-depth statistics on BAM files, based on the raw alignment information file obtained above, the alignment rate, coverage, and T 20% x coverage rate (%) of the samples are statistically analyzed, along with the Fold 80 value for each sample, to obtain the overall uniformity of the amplicon sequencing region. Then, based on the comparison between the obtained overall uniformity and a predetermined uniformity threshold, the remaining alignment results are filtered to obtain uniformity-filtered alignment results. For example, if the Fold 80 value of the current sample is >3.3 (e.g., 3.3 is one of several predetermined uniformity thresholds corresponding to the Fold 80 value), the current sample is filtered out.

[0111] In step 512, the target analysis region for stable amplification is determined.

[0112] It should be understood that in actual experiments, due to the influence of multiple factors, some primers may amplify unstablely or there may be cross-amplification between primers in multiplex PCR. In order to ensure the accuracy of CNV detection, it is necessary to determine the region of stable amplification.

[0113] Methods for determining stable amplification target regions include, for example, using the start and end positions of sequence alignment to the reference genome as a bed region, counting the sequence coverage within each bed region (equivalent to one read alignment region being one bed region), correcting the sequence coverage within each bed region using the number of aligned sequences in the test sample to obtain a corrected sequence coverage within each bed region (equivalent to correction using the number of reads in a fully aligned sequence); based on the alignment results data filtered for uniformity, extracting regions where all samples in the negative control set have sequence coverage and where the difference in corrected sequence coverage is less than a predetermined correction threshold, thus identifying these as stable amplification target analysis regions. By employing the above methods, the differences in data volume between individual samples can be eliminated, ensuring overall data consistency. It should be understood that this invention, through a target region screening strategy based on negative control samples, provides a stable foundation for downstream CNV detection.

[0114] In step 514, sequencing depth correction is performed on the identified stable amplified target analysis region.

[0115] For example, using the filtered alignment information file, the sequence coverage depth of each bed region within the stably amplified target analysis region (CNV target bed region) of the test sample and control set samples is statistically analyzed; then, the average depth of all samples is corrected using the method 400 mentioned above to obtain the corrected sequencing depth of the stably amplified target analysis region.

[0116] At step 516, based on the sequencing depth of the targeted analysis region via corrected stable amplification, the mean and standard deviation of each bed region in the control set samples are calculated to construct the control set baseline.

[0117] For example, based on the depth of the corrected target region obtained in step 512, a control set training set is established using the reference module of cnvkit. In some embodiments, control sets are established separately for each sex, targeting sex chromosomes.

[0118] At step 518, based on the prediction model used to predict the breakpoint location of copy number variation in the test sample and the constructed control set baseline, the detection result of copy number variation in the test sample is generated.

[0119] For example, the ratio (i.e., the mean depth of each bed region in the test sample to the baseline bed region in the control set) is calculated. If the ratio is between 0.2 and 0.7, the corresponding bed region is heterozygous deletion; if the ratio is < 0.2, the corresponding bed region is homozygous deletion; if the ratio is between 0.7 and 1.3, the corresponding bed region is diploid; and if the ratio is between 1.3 and 1.7, the corresponding bed region is triploid. Then, the proportion of each ploidy bed region within the breakpoint region is statistically analyzed. If the proportion of a certain ploidy bed region is higher than 50% of the total number of beds in the breakpoint region, then that ploidy is considered the ploidy of the entire region.

[0120] In the above scheme, the present invention can systematically provide a copy number variation detection and analysis method based on targeted amplicon sequencing data, which can significantly improve the accuracy of copy number variation detection results, taking into account the characteristics of multiplex amplicon sequencing kits.

[0121] In some embodiments, the present invention also provides an apparatus for copy number variation detection and analysis based on targeted amplicon sequencing data. This apparatus includes, for example, a data quality control module, a targeted region screening module, a targeted region depth homogenization module, and a CNV detection module.

[0122] In terms of quality control of sequencing data after sequencing, this invention filters out sequencing data that meets predetermined filtering conditions from the alignment results data, and obtains the overall uniformity of the amplicon sequencing region based on the remaining alignment results data after filtering. This allows for further filtering of the remaining alignment results data based on the obtained overall uniformity to obtain alignment results data filtered by uniformity. Thus, this invention adds a means of evaluating coverage uniformity in terms of quality control of sequencing data after sequencing.

[0123] In the targeted region screening module, this invention uses the start and end positions of sequence alignment to the reference genome as a bed region. For the sequence coverage of each bed region, the number of aligned sequences of the test sample is used for correction to obtain the corrected sequence coverage of each bed region. Based on the alignment results data filtered for uniformity, regions in the negative control set where all samples have sequence coverage and the difference in corrected sequence coverage is less than a predetermined correction threshold are extracted to determine the stable amplification targeted analysis regions. Therefore, this invention proposes a targeted region screening strategy based on negative control samples in the targeted region screening module, providing a stable foundation for downstream CNV detection.

[0124] Regarding the target region depth homogenization module, this invention uses alignment result data filtered for homogenization to statistically analyze the sequence coverage depth of each bed region in the stable amplified target analysis region of the test sample and the control set samples; it corrects based on the number of aligned sequences, the number of bed regions, and the average sequencing depth of the test sample and the control set samples; and it uses a local multinomial regression equation to correct GC bias and amplicon length bias to obtain the standardized depth of the stable amplified target analysis region after correction.

[0125] In CNV module detection, this invention employs a prediction model to predict the breakpoint location of copy number variants in the test sample; calculates the ratio of the average sequencing depth of each bed region in the test sample to the mean of the baseline bed regions in the constructed control set; compares the calculated ratio with a predetermined ratio threshold to determine the ploidy of each bed region; statistically analyzes the proportion of bed regions of each ploidy within the breakpoint region; and determines the ploidy of the target region based on the comparison between the statistically analyzed proportion of bed regions of each ploidy and a predetermined proportion threshold. Therefore, this invention provides a correction and optimization method for copy number variant detection, significantly improving overall detection accuracy.

[0126] Figure 6 A block diagram schematically illustrates an electronic device 600 suitable for implementing embodiments of the present invention. The electronic device 600 may be used to implement... Figures 2 to 5 The methods shown range from 200 to 500. (For example...) Figure 6As shown, the electronic device 600 includes a central processing unit (i.e., CPU 601), which can perform various appropriate actions and processes according to computer program instructions stored in read-only memory (i.e., ROM 602) or loaded from storage unit 608 into random access memory (i.e., RAM 603). The RAM 603 may also store various programs and data required for the operation of the electronic device 600. The CPU 601, ROM 602, and RAM 603 are interconnected via bus 604. An input / output interface (i.e., I / O interface 605) is also connected to bus 604.

[0127] Multiple components in electronic device 600 are connected to I / O interface 605, including: input unit 606, output unit 607, and storage unit 608. CPU 601 executes the various methods and processes described above, such as executing methods 200 to 500. For example, in some embodiments, methods 200 to 500 may be implemented as computer software programs stored in a machine-readable medium, such as storage unit 608. In some embodiments, part or all of the computer program may be loaded and / or installed on electronic device 600 via ROM 602 and / or communication unit 609. When the computer program is loaded into RAM 603 and executed by CPU 601, one or more operations of methods 200 to 500 described above may be performed. Alternatively, in other embodiments, CPU 601 may be configured to execute one or more actions of methods 200 to 500 by any other suitable means (e.g., by means of firmware).

[0128] It should be further noted that the present invention can be a method, apparatus, system, and / or computer program product. A computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the present invention.

[0129] Computer-readable storage media can be tangible devices capable of holding and storing instructions for use by an instruction execution device. Computer-readable storage media can be, for example, but not limited to, electrical storage devices, magnetic storage devices, optical storage devices, electromagnetic storage devices, semiconductor storage devices, or any suitable combination thereof. More specific examples (a non-exhaustive list) of computer-readable storage media include: portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), static random access memory (SRAM), portable compact disc read-only memory (CD-ROM), digital multifunction disc (DVD), memory sticks, floppy disks, mechanical encoding devices, such as punch cards or recessed protrusions storing instructions thereon, and any suitable combination thereof. The computer-readable storage media used herein are not to be construed as transient signals themselves, such as radio waves or other freely propagating electromagnetic waves, electromagnetic waves propagating through waveguides or other transmission media (e.g., light pulses through fiber optic cables), or electrical signals transmitted through wires.

[0130] The computer-readable program instructions described herein can be downloaded from computer-readable storage media to various computing / processing devices, or downloaded via a network, such as the Internet, local area network, wide area network, and / or wireless network, to an external computer or external storage device. The network may include copper transmission cables, fiber optic transmission, wireless transmission, routers, firewalls, switches, gateway computers, and / or edge servers. A network adapter card or network interface in each computing / processing device receives the computer-readable program instructions from the network and forwards them to the computer-readable storage media in the respective computing / processing device.

[0131] The computer program instructions used to perform the operations of this invention may be assembly instructions, instruction set architecture (ISA) instructions, machine instructions, machine-dependent instructions, microcode, firmware instructions, state setting data, or source code or object code written in any combination of one or more programming languages, including object-oriented programming languages ​​such as Smalltalk, C++, etc., and conventional procedural programming languages ​​such as the "C" language or similar programming languages. The computer-readable program instructions may be executed entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving a remote computer, the remote computer may be connected to the user's computer via any type of network—including a local area network (LAN) or a wide area network (WAN)—or may be connected to an external computer (e.g., via the Internet using an Internet service provider). In some embodiments, electronic circuitry, such as programmable logic circuitry, field-programmable gate arrays (FPGAs), or programmable logic arrays (PLAs), is personalized by utilizing state information from the computer-readable program instructions. This electronic circuitry can execute the computer-readable program instructions to implement various aspects of the invention.

[0132] Various aspects of the present invention are described herein with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It should be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer-readable program instructions.

[0133] These computer-readable program instructions can be provided to a processor in a voice interaction device, a general-purpose computer, a special-purpose computer, or a processing unit of another programmable data processing device, thereby producing a machine such that, when executed by the processing unit of the computer or other programmable data processing device, these instructions create means for implementing the functions / actions specified in one or more blocks of the flowchart and / or block diagram. These computer-readable program instructions can also be stored in a computer-readable storage medium that causes a computer, programmable data processing device, and / or other device to operate in a particular manner; thus, the computer-readable medium storing the instructions comprises an article of manufacture that includes instructions for implementing aspects of the functions / actions specified in one or more blocks of the flowchart and / or block diagram.

[0134] Computer-readable program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be performed on the computer, other programmable data processing apparatus, or other device to produce a computer-implemented process, thereby causing the instructions executed on the computer, other programmable data processing apparatus, or other device to perform the functions / actions specified in one or more boxes of a flowchart and / or block diagram.

[0135] The flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of devices, methods, and computer program products according to various embodiments of the present invention. In this regard, each block in a flowchart or block diagram may represent a module, segment, or portion of an instruction containing one or more executable instructions for implementing a specified logical function. In some alternative implementations, the functions marked in the blocks may occur in a different order than those shown in the drawings. For example, two consecutive blocks may actually be executed substantially in parallel, and they may sometimes be executed in reverse order, depending on the functions involved. It should also be noted that each block in the block diagrams and / or flowcharts, and combinations of blocks in the block diagrams and / or flowcharts, may be implemented using a dedicated hardware-based system that performs the specified function or action, or using a combination of dedicated hardware and computer instructions.

[0136] The various embodiments of the present invention have been described above. These descriptions are exemplary and not exhaustive, nor are they limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments. The terminology used herein is chosen to best explain the principles, practical application, or technical improvements to the embodiments in the market, or to enable others skilled in the art to understand the embodiments disclosed herein.

[0137] The above are merely optional embodiments of the present invention and are not intended to limit the present invention. For those skilled in the art, the present invention can have various modifications and variations. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method of detecting copy number variations based on amplicon sequencing data, characterized in that, The method comprises: aligning sequencing data obtained via a pre-processing with sequencing data of a reference genome, so as to obtain alignment result data, the sequencing data being obtained via an amplicon sequencing technology; filtering out, for the alignment result data, sequencing data satisfying a predetermined filtering condition; obtaining overall uniformity of an amplicon sequencing region based on the alignment result data left after filtering, so as to determine a stably amplified targeted analysis region; constructing a control set baseline based on sequencing depth of the stably amplified targeted analysis region; generating a detection result about copy number variation of the to-be-tested sample based on a breakpoint position of the copy number variation of the to-be-tested sample determined by a prediction model and the constructed control set baseline, wherein generating the detection result about copy number variation of the to-be-tested sample comprises: determining, via the prediction model, the breakpoint position of the copy number variation of the to-be-tested sample; calculating a ratio of sequencing depth of each bed region of the to-be-tested sample to a mean of sequencing depth of bed regions of the constructed control set baseline; and comparing the calculated ratio with a predetermined ratio threshold, so as to determine a ploidy of each bed region, wherein determining the stably amplified targeted analysis region comprises: filtering, based on the obtained overall uniformity, the alignment result data left after filtering, so as to obtain alignment result data filtered via uniformity; and obtaining the alignment result data filtered via uniformity comprises: counting T 20% X coverage and Fold 80 value of the to-be-tested sample, so as to obtain overall uniformity of the amplicon sequencing region of the to-be-tested sample; and filtering, based on a comparison result of the obtained overall uniformity and a predetermined uniformity threshold, the alignment result data left after filtering, so as to obtain alignment result data filtered via uniformity, determining the stably amplified targeted analysis region further comprises: counting a sequence coverage in each bed region, taking a start position and an end position of sequence alignment to the reference genome as one bed region; correcting, for the sequence coverage of each bed region, the number of aligned sequences of the to-be-tested sample, so as to obtain a corrected sequence coverage in each bed region; and extracting, based on the alignment result data filtered via uniformity, a region in which all samples in the negative control set are covered by sequences and a difference in corrected sequence coverage is less than a predetermined correction threshold, so as to determine the region as the stably amplified targeted analysis region.

2. The method of claim 1, wherein, Generating the detection result about copy number variation of the to-be-tested sample further comprises: counting a proportion of each-ploidy bed region in the breakpoint region; determining whether the counted proportion of each-ploidy bed region is greater than or equal to a predetermined proportion threshold; and in response to determining that the proportion of the current-ploidy bed region is greater than or equal to the predetermined proportion threshold, determining that the ploidy of the stably amplified targeted analysis region is the current ploidy.

3. The method of claim 1, wherein, The sequencing data satisfying the predetermined filtering condition comprises: sequence data of multi-position alignment and sequence data with an alignment quality value lower than a predetermined quality threshold.

4. The method of claim 1, wherein, Constructing the control set baseline based on the sequencing depth of the stably amplified targeted analysis region comprises: performing sequencing depth correction for the determined stably amplified targeted analysis region; and Based on the sequencing depth via the corrected targeted regions, the mean and standard deviation of each bed region of the control set samples are calculated in order to construct the control set baseline.

5. The method of claim 4, wherein, The sequencing depth correction for the determined stably amplified targeted analysis regions comprises: Using the aligned data via uniformity filtering, the sequence coverage depth of each bed region of the stably amplified targeted analysis regions of the test sample and the control set samples are counted; Based on the number of aligned sequences, the number of bed regions, the average sequencing depth of the test sample and the control set samples are corrected; and Using a local multinomial regression equation to correct GC bias, amplicon length bias to obtain the normalized depth of the stably amplified targeted analysis regions after correction.

6. A computing device, comprising: Comprise: at least one processing unit; at least one memory coupled to the at least one processing unit and storing instructions for execution by the at least one processing unit, the instructions, when executed by the at least one processing unit, cause the apparatus to perform steps of the method according to any one of claims 1 to 5.

7. A computer-readable storage medium, characterized in that, A computer program is stored on a computer readable storage medium, the computer program, when executed by a machine, implements the method according to any one of claims 1 to 5.

Images