A mutation noise filtering calculation method based on second-generation sequencing data

Abstract

The application discloses a mutation background noise filtering calculation method based on second-generation sequencing data, and relates to the technical field of high-throughput sequencing data analysis, to solve the problems that the prior art depends on a baseline, has poor adaptability and high false positives. The algorithm first preprocesses sequencing data to obtain a pileup file, classifies sites according to sequence context and calculates error rates, constructs a hierarchical background noise model to test mutations, combines sample types to qualitatively determine real mutations and background noise, and simultaneously evaluates sample quality. The algorithm does not need to construct a baseline in advance, is suitable for hybrid capture and amplicon sequencing, can improve mutation detection specificity and consistency of repeated samples, guarantees reliable results, and is suitable for tumor gene detection and other scenes.

Description

Technical Field

[0001] This invention relates to the field of high-throughput sequencing data analysis technology, and in particular to a mutation background noise filtering calculation method based on second-generation sequencing data. Background Technology

[0002] Existing sequencing data background noise filtering techniques mainly fall into two categories: one is through experimental workflow optimization, such as using molecular tags to analyze reads with the same molecular tag and remove some noise during read deduplication. However, molecular tags have extremely high requirements for sequencing depth, resulting in high experimental and sequencing costs. The other approach involves estimating the noise level during the post-deduplication analysis and using models to filter single-base mutations. For model-based filtering methods, fitting the background noise pattern is a challenge. Since different sequencers have different background noise patterns, sample quality and experimental procedures also affect the background noise pattern.

[0003] One approach to model filtering methods involves constructing a baseline. This involves pre-estimating the background noise error rate (PRR) from a batch of baseline samples and using this estimated PRR to filter single-base mutations in the sequencing samples. This method requires a certain number of baseline samples and is costly. Whenever the sequencer or experimental workflow is changed, baseline samples need to be collected and reconstructed; otherwise, the PRR may not accurately reflect the actual situation of the sequencing samples, leading to overly lenient or overly strict mutation filtering, thus affecting the accuracy of the detection results. Even within the same sample, the PRR may differ depending on the sequence context and depth.

[0004] Other methods do not rely on a baseline but directly build models within the sample, learning noise patterns under different sequence contexts and filtering sites based on these learned noise patterns. For this approach, currently published noise estimation and modeling methods cannot comprehensively handle these differences, resulting in poor performance of filtering algorithms at some sites. Furthermore, due to the characteristics of amplicon sequencing, the amplification efficiency of different PCR templates varies, causing some reads containing erroneous signals to be amplified more easily than others. This phenomenon makes the background noise variation at different sites in amplicon sequencing much greater than in capture sequencing, rendering some model filtering methods unusable in amplicon sequencing environment. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of existing technologies and provide a single-base mutation (SNV) background noise filtering model algorithm that does not require pre-construction of a baseline, is unrestricted by sample type, and can be applied to hybridization capture sequencing or amplicon sequencing data. While maintaining a certain level of sensitivity, it improves the specificity of SNV reporting results and the consistency among duplicate samples. Furthermore, this algorithm can also be extended to complex mutation types.

[0006] To achieve the above objectives, this invention provides a mutation back noise filtering calculation method based on second-generation sequencing data, the steps of which are as follows:

[0007] S1. Preprocess the sequencing data of the sequencing sample and obtain the pileup file of the sample at each capture interval or amplicon interval site;

[0008] S2. Based on the base type of the current site, the sequence context information of the mutation location, the site depth, and the sequencing direction, the sites are classified, and the global background noise error rate and the error rate of each type of site at a specific depth are calculated.

[0009] S3. Construct a first-layer background noise model for each error type of each base, estimate the first-layer error rate, and perform a first-layer test on sample mutations based on the first-layer background noise model.

[0010] S4. Count the number of potential noise sites that pass the first layer test in each error type. If the number of potential noise sites is greater than or equal to the preset threshold v, then construct a second layer background noise model for the error type and perform a second layer test on the matched mutations based on the second layer background noise model.

[0011] S5. Based on the test results, qualitative analysis is performed on single-base mutations and multi-base mutations to distinguish between real mutations and background noise.

[0012] Preferably, in S1, the specific process of the preprocessing is as follows:

[0013] S1.1 Align the sequenced bases to the human genome reference sequence to obtain the original BAM file;

[0014] S1.2 If the sequencing sample is a capture sequencing sample, deduplication is performed on the original BAM file; if the sequencing sample is an amplicon sequencing sample, deduplication is not performed.

[0015] S1.3 Perform a local re-comparison of the BAM file to obtain the processed BAM file;

[0016] S1.4 Using the processed BAM file, count the number of positive strand reads with base quality > BQ and sequencing quality > MQ at each base site within the capture interval, and the number of negative strand reads (D).+ and D - ), and the positive and negative support numbers (SRs) of various base types at that position that satisfy base quality > BQ and sequencing quality > MQ. + and SR - The output is a pileup file.

[0017] For samples prepared using strand-specific library construction, a positive strand read (+) indicates that the original DNA molecule originated from the positive strand, and a negative strand read (-) indicates that the original DNA molecule originated from the negative strand.

[0018] For samples constructed using a non-strand-specific library method, a positive read (+) means that read1 is aligned to the positive strand of the reference genome sequence, and a negative read (-) means that read1 is aligned to the complementary strand of the reference genome sequence.

[0019] Preferably, S2 specifically includes:

[0020] S2.1, Base Classification:

[0021] The base types at the current site include four types: A, T, C, and G.

[0022] The sequence context information includes bases at n upstream positions and m downstream positions. For each of the n upstream bases, there are four possibilities: A, T, C, and G, for a total of 4... n Types; for the m downstream bases, each position has four possibilities: A, T, C, and G, for a total of 4. m Types;

[0023] Arrange and combine the base types of the current site, upstream of the current site, and downstream of the current site, for a total of 4. n+m+1 There are three base types, where n and m are preset positive integers;

[0024] S2.2 The screening rule for potential noise sites is: merge positive and negative chain reads, and calculate the total number of supporting reads SR for each point in the current classification. total =SR + +SR - mutation frequency Freq total =(SR + +SR - ) / (D + +D - ) and SR total If a certain mutation occurs at a site, the Freq... total ≥f and SR totalIf the value is ≥s, the mutation is considered to be highly likely to be a true mutation; otherwise, it is considered a potential noise site. Subsequently, when calculating the error rate and estimating the error rate formula parameters, these sites are removed, and only potential noise sites are retained.

[0025] S2.3, regarding the 4 n+m+1 For each base in the base type, calculate the global background noise error rate for different error types. The global background noise error rate is the ratio of the total number of erroneous supports to the total depth; specifically:

[0026] i. For each base context sequence, the error rate for three error types can be calculated. For sites with base A, the error rates are A→T, A→C, and A→G; for sites with base C, the error rates are C→A, C→T, and C→G; for sites with base G, the error rates are G→A, G→T, and G→C; and for sites with base T, the error rates are T→A, T→C, and T→G.

[0027] ii. If forward and reverse reads are not distinguished: Reverse reads are categorized according to the base type after reverse complementation, and the upstream and downstream base types, and are counted together with the corresponding forward reads after reverse complementation. Extract all sites in the corresponding context and calculate the total depth ∑D of these sites. context And the total number of supported reads ∑SR for each error type. context,subtype The global error rate g for each error type in this context. context,subtype =∑SR context,subtype / ∑D context ;

[0028] iii. If distinguishing between forward and reverse reads: for each error type within each context, the global error rates for forward and reverse reads are calculated separately, and are respectively g context,subtype,+ and g context,subtype,- ;

[0029] S2.4, regarding the 4 n+m+1 For each base in the base type, calculate the error rate of different error types at a specific depth d. The error rate at a specific depth is the ratio of the total number of error supports to the total depth within a specified depth interval, specifically:

[0030] i. For each base context sequence, the error rate for three error types can be calculated;

[0031] ii. Select sites in the base context with depths in the range of (dw / 2, d+w / 2) to calculate the error rate at depth d, where w is the pre-defined width of the depth interval;

[0032] iii. If forward and reverse reads are not distinguished: Reverse reads are categorized according to the base type after reverse complementation, and the upstream and downstream base types, and are counted together with the corresponding forward reads after reverse complementation. Calculate the total depth ∑D of these sites. context,(d-w / 2,d+w / 2] Calculate the total number of supported reads ∑SR for each error type. context,subtype,(d-w / 2,d+w / 2] At a specific depth d, the error rate e for each error type in this context. context,subtype,d =∑SR context,subtype,(d-w / 2,d+w / 2] / ∑D context,(d-w / 2,d+w / 2] ;

[0033] iv. If distinguishing between forward and reverse reads: for each error type within each context, the error rates at a specific depth d for forward and reverse reads are calculated separately, and are respectively e context,subtype,d,+ and e context,subtype,d,- .

[0034] Preferably, S3 specifically includes selecting one from a depth-independent constant model, a depth-dependent power function model, a depth-dependent exponential model, or a threshold model as the first-layer background noise model, and estimating the first-layer error rate formula e. T1,context,subtype (d): Estimate the parameters of each model formula by least squares method or maximum likelihood method, and select the optimal formula and parameters from them;

[0035] The first-level error rate is calculated for all mutations in the sample using the selected optimal formula and parameters. The first-level test was then performed. This first-level test employed a binomial test, with the null hypothesis that the mutation support number follows a binomial distribution (Binomial(D, e)). The p-values ​​for the positive and negative strands were calculated. Specifically, this included:

[0036] a) Select formula {A, B, C, ...} as candidates from constant models that are independent of depth, power function models that are dependent on depth, exponential models that are dependent on depth, or threshold models. Estimate the parameters contained in each candidate formula and select the optimal formula and parameters from them.

[0037] b) The data points used for parameter estimation can be any of the following:

[0038] i. The depth of potential noise sites {d1, d2, d3, ..., d} after filtering by the first layer of the background noise model. x Supported read counts {SR1, SR2, SR3, ..., SR} x} and error rates {e1, e2, e3, ..., e x};

[0039] ii. The specific depth {d1, d2, d3, ..., d...} calculated according to the method described in S2.4. x}, the total number of supported reads at the corresponding depth {SR1, SR2, SR3, ..., SR x} and error rates {e1, e2, e3, ..., e x};

[0040] c) The estimation method can be chosen from the least squares method and the maximum likelihood method:

[0041] i. Least squares method:

[0042] 1) Using different parameters for the candidate formulas, calculate the error rate for each data point estimated according to the formula. , , , ..., };

[0043] 2) Calculate the error rate estimate { , , , ..., } and the actually observed {e1, e2, e3, ...,e x The sum of squared residuals between} ;

[0044] 3) When using latent noise sites for parameter estimation, the sites with the largest q% of residual squares can be removed each time the RSS is calculated to avoid the influence of outliers on the estimation results;

[0045] 4) For each candidate formula, select the parameter that minimizes the RSS.

[0046] 5) When trying multiple formula forms simultaneously, perform parameter estimation for each tried formula {A, B, C, ...}. After completing the parameter estimation for each formula, calculate the sum of squared residuals {RSS} between the error rate estimates and the actual observed values ​​for each formula. A RSS B RSS C , ...}, select the candidate formula with the smallest RSS as the final error rate estimation formula;

[0047] ii. Maximum Likelihood Method:

[0048] 1) Using different parameters for the candidate formulas, calculate the error rate for each data point estimated according to the formula. , , , ..., };

[0049] 2) According to the likelihood formula, likelihood = SR * ln( )+(D-SR)*ln(1- ), calculate the error rate for each observation point { , , , ..., At depths {D1, D2, D3, ..., D} x In}, we observe {SR1, SR2, SR3, ...,SR} x The likelihood values ​​of the support reads are {l1, l2, l3, ..., l}. x};

[0050] 3) Calculate the sum of likelihood values ​​after log processing: logsum = ∑(log(l));

[0051] 4) When using latent noise sites for parameter estimation, the sites with the smallest likelihood value (q%) can be removed each time the logsum is calculated to avoid the influence of outliers on the estimation results.

[0052] 5) For each candidate formula, select the parameter that maximizes logsum;

[0053] 6) When trying multiple formula forms simultaneously, perform parameter estimation for each tried formula {A, B, C, ...}. After completing the parameter estimation for each formula, calculate the log likelihood and {logsum} for different formulas. A , logsum B ,logsum C , ...}, select the candidate formula with the largest logsum as the final error rate estimation formula;

[0054] d) After obtaining formula e T1,context,subtype (d) After that, calculate the first-level error rate for all mutations in the sample and perform the first-level error check:

[0055] i. Extract the positive and negative read depths D for each mutation site from the pileup file. + and D - Positive and negative support for reads SR + and SR - ;

[0056] ii. For forward reads, find the first-layer error rate formula for the context and subtype to which the site belongs, and calculate the first-layer error rate e of the site at the corresponding depth. +=e T1,context,subtype (D + Perform a binomial test on the positive reads, according to D. + SR + and e + Calculate the first layer p value (p) of the forward reads. T1,+ The null hypothesis of the binomial distribution is: SR + ≤X,X~Binomial(D + , e + );

[0057] iii. For reverse reads, find the corresponding first-layer error rate formula for that type, and calculate the first-layer error rate e at the corresponding depth for the site. - =e T1,context,subtype (D - Perform a binomial test on the reverse reads, based on D. - SR - and e - Calculate the first layer p value (p) of the reverse reads. T1,- The null hypothesis of the binomial distribution is:

[0058] iv.SR - ≤X, X~Binomial(D - , e - );

[0059] Set the threshold for the first layer's p value to p. T1 If the site satisfies p T1,+ <p T1 or p T1,- <p T1 If the mutation frequency of the current mutation is significantly higher than the first layer background noise error rate in a certain direction of reads, then it is considered that the mutation frequency of the current mutation is significantly higher than the first layer background noise error rate in a certain direction of reads.

[0060] Preferably, in S3, the candidate formula for the first layer background noise model is one of the following: a constant model independent of depth, a power function model dependent on depth, an exponential model dependent on depth, or a threshold model, wherein:

[0061] The depth-independent constant model is: ; g is the value calculated in S2.3 context,subtype g context,subtype,+ or g context,subtype,- ;;

[0062] The power function model related to depth is: At this point, the error rate of the first layer increases as the depth decreases; or At this point, the error rate of the first layer decreases as the depth decreases; a and b are relevant parameters;

[0063] The depth-related exponential model is: , where the error rate of the first layer increases as the depth decreases; or , where the error rate of the first layer decreases as the depth decreases; a and b are related parameters, g is the global background noise error rate calculated in S2.3, and D is the median depth of the sample;

[0064] The form of the threshold model is related to the site depth: when the site depth ≥ maxD, ; is the global background noise error rate calculated in S2.3; when the site depth < maxD, the depth-related power function model or the depth-related exponential model is adopted. The power function model is: , where the error rate of the first layer increases as the depth decreases; or , where the error rate of the first layer decreases as the depth decreases. The form of the exponential model is: , where the error rate of the first layer increases as the depth decreases; or , where the error rate of the first layer decreases as the depth decreases, b is a related parameter, and g and maxD are constant terms.

[0065] Preferably, in S4, the second-layer background noise model is constructed based on the potential noise sites passing the first-layer test. The model form is consistent with the first-layer background noise model. The parameters are estimated by the maximum likelihood method and iteratively optimized until the number of potential noise sites passing the test is lower than the threshold v; the second-layer test uses a binomial test to calculate the p-values of the positive and negative strands, specifically including:

[0066] Merge all the base context sequences within the same error type and confirm whether it is necessary to perform second-layer filtering on a certain error type among the 12 error types. If necessary, a second-layer background noise model can be constructed for this error type, and the sites are screened to estimate the second-layer error rate formula e T2,subtype (d):

[0067] a) Statistically count the number of potential noise sites satisfying p T1,+ < p T1 or p T1,- < p T1 in each subtype. If it is greater than or equal to v, the second-layer error rate needs to be calculated for this subtype and background noise filtering is performed. v is a preset upper limit of the quantity;

[0068] b) Use the potential noise sites {pos1, pos2, pos3,..., pos zUsing data, we select formulas {A, B, C, ...} as candidates from depth-independent constant models, depth-dependent power function models, depth-dependent exponential models, or threshold models, estimate the formula parameters, and then select the optimal formula and parameters.

[0069] i. First, for {pos1, pos2, pos3, ..., pos...} z Using the maximum likelihood method, we determine the temporary formula and parameters for the second-level error rate of this subtype, and obtain the temporary formula e. t,subtype (d);

[0070] ii. According to the temporary formula e t,subtype (d) For {pos1, pos2, pos3, ..., pos} z}, calculate the temporary error rate and perform a test to obtain the temporary p-value (p) of the site. t,+ and p t,- );

[0071] iii. Set the threshold for the temporary p-value to p. t Statistical analysis of {pos1, pos2, pos3, ..., pos} y In}, p satisfies t,+ <p t or p t,- <p t Total number of sites:

[0072] iv. If the number of sites satisfying the condition is still greater than or equal to v, then repeat steps i-iii in the current step, using the currently satisfying p t,+ <p t or p t,- <p t The parameters of the site are re-estimated, and the temporary formula is updated;

[0073] v. If the number of sites meeting the condition is less than v, fix the current formula and use this temporary formula as the final second-level error rate formula e for the subtype. T2,subtype (d);

[0074] c) Obtaining formula e T2,subtype (d) After that, for all mutations in the sample that can match the second-layer error rate formula, calculate the second-layer error rate and perform second-layer error testing:

[0075] i. For each site capable of performing a second-level error check, calculate the second-level error rate and the second-level p-value (p0.05) for both forward and reverse reads. T2,+ and p T2,- );

[0076] ii. Set the threshold for the second layer p value to p. T2 If the site satisfies p T2,+ <p T2 or p T2,- <p T2 If the mutation frequency in a certain direction of reads is significantly higher than the second-layer background noise error rate, then it is considered that the mutation frequency of the current mutation is significantly higher than the second-layer background noise error rate.

[0077] Preferably, in S5, based on the calculated p T1,+ p T1,- p T2,+ and p T2,- Each single-base mutation is then given a final qualitative analysis:

[0078] For amplicon sequencing samples, a mutation is considered a true mutation if it meets the following conditions; otherwise, it is considered background noise:

[0079] If the error type of the mutation does not constitute a second-layer model, the following requirements must be met: the p-value of the positive strand must be less than the threshold, or the p-value of the negative strand must be less than the threshold, i.e., p T1,+ <p T1 , or p T1,- <p T1 ;

[0080] If the error type of the mutation constructs a second-layer model, it requires that: the positive strand p-value < threshold and the positive strand second-layer p-value < threshold, or the negative strand p-value < threshold and the negative strand second-layer p-value < threshold, that is: p T1,+ <p T1 And p T2,+ <p T2 , or p T1,- <p T1 And p T2,- <p T2 ;

[0081] For hybridization capture sequencing samples, a mutation is considered a true mutation if it meets the following conditions; otherwise, it is considered background noise:

[0082] If the error type of the mutation does not construct a second-layer model, the following requirements must be met: the p-value of the positive strand must be less than the threshold and the p-value of the negative strand must be less than the threshold, i.e., p T1,+ <p T1 , and p T1,- <p T1 ;

[0083] If the error type of the mutation constructs a second-layer model, it requires that: the positive strand p-value < threshold and the positive strand second-layer p-value < threshold, and the negative strand p-value < threshold and the negative strand second-layer p-value < threshold, i.e., p T1,+ <p T1 And p T2,+ <p T2, and p T1,- <p T1 And p T2,- <p T2 .

[0084] Preferably, in S5, the rule for qualitative analysis of multi-base mutations is as follows: the multi-base mutation is broken down into multiple single-base mutations; if all single-base mutations are determined to be background noise, then the multi-base mutation is background noise; otherwise, it is a true mutation, specifically including:

[0085] b) For any multi-base substitution mutation in the sample (the number of substituted bases is S), decompose the mutation into S single-base mutations;

[0086] c) After disassembly, for each mutation s in the S single-base mutations:

[0087] i. Find the corresponding first-level error rate formula, perform a binomial distribution test, and obtain the test values ​​(p) for the forward and backward reads. s,T1,+ and p s,T1,- );

[0088] ii. For s, if the corresponding subtype has a second-layer model, find the corresponding second-layer error rate formula and obtain the p-values ​​of the positive and negative reads (p... s,T2,+ and p s,T2,- )

[0089] iii. Compare the p-value with the threshold to determine whether s belongs to background noise;

[0090] If all single-base mutations identified from the S single-base mutations are determined to be background noise, then the multi-base substitution type mutation is background noise; otherwise, the mutation is determined to be a real mutation.

[0091] Preferably, S5 also includes a sample quality assessment step, which can estimate the overall noise level of the sample using formulas in the model and the median depth of the sample, calculate corresponding indicators, and provide feedback on the sample quality: calculate the overall background noise error rate and detection limit for each error type; overall background noise error rate = max(average first layer error rate, second layer error rate); detection limit = max(minimum supported reads in the first layer, minimum supported reads in the second layer) / median depth of the sample; combine complementary error types for evaluation, and if the detection limit exceeds a preset threshold, the sample quality is determined to be abnormal. Specifically, this includes:

[0092] a) Calculate the overall noise error rate for each subtype:

[0093] i. Calculate the median depth D within the sample amplification interval. If the sample is a hybridization capture sample, calculate the deduplicated D;

[0094] ii. Based on the first-layer model, the error rate e at depth D can be calculated for different contexts and subtypes. T1,context,subtype (D). According to the 12 error types (A>T, A>C, A>G, C>A, C>T, C>G, G>A, G>T, G>C, T>A, T>C, T>G), classify e... T1,context,subtype (D) Merge the results and calculate the geometric mean of the Level 1 error rates of all contexts corresponding to each subtype, which will be used as the average Level 1 error rate e for this subtype. T1,subtype ;

[0095] iii. If the subtype has a second-layer model, calculate the error rate e of this subtype at depth D based on the second-layer model obtained in S4. T2,subtype ;

[0096] iv. For each subtype, calculate the overall noise error rate, e subtype =max(e T1,subtype , e T2,subtype );

[0097] b) Calculate the detection limit (LOD) for each subtype of the sample. subtype :

[0098] i. Based on the subtype e T1,subtype After deduplication of the samples, find the median depth D that allows the null hypothesis to be SR. T1,subtype ≤X, X~Binomial(D, p=e T1,subtype The test value p) T1,subtype <p T1 Minimum SR T1,subtype value sigSR T1,subtype ;

[0099] ii. Similarly, if the subtype has a second-level model, according to e T2,subtype And D, find the second-level test value p that makes the subtype satisfy. T2,subtype <p T2 sigSR T2,subtype ;

[0100] iii.LOD subtype =max(sigSR T1,subtype ,sigSR T2,subtype ) / D;

[0101] c) Since each mutation type has both forward and reverse support reads, reads in either direction can pass the background noise test to detect the mutation. The LOD of a subtype is actually determined by the LODs of that subtype and its complementary subtype. Therefore, for each complementary error type (subtype-pair) subtype A and subtype B (AB-pair), they share the error rate and detection limit, and their overall noise error rate and detection limit (LOD) can be combined:

[0102] i. The 12 error types can be combined into 6 pairs according to complementary types (e.g., C>T and G>A).

[0103] ii. Overall error rate of AB-pair: e AB-pair =min(e A , e B );

[0104] iii. Overall LOD of AB-pair: LOD AB-pair =min(LOD A LOD B );

[0105] d) Based on the LOD of the subtype or subtype-pair, if the following conditions are met, the current sample is considered to have a quality problem and cannot obtain the expected limit of detection for all mutation types:

[0106] i. For amplicon sequencing samples designed with single-end primers, if a certain subtype has a LOD... subtype The detection limit is higher than the preset detection limit L;

[0107] ii. For amplicon sequencing samples designed with paired primers, if a certain subtype-pair has a LOD... subtype-pair The detection limit is higher than the preset detection limit L;

[0108] e) For each complementary AB-pair, calculate the difference between their totalities and the difference between their LODs:

[0109] i. Overall error rate difference Δe between AB pairs AB-pair =abs(e A , e B );

[0110] ii. LOD difference ΔLOD of AB-pair AB-pair =abs(LOD A LOD B );

[0111] f) Based on the magnitude of the difference, determine whether there is an imbalance in the subtype error rate among the samples. If the Δe of a certain subtype-pair... subtype-pair Errors exceeding the preset error rate difference threshold Δe, or ΔLOD subtype-pair If the error rate difference exceeds the preset LOD difference threshold ΔLOD, it is considered that the sample has an imbalance in error rate among complementary subtypes, and the sample may have quality problems.

[0112] According to specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0113] It eliminates the need to build a baseline in advance, avoiding the high costs and adaptation challenges of baseline-dependent technologies, and lowers the application threshold.

[0114] Wildcard hybridization capture sequencing and amplicon sequencing are not limited by sequencing platforms or sample quality, and have wide applicability;

[0115] The hierarchical background noise model accurately handles noise differences in different sequence contexts and sequencing depths, significantly improving the specificity of single / multi-base mutation detection and reducing false positives;

[0116] It improves the consistency of mutation reporting among duplicate samples, and can calculate the overall noise error rate and detection limit of the sample, evaluate sample quality in real time, and ensure reliable results;

[0117] It expands its capabilities to include multi-base substitution mutation filtering, making it applicable to core areas such as tumor gene detection.

[0118] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0119] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0120] Figure 1 This is a graph showing the number of SNVs reported by hybridization capture sequencing plasma samples before and after filtering, based on an embodiment of a mutation background noise filtering calculation method for second-generation sequencing data according to the present invention.

[0121] Figure 2 This is a graph showing the overall error rate of each error type in hybridization capture sequencing plasma samples at different median depths, according to an embodiment of the present invention.

[0122] Figure 3This is a LOD diagram of different error types in hybridization capture sequencing plasma samples at different median depths, according to an embodiment of the present invention. Detailed Implementation

[0123] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0124] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0125] Example 1

[0126] This embodiment uses plasma and blood cell paired sample hybridization capture sequencing data as the experimental object to verify the performance of the first layer background noise model and the sample quality assessment effect described in S1-S3 and S5 of this patent.

[0127] A total of 976 pairs of clinical plasma and blood cell matching samples were used for validation. Library construction experiments were performed on these 976 pairs of clinical samples, resulting in 976 blood cell libraries and 976 plasma libraries.

[0128] S1. Obtain sequencing data from hybridization capture samples and perform data preprocessing.

[0129] Hybridization capture was performed on 976 library pairs. The captured libraries were then sequenced at both ends using a T7 sequencing platform to obtain raw bcl files. All bcl files were split using bcl2fastq to obtain sequencing fastq data for each sample.

[0130] For each sample's FastQ data, the adapter sequences of each read were removed using Trimmomatic, and the remaining sequences of the read were aligned to the human reference genome hg19 using BWA to obtain the original BAM file. After deduplication, these reads were re-aligned using GATK Realigner to obtain more accurate alignment positions, and stored in the preprocessed BAM file.

[0131] S2, Mutation Detection and Annotation

[0132] Using Vardict software and Mutect, paired mutation detection was performed on the BAM files of each pair of preprocessed plasma and blood cell samples. After obtaining the mutation list, ANNOVAR was used to annotate the mutations detected by the two software programs. The annotated mutations were then combined and stored in a file for later filtering.

[0133] S3. Construct a first-layer background noise model within each sample and test for sample mutations based on the first-layer background noise model.

[0134] Following the method described in step S1.4, set BQ=25 and MQ=0 to obtain the pileup file for each plasma sample. Following the method described in step S2.2, remove sites that may contain actual mutations. When filtering sites, set f=0.03 and s=6 to obtain the filtered pileup site set.

[0135] Set the parameters n=1 and m=1 in step S2.1. Based on the base type of each site and the bases upstream and downstream of the site, the filtered sites are divided into 64 base classes. Each base class contains 3 error types, distinguishing between positive and negative strands, resulting in a total of 384 combinations of bases and error types. Following the method described in step S2.3, calculate the global background noise error rate g for each base + error type combination on both the positive and negative strands. context,subtype,+ and g context,subtype,- ;

[0136] Following the method described in step S2.4, calculate the background noise error rate e on the positive and negative strands for each combination of base + error types at depths of 50x, 100x, 150x, 200x, ..., 1900x, 1950x, and 2000x. context,subtype,d,+ and e context,subtype,d,- ;

[0137] The candidate formulas for all classifications are set to be either depth-independent constant models or depth-dependent power function models whose error rate decreases with depth. Using e... context,subtype,d,+ and e context,subtype,d,- As training data, the optimal first-layer back noise rate formula e for each of the 384 types of sites was obtained using the least squares method. T1,context,subtype,strand (d);

[0138] Set the first-level test threshold p T1 =0.001, for each single-base mutation after annotation, find the corresponding first-layer background noise model, and calculate the first-layer test value (p) for positive and negative reads at each site. T1,+ and p T1,- Following the rules in S5, each mutation is qualitatively analyzed to determine whether the site is a true mutation.

[0139] S4. Filter the mutations detected by Vardict and Mutect and statistically compare the number of reported sites before and after filtering.

[0140] The number of SNVs in each clinical plasma sample before and after filtration was recorded. The specific changes in the reported SNV count for each sample are shown below. Figure 1 As shown, after filtering, the overall number of SNVs decreased, and the difference in the number of reported SNVs between samples was reduced, indicating that the present invention can reduce the number of noise sites reported in abnormal samples.

[0141] S5. Estimate the overall noise level of the sample and verify the sample quality.

[0142] For each sample, calculate the median depth, and according to the methods described in steps a) and b) of the sample quality assessment in S5, calculate the overall noise error rate and the median depth for each error type of the sample, as well as p. T1 LOD under the condition of 0.001. Compare the relationship between the overall noise error rate and the median depth of the samples in different samples, such as... Figure 2 As shown, the overall error rate of the sample shows a slight negative correlation with the median depth. Comparing the relationship between LOD and median depth across different samples, as shown... Figure 3 As shown, there is a significant negative correlation between median depth and LOD: the lower the median depth, the higher the LOD. When L=0.5% in step S5 d), 7 out of 976 plasma samples had LOD>L for certain subtype pairs and were judged to have quality problems. The median depth of these 7 samples after deduplication was all below 550x.

[0143] Example 2

[0144] The experimental subjects in this embodiment are amplicon sequencing data of cell lines and tissue samples, used to verify the effectiveness of the first and second layer background noise models described in S1-S5 of this invention in amplicon sequencing.

[0145] The sample information is as follows:

[0146] a) Cell line simulation validation samples

[0147] Using cell line samples, a positive cell line sample with mutations was mixed with a negative cell line sample to obtain a cell line simulation validation sample. Eleven known single-base mutations within the amplicon amplification region were identified, as shown in Table 1. Three replicate library construction experiments were performed on the cell line simulation validation sample, yielding three libraries: L1, L2, and L3.

[0148] Table 1. Mutation information of amplicon sequencing cell line simulation validation samples

[0149]

[0150] b) Clinical validation samples

[0151] A total of 60 clinical tumor tissue samples were collected. Library construction experiments were performed on these 60 clinical samples, resulting in 60 libraries.

[0152] S1. Obtain amplicon sequencing data and perform data preprocessing.

[0153] Amplicons designed with paired-end primers were used to enrich libraries from three cell line mock validation samples and 60 clinical sample libraries. For the three cell line mock validation libraries, single-end sequencing was performed using the MGI-200 sequencing platform to obtain raw bcl files. For the 60 clinical sample libraries, each library was sequenced once on seven sequencing platforms: Miseq, Miniseq, FASTASeq300, GenoLab M, DNBSEQ-E25, and DNBSEQ-G99, to obtain raw bcl files. Seven sequencing data points were available for each clinical sample, which could be used to compare their consistency.

[0154] All BCL files were split using bcl2fastq to obtain sequencing FastQ data for each sample. For each sample's FastQ data, Trimmomatic was used to remove adapter sequences from each read, and the remaining sequences were aligned to the human reference genome hg19 using bwa to obtain the raw BAM file. These reads were then re-aligned using GATK Realigner to obtain more accurate alignment positions, and stored in the preprocessed BAM file.

[0155] S2, Mutation Detection and Annotation

[0156] Using Vardict software, each preprocessed BAM file is subjected to mutation detection to obtain a mutation list. ANNOVAR is then used to annotate the detected mutations. The annotated mutations are stored in a file for later filtering.

[0157] S3. Construct the first layer of back noise model within each sample.

[0158] Following the method described in step S1.4, set BQ=25 and MQ=10 to obtain the pileup file for each sample. Following the method described in step S2.2, remove sites that may contain true mutations. When filtering sites, set f=0.1 and s=6 to obtain the filtered pileup site set.

[0159] Set the parameters n = 2 and m = 0 in S2.1. According to the base types of each site and the 2 bases upstream of the site, the filtered sites are divided into 64 categories of context. Each category of context contains 3 types of errors, resulting in a total of 192 combinations of base context sequences + error types. The types of sites on the reverse reads are classified according to the reverse complemented situation and merged with the sites on the forward reads for analysis. Calculate the global background error rate g of each combination of base context sequence + error type according to the method described in S2.3 context,subtype ;

[0160] Set the candidate formulas for all base context sequences of the A→G and T→C error types as a constant model independent of depth and a threshold model that is an exponential model when the site depth < maxD. The threshold model tries two forms where the error rate rises and falls as the depth decreases respectively; Set the candidate formulas for all base context sequences of the remaining error types as a constant model independent of depth and a threshold model that is an exponential model when the site depth < maxD. The threshold model only tries the form where the error rate rises as the depth decreases. Use the potential noise sites in the filtered pileup site set as training data, and obtain the optimal first-layer background error rate formula e T1,context,subtype (d). When calculating logsum in the maximum likelihood method each time, remove the 5% of the sites with the smallest likelihood values to avoid the influence of outliers on the estimation results;

[0161] S4. Construct the second-layer background model within each sample

[0162] Merge all contexts within the same error type, and set the first-layer test threshold p according to the method described in step S4 T1 = 0.001, and count the number of sites passing the test within each error type. Set the maximum value v of the allowed number of sites passing the test to 10, and find out the error types that need to construct the second-layer background model;

[0163] For the error types that need to construct the second-layer background model, set the candidate formula form to be consistent with the first-layer background model. Set the threshold p of the temporary p-value according to the description in S4 b) t = 0.1, use the potential noise sites passing the first-layer test, estimate the parameters by the maximum likelihood method, and obtain the optimal second-layer background error rate formula e T2,subtype (d).

[0164] S5. Filter the mutations detected by Vardict, qualitatively analyze single-base mutations and multi-base mutations, and distinguish true mutations from background noise

[0165] Set the first-layer test threshold p T1=0.001, second-level test threshold p T2 =0.001. For each single-base mutation after annotation, find the corresponding first-layer background noise model and second-layer background noise model (if any), and calculate the first-layer test value (p) for positive and negative reads at each site. T1,+ and p T1,- ) and second-level test value (p T2,+ and p T2,- Following the rules in S5, each mutation is qualitatively analyzed to determine whether the site is a true mutation.

[0166] For each polybase substitution mutation after annotation, the mutation is split according to the polybase mutation typing rules, and the test value (p) of each single base after splitting is calculated. s,T1,+ and p s,T1,- p s,T2,+ and p s,T2,- The mutations were then qualitatively analyzed to determine whether the site was a true mutation.

[0167] Then perform the following verification:

[0168] 1. Verify the filtering effect of SNVs in amplicon sequencing cell line simulation sample.

[0169] This step primarily aims to verify whether the present invention filters out genuine SNV mutations, leading to false negatives and affecting the sensitivity of mutation detection. The examination revealed that all 11 known SNV sites in the three replicate samples were detected by Vardict in all three replicates, and all 11 sites passed the two-layer model test, ultimately being retained and reported. Detailed information on the depth, mutation frequency, and p-value of these 11 mutations in the three replicates is provided. s,T1,+ p s,T1,- p s,T2,+ and p s,T2,- As shown in Table 2, some of these mutations have a mutation frequency of around 1%. Comparing the p-values ​​in Table 2 reveals that some sites only passed the test in one direction of reads, while the p-value in the other direction exceeded the threshold. This phenomenon indicates that in amplicon sequencing samples, when using bidirectional primer design, establishing models for reads in both directions separately can achieve higher sensitivity in low-frequency mutations.

[0170] Table 2. Mutation frequency, depth, and p-value of known sites in amplicon sequencing cell line simulation validation samples.

[0171]

[0172] 2. Verify the consistency between the detection of SNVs and Complex sites in clinical repeat sequencing samples.

[0173] Mutation consistency was compared between 60 clinical samples obtained from data analysis on seven different sequencing instruments. Table 3 shows the consistency statistics of pre- and post-filtered sites across different platforms. The first column in the table represents the consistency of mutation reporting across platforms (e.g., 1 / 7, meaning the mutation was reported only on one of the seven platforms). Each subsequent column shows the number of sites in the 60 clinical samples that met the reported consistency criteria, as well as the proportion of these sites among all reported sites.

[0174] Table 3. Consistency statistics of mutations reported by different sequencing platforms before and after filtering in amplicon sequencing clinical samples.

[0175]

[0176] Of the 60 samples, Vardict reported a total of 134,898 SNVs, of which only 1,589 (1.18% of the total) were reported simultaneously by 7 platforms; Vardict reported 285 Complexes, of which only 2 (0.7% of the total) were reported simultaneously by 7 platforms. After applying the first-layer model filtering of this invention, a total of 31,236 SNVs remained, of which 1,140 (3.65%) were reported simultaneously by 7 platforms; a total of 175 Complexes remained, of which 2 (1.14%) were reported simultaneously by 7 device types. After applying the two-layer model filtering, a total of 3,151 SNVs remained, of which 917 (29.1%) were reported simultaneously by 7 platforms; a total of 88 Complexes remained, of which 2 (2.27%) were reported simultaneously by 7 device types.

[0177] After applying the two-layer model filtering of this invention, a total of 3,151 SNVs remain, of which 917 (29.1%) are reported by 7 platforms simultaneously; and a total of 88 Complexes remain, of which 2 (2.27%) are reported by 7 models simultaneously.

[0178] Before filtering, sites reported by only one platform can be considered noise. Of these 98,591 SNVs and 265 Complex sites, the first-layer model of this invention can filter out 73.8% of the SNVs and 40.8% of the Complex sites, while the two-layer model can filter out 98.3% of the SNVs and 73.6% of the Complex sites.

[0179] The above results demonstrate that the filtering method described in this invention can significantly reduce the number of reported sites, and most of the reduced sites are noise sites detected only on a certain platform, thereby improving consistency. Furthermore, in amplicon sequencing samples, using a two-layer model outperforms using only a single-layer model.

[0180] The remaining technical features in the above embodiments can be flexibly selected by those skilled in the art to meet different specific practical needs according to actual circumstances. Modifications and variations made by those skilled in the art that do not depart from the spirit and scope of the present invention should be within the protection scope of the appended claims. In the above description, numerous specific details have been set forth to provide a thorough understanding of the present invention. However, it will be apparent to those skilled in the art that these specific details are not necessary to implement the present invention. In other instances, to avoid obscuring the present invention, well-known techniques, such as specific construction details, operating conditions, and other technical conditions, have not been specifically described.

[0181] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A mutation background noise filtering calculation method based on second-generation sequencing data, characterized in that, The steps are as follows: S1. Preprocess the sequencing data of the sequencing sample and obtain the pileup file of the sample at each capture interval or amplicon interval site; S2. Based on the base type of the current site, the sequence context information of the mutation location, the depth of the site, and the sequencing direction, potential noise sites are classified, and the global background noise error rate and the error rate of each type of site at a specific depth are calculated. The specific depth interval is (dw / 2, d+w / 2], where w is the pre-set depth interval width. S3. Construct a first-layer background noise model for each error type of each base, estimate the first-layer error rate, and perform a first-layer test on sample mutations based on the first-layer background noise model. S4. Count the number of potential noise sites that pass the first layer test in each error type. If the number of potential noise sites is greater than or equal to the preset threshold v, then construct a second layer background noise model for the error type and perform a second layer test on the matched mutations based on the second layer background noise model. S5. Based on the test results, qualitatively identify single-base mutations and multi-base mutations, distinguish between true mutations and background noise, and evaluate the sample quality.

2. The mutation background noise filtering calculation method based on second-generation sequencing data according to claim 1, characterized in that: In S1, the specific process of preprocessing is as follows: S1.1 Align the sequenced bases to the human genome reference sequence to obtain the original BAM file; S1.2 If the sequencing sample is a capture sequencing sample, deduplication is performed on the original BAM file; if the sequencing sample is an amplicon sequencing sample, deduplication is not performed. S1.3 Perform a local re-comparison of the BAM file to obtain the processed BAM file; S1.4 Using the processed BAM file, count the number of positive and negative reads at each base site within the capture interval where the base quality is greater than BQ and the sequencing quality is greater than MQ, as well as the number of positive and negative reads for each base type at that base site where the base quality is greater than BQ and the sequencing quality is greater than MQ. Output the results as a pileup file.

3. The mutation background noise filtering calculation method based on second-generation sequencing data according to claim 2, characterized in that: S2 specifically includes: S2.1, Base Classification: The base types at the current site include four types: A, T, C, and G. The sequence context information includes bases at n upstream positions and m downstream positions. For each of the n upstream bases, there are four possibilities: A, T, C, and G, for a total of 4... n Types; for the m downstream bases, each position has four possibilities: A, T, C, and G, for a total of 4. m Types; Arrange and combine the base types of the current site, upstream of the current site, and downstream of the current site, for a total of 4. n+m+1 There are three base types, where n and m are preset positive integers; S2.2 The screening rules for potential noise sites are as follows: merge positive and negative strand data, calculate mutation frequency and support number; if the mutation frequency is ≥ f and the support number is ≥ s, where f is the preset mutation frequency threshold and s is the preset mutation support number threshold, then it is determined to be a potential real mutation site and excluded; the remaining sites are potential noise sites. S2.3, regarding the 4 n+m+1 For each base in the base type, calculate the global background noise error rate for different error types. The global background noise error rate is the ratio of the total number of error supports to the total depth. S2.4, regarding the 4 n+m+1 For each base in the base type, calculate the error rate of different error types at a specific depth d. The error rate at a specific depth is the ratio of the total number of error supports to the total depth within a specified depth interval. The specific depth interval is (dw / 2, d+w / 2], where w is the pre-set width of the depth interval.

4. The mutation background noise filtering calculation method based on second-generation sequencing data according to claim 3, characterized in that: S3 specifically includes selecting one from depth-independent constant model, depth-dependent power function model, depth-dependent exponential model or threshold model as the first layer background noise model, estimating the parameters of each model formula by least squares method or maximum likelihood method, and selecting the optimal formula and parameters from them. The first-level error rate is calculated for all mutations in the sample using the selected optimal formula and parameters. And perform the first layer of verification; The first layer test uses a binomial test, with the null hypothesis that the number of mutation support numbers follows a binomial distribution Binomial(D, e). The p-values ​​of the positive and negative strands are calculated, where e is the first layer error rate estimated by the first layer background noise model, and D is the site depth of the corresponding sequencing strand.

5. The mutation background noise filtering calculation method based on second-generation sequencing data according to claim 4, characterized in that: In S3, the first layer background noise model is one of the following: a depth-independent constant model, a depth-dependent power function model, a depth-dependent exponential model, or a threshold model, wherein: The depth-independent constant model is: ; The global background noise error rate calculated in S2.3; The power function model related to depth is: At this point, the error rate of the first layer increases as the depth decreases; or At this point, the error rate of the first layer decreases as the depth decreases; a and b are relevant parameters; The depth-related index model is as follows: At this point, the error rate of the first layer increases as the depth decreases; or At this point, the error rate of the first layer decreases as the depth decreases; a and b are relevant parameters, g is the global back noise error rate calculated in S2.3, and D is the median depth of the sample; The form of the threshold model is related to the site depth: when the site depth ≥ maxD, ; is the global background noise error rate calculated in S2.3; when the site depth < maxD, the depth-related power function model or the depth-related exponential model is adopted; the power function model is: , at this time the first-layer error rate increases as the depth decreases, or , at this time the first-layer error rate decreases as the depth decreases; the form of the exponential model is: , at this time the first-layer error rate increases as the depth decreases, or , at this time the first-layer error rate decreases as the depth decreases, b is the relevant parameter, and g and maxD are constant terms.

6. The mutation background noise filtering calculation method based on second-generation sequencing data according to claim 5, characterized in that: In S4, the second-layer background noise model is constructed based on the potential noise sites that have passed the first-layer test. The model form is consistent with the first-layer background noise model. The parameters are estimated using the maximum likelihood method and iteratively optimized until the number of potential noise sites that have passed the test is lower than the threshold v. The second layer of testing uses a binomial test to calculate the p-values ​​for the positive and negative chains.

7. The mutation background noise filtering calculation method based on second-generation sequencing data according to claim 6, characterized in that: In S5, the rules for qualitative analysis of single-base mutations are as follows: For amplicon sequencing samples, a mutation is considered a true mutation if it meets the following conditions; otherwise, it is considered background noise: If the error type of the mutation does not build a second-layer model, the following requirements must be met: positive chain p-value < threshold or negative chain p-value < threshold; If the error type of the mutation constructs a second-layer model, the following requirements must be met: the positive chain p-value < threshold and the positive chain second-layer p-value < threshold, or the negative chain p-value < threshold and the negative chain second-layer p-value < threshold; For hybridization capture sequencing samples, a mutation is considered a true mutation if it meets the following conditions; otherwise, it is considered background noise: If the error type of the mutation does not build a second-layer model, the following requirements must be met: the p-value of the positive strand is less than the threshold and the p-value of the negative strand is less than the threshold. If the error type of the mutation constructs a second-layer model, the following requirements must be met: the positive chain p-value < threshold and the positive chain second-layer p-value < threshold, and the negative chain p-value < threshold and the negative chain second-layer p-value < threshold.

8. The mutation background noise filtering calculation method based on second-generation sequencing data according to claim 7, characterized in that: In S5, the rule for qualitative analysis of polybase mutations is as follows: break down polybase mutations into multiple single-base mutations; if all single-base mutations are determined to be background noise, then the polybase mutation is background noise; otherwise, it is a true mutation.

9. The mutation background noise filtering calculation method based on second-generation sequencing data according to claim 8, characterized in that: S5 also includes a sample quality assessment step, specifically including: Calculate the overall background noise error rate and detection limit for each error type; Overall background noise error rate = max(average first layer error rate, second layer error rate); Detection limit = max(minimum number of reads supported by the first layer, minimum number of reads supported by the second layer) / median depth of the sample; Complementary error types are evaluated together, and if the detection limit exceeds a preset threshold, the sample quality is deemed abnormal.

Images