CN111243664B - Gene variation detection method based on high-throughput sequencing - Google Patents

Abstract

The invention discloses a gene variation detection method based on high-throughput sequencing. The method comprises the following steps: preprocessing high-throughput sequencing data, mainly comprising base quality correction and read de-duplication; calculating the preference generated by sequencing, and reducing the number of effective reads supporting mutation through the preference to obtain the corrected read number supporting mutation; and calculating the correlation quality of all single samples by using the corrected read number of the supporting mutant types, then obtaining the variation quality of the single samples, and judging according to the variation quality. Experiments prove that the method improves the true positive rate and reduces the false positive rate aiming at the WGS data of the tissue sample, and possibly reduces biological tests.

Description

A Gene Variation Detection Method Based on High-throughput Sequencing

technical field

The invention belongs to the technical field of biological information, and in particular relates to a gene variation detection method based on high-throughput sequencing.

Background technique

With the development of sequencing technology, the throughput of high-throughput sequencing is increasing, so the data generated is also increasing. At the same time, the application of high-throughput sequencing is becoming more and more extensive and more and more important. Nowadays, the application of high-throughput sequencing in cancer has been very extensive and very important. For example, high-throughput sequencing is used in early tumor screening, tumor prognosis, tumor classification and tumor drug guidance. At the same time, the large amount of high-throughput sequencing data brings many computational challenges. In the high-throughput sequencing data analysis process, the most important step is variant call. Before the clinical report is issued, it is generally necessary to manually review the variant detection results reported by the variant caller.

The current methods of variant detection based on high-throughput sequencing data all have problems of missed detection and false positives in the detection results. Missed detection may lead to wrong diagnosis and cause serious harm to the patient's life (for example, reduce the probability of survival). However, a large number of false positive mutations need to increase the work of manual review, resulting in very high time and cost for manual review. All in all, the consistency between the results of automatic detection by existing software and the results of manual review is relatively poor. In addition, some variant detection software runs slowly, which consumes more computing resources and causes doctors to wait longer to get the patient's genetic test report.

Contents of the invention

In view of the above problems, the object of the present invention is to provide a gene variation detection method based on high-throughput sequencing.

In the first aspect, the present invention claims a method for detecting gene variation based on high-throughput sequencing.

The gene variation detection method based on high-throughput sequencing provided by the present invention supports single-end sequencing and paired-end sequencing data.

The gene variation detection method based on high-throughput sequencing provided by the present invention may include the following steps:

(A) Preprocessing the high-throughput sequencing data, mainly including base quality correction and read deduplication;

(B) Calculating the bias generated by sequencing, reducing the number of effective reads that support mutations through biases, and obtaining the corrected number of reads that support mutations;

(C) Using the corrected number of reads supporting the mutant type, calculate the relevant quality of all individual samples, and then obtain the variation quality of a single sample, and judge by the variation quality.

In the method, the detected genetic variation can be SNV variation or InDel variation. If there is no normal control sample available, the data corresponding to the normal control sample is considered to be a virtual file in BAM format that does not contain any reads, and this virtual file is equivalent to an empty file.

All the qualities mentioned in this invention are adopted "Ewing B, Green P. Base-calling of automated sequencer traces using phred. II. Error probabilities [J]. Genomeresearch, 1998, 8 (3): 186-194." The first The definition of quality in this formula, so any quality (for example: base quality and mutation quality) has a one-to-one correspondence with probability. So mass and probability are interchangeable.

When the detected genetic variation is an InDel variation, in step (A), the step of base quality estimation is also included before performing the base quality correction; the base quality estimation is performed as follows:

InDelQual(x, y, a, b, c) = min(x+10, y+10, STRQual(a, b))+10×log ₁₀ (c)

Among them, InDelQual is the InDel quality that has the same purpose as the base quality, x and y are the base quality of the previous base and the next base where InDel occurs, respectively, and a and b are the units of the short tandem repeat (STR) The size (repeat size) and the number of units (repeat number), c is the number of bases of InDel. The STRQual function is defined as follows:

STRQual(a,b)=44-10×log ₁₀ (1+8×softplus(a×b-8)/(θ×a ² ));

The softplus function is defined as follows:

softplus(v)=log(1+exp(v)), where log refers to the natural logarithm.

If InDel is a deletion of a STR unit, then θ=0.25, otherwise θ=1.

When the detected genetic variation is SNV or InDel variation, in step (A), base quality correction can be performed as follows:

If the high-throughput sequencing is paired-end sequencing, merge before deduplication and filtering, merge the R1 and R2 ends, and then use the base quality of the merged reads.

For the overlapping region between the R1 and R2 ends, the calculation method of the combined base quality is as follows:

When the bases corresponding to p and q are the same, calculate according to the following formula:

max(p,q);

When the bases corresponding to p and q are different, calculate according to the following formula:

max(p,q)-min(p,q);

Among them, p is the base quality before R1 end correction, and q is the base quality before R2 end correction.

For the base quality in the non-overlapping region of R1 and R2 ends, the combined base quality is equal to the base quality before correction.

If it is single-end sequencing, no merging is required (the base quality after correction is equal to the base quality before correction).

When the detected gene variation is SNV or InDel variation, in step (A), the reads can be deduplicated as follows:

Put the reads of PCR duplication or optical duplication into a family (cluster) to remove duplication.

Specifically, the reads of PCR duplication or optical duplication can be put into a family (cluster) according to the method including the following steps to remove duplication:

(a1) For each locus (every locus on each chromosome, e.g.: chr1:1, chr1:2, ..., chr1:249250621, chr2:1, ...), find in this gene The number of reads starting or ending at the seat, referred to as the end depth.

(a2) Each locus x attracts other loci y, and calculates the attraction strength s of x to y. The attraction strength s is proportional to the depth of the end of x (the number of reads starting/terminating at this locus) (proportionality constant =1), which is inversely proportional to the depth of the end of y (proportionality constant=1) and has an exponential decay relationship with the number of bases z between x and y (attenuation coefficient=5);

The formula for the strength of attraction s is as follows;

s=x/y×5 ^-z ;

If s > 1 (y is strongly attracted to x), then y is considered to be a PCR or optical repeat of x, thus putting x and y into a family.

(a3) Filtering: calculate the threshold value (threshold) of base quality, the formula is as follows:

Threshold(P)=argmax _p∈P (p ^{-|{q∈P|q≤p}|} );

Among them, P is the set of sequencing error probabilities corresponding to all base qualities (one base quality corresponds to one sequencing error probability); the mathematical symbol ∈ represents the set belonging, so p∈P means that p is any element in P , and q∈P means that q is any element in P.

Filter the bases on the reads according to the base quality threshold: filter out the bases whose corrected base quality is lower than the base quality threshold; if a base is filtered out, then it is considered that the base is not Existence, it can be understood that no base has been detected.

It should be noted that filtering bases and filtering reads are two concepts.

After filtering, for each family, take the most common base of the family (the number of times the base appears the most) as the consensus family base after deduplication, that is, the base after deduplication; here the family can be either for the base or Can be for reads.

Then, integrate the repeated reads into a consensus read (consensus read); here it needs to be mentioned that no matter whether the deduplication is reads or bases, the final result is the same after deduplication. Therefore, it can be understood as deduplication of reads, or as filtering of bases.

In step (B), the four preferences of deduplication bias, position bias, strand bias and mismatch bias are calculated for each chromosome chain; then the highest bias is selected as the overall preference of the chromosome chain; the number of reads supporting the mutant type is corrected , the number of reads after correction is equal to the number of reads before correction divided by the overall preference.

When the detected gene variation is SNV or InDel variation, the four preferences of deduplication bias, position bias, strand bias and mismatch bias can be calculated according to the following formula:

Bias(a ₁ , a ₂ , b ₁ , b ₂ )=max(0, OR(LS(a ₁ , a ₂ , b ₁ , b ₂ ))-1)*(b ₂ ÷(b ₂ +b ₁ ))+1;

In the formula, OR stands for odds ratio; OR is calculated as follows:

OR(a ₁ , a ₂ , b ₁ , b ₂ )=(a ₁ ×b ₂ )÷(a ₂ ×b ₁ ).

In the formula, LS represents Laplace smoothing; LS formula is as follows: LS (a, b, c, d)=Laplace ₂ (Laplace ₁ (a, b, c, d)); Wherein, Laplace ₁ formula is as follows: Laplace ₁ ( a, b, c, d)=(a+p, b+p, c+p, d+p); Laplace ₂ formula is as follows: Laplace ₂ (a, b, c, d)=(a+p×α , b+p×α, c+p×β, d+p×β);

In the Laplace ₁ formula and Laplace ₂ formula, the default value of p is 0.5;

In the Laplace ₂ formula, the formulas of α and β are as follows: α=max(0, p×((a ₁ +a ₂ )/(b ₁ +b ₂ )-1)); β=max(0, p×( b ₁ +b ₂ )/(a ₁ +a ₂ )-1);

For deduplication bias, the meaning of each parameter is as follows:

a ₁ : the number of reads that support any alleles removed by deduplication;

a ₂ : After deduplication, the number of reads supporting any allele;

b ₁ : the number of reads supporting mutant alleles removed by deduplication;

b ₂ : After deduplication, the number of reads supporting the mutant allele.

If there are 7 reads before deduplication and 2 reads are left after deduplication, then 7-2=5 reads are removed after deduplication.

For position bias, the meaning of each parameter is as follows:

For each direction (a total of left and right directions), extend 1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66 bases according to the direction; then substitute the following four values:

a ₁ : Among the reads that have spanned the extended locus, the number of reads supporting any allele;

a ₂ : Among the reads that do not span the extended locus, the number of reads supporting any allele;

b ₁ : among the reads that have spanned the extended locus, the number of reads supporting the mutant allele;

b ₂ : Number of reads supporting the mutant allele among reads that do not span the extended locus.

Substituting 11 extended units in turn produces a sequence consisting of 11 values; taking every 3 adjacent values in the sequence as an array, thereby generating 9 arrays; taking the minimum value for each array, thereby generating 9 minimum values; Take the maximum value of the 9 minimum values as position bias.

For strand bias, substitute the following values:

a ₁ : In the forward direction of the chromosome strand S, the number of reads supporting any allele;

a ₂ : In the reverse direction of the chromosome strand S, the number of reads supporting any allele;

b ₁ : in the forward direction of the chromosome strand S, the number of reads supporting the mutant allele;

b ₂ : In the reverse direction of chromosome strand S, the number of reads supporting the mutant allele.

Simultaneously calculate base quality imbalance (abbreviated as BQI, represented by I in the formula, with 1 I) and strandposition imbalance (SPI for short, represented by I in the formula, with 2 Is), and calculate the values of 3 Is by the following method, Then take the minimum value; then divide the strand bias by this minimum value to normalize the strand bias; finally ensure that the strand bias is at least 1;

For BQI, c is the average base quality of all alleles of another chain, d is the average base quality of the mutant alleles of this chain, and I=10 ^dc ;

每一条read有两个末端：左边和右边。因此左末端和右末端分别产生两个I。For SPI, c is the average base distance to the end of the read for all alleles of the other strand, d is the average base distance to the end of the read for the mutant alleles of the suspected biased strand, and

Each read has two ends: left and right. So two I's are generated at the left and right ends respectively.

For the mismatch bias, let x=1, 2, 3, ..., 11; x is the threshold of the number of mismatches (can be inferred by the number of mismatches < x). 1, 2, ..., 11 are thresholds for different attempts.

For each value of x, substitute the following four values:

a ₁ : the number of mismatches < x, the number of reads supporting any allele;

a ₂ : The number of mismatches ≥ x, the number of reads supporting any allele;

b ₁ : the number of mismatches < x, the number of reads supporting the mutant allele;

b ₂ : the number of mismatches ≥ x, the number of reads supporting the mutant allele;

Substituting 11 x values in turn produces a sequence consisting of 11 values; taking every 5 adjacent values in the sequence as an array, thus generating 7 arrays; taking the minimum value for each array, thereby generating 7 minimum values; Take the maximum value of the 7 minimum values as the mismatch bias.

When the detected gene variation is SNV or InDel variation, the step (C) can be carried out as follows:

(C1) For each site that needs to calculate the quality, use the following formula to calculate family quality, single-strand quality, double-strand quality;

Quality (p, B, b, r) = log _r (BIAS (B, p × B, B, b) × (r-1) + 1) × 10 × log ₁₀ (OR (B, p × B, B , b)).

in:

(1) For the reads in a family, let r=2 and use the following definition to calculate the family quality:

p: error probability corresponding to the threshold value of base quality calculated in step (a3) reads deduplication;

B: The number of reads that support any allele from the actual observation of the data;

b: Number of reads actually observed from the data supporting the mutant allele.

(2) If the high-throughput sequencing experiment uses the PCR amplification method to build the library, then r=1.5, and if the high-throughput sequencing experiment uses the capture or WGS method to build the library, then r=1.25. The families compared to a chromosome chain are all on a single-strand; for a single-strand, use the following definition to calculate the single-strand quality (the maximum family quality is as follows: for C: G > T: A mutant form is 44, For T: A>C:G mutant form is 48, for other SNV mutant forms is 52, for all InDel mutant forms is 60). If the single-strand quality is greater than the maximum family quality value, then let the single-strand quality be equal to the maximum family quality value;

p: The error probability corresponding to each possible family quality (ie try all family qualities);

B: The number of families actually observed to support any allele;

b: The number of families actually observed to support the mutant allele;

After trying all p, take the largest single-strand quality as the final single-strand quality.

(3) The forward or reverse direction of each single-strand on a double-strand; for a single-strand, use the following definition to calculate the double-strand quality:

v ₁ '=v ₁ +v ₂ ×min(1, (min(w ₁ +w ₂ ,w ₀ )-w ₁ )/w ₁ );

Among them, v ₁ ′ is the final mutation quality on a certain chain, v ₁ is the mutation quality of this chain (that is, the final single-strand quality), v ₂ is the mutation quality on another chain, and w ₁ is the mutation quality of this chain. The average value of family quality, w ₂ is the average value of the family quality of another chain, w ₀ is equal to the maximum family quality (the maximum family quality is as follows: for C: G>T: A mutant form is 44, for T: A>C: G mutant form is 48, for other SNV mutant forms is 52, for all InDel mutant forms is 60) plus 10;

Exchange v ₁ and v ₂ and exchange w ₁ and w ₂ to calculate v ₂ ′ according to the symmetry of the formula for calculating v ₁ ′;

v ₂ ′=v ₂ +v ₁ ×min(1, (min(w ₁ +w ₂ , w ₀ )−w ₂ )/w ₂ ).

(C2) Define double-strand quality as the quality of variation of a single sample.

(C3) According to the above (C1) and (C2), calculate the single-sample mutation quality of the tumor sample and the normal control sample, respectively. Then calculate TLOD. TLOD represents the probability that the variation is neither derived from the body (soma) nor from the embryo (germline), as follows:

First, calculate the output value of the following function.

TNreward(a ₁ , a ₂ , b ₁ , b ₂ )=maxH(TNrewardBinomial(a ₁ , a ₂ , b ₁ , b ₂ ), PLQTN(b ₂ /b ₁ , a ₂ /a ₁ ));

Among them, TNrewardBinomial(a ₁ , a ₂ , b ₁ , b ₂ )=10/log(10)×(a1×KL Bernoulli(b ₂ /b ₁ , a ₂ /a ₁ ))

If max(a,-a)<max(b,-b), then maxH(a,b)=a, otherwise maxH(a,b)=b;

PLQTN(t, n)=3×10×log ₁₀ (min(t/n, 10 ^2.5/3 )).

TNpenal(a ₁ , a ₂ , b ₁ , b ₂ )=max(0, min(DSVQnormal, PLQ(a ₂ /a ₁ ))-12.5×(max(0, OR(a ₁ , a ₂ , b ₁ , b ₂ )-1)) ² ); where DSVQnormal refers to the single-sample mutation quality of the normal sample, that is, the double-strand quality of the normal sample.

a ₁ , a ₂ , b ₁ and b ₂ have the following meanings:

a ₁ : The sum of the family quality of the reads supporting any allele in the paired normal control sample;

a ₂ : The sum of the family quality of the reads supporting the mutant allele in the paired normal control sample;

b ₁ : the sum of family quality of reads supporting any allele in the tumor sample;

b ₂ : The sum of the family quality of the reads supporting the mutant allele in the tumor sample.

Then, calculate TLOD, the calculation formula is as follows:

TLOD(a ₁ ,a ₂ ,b ₁ ,b ₂ )=min(DSVQtumor,PLQ(b ₂ /b ₁ ))+TNreward(a ₁ ,a ₂ ,b ₁ ,b ₂ )-TNpenal(a ₁ ,a ₂ , b ₁ , b ₂ ).

Wherein, PLQ(f)=90+3×10×log ₁₀ (f). DSVQtumor refers to the single-sample mutation quality of the tumor sample, that is, the double-strand quality of the tumor sample.

(C4) Calculate NLOD. NLOD represents the probability that the variation originates from the embryo (germline). details as follows:

First, for each sample (there are two samples of tumor and paired normal in total), the genotype likelihood (GL) of each genotype (GT) is calculated. The definitions and descriptions of GT and GL are in https://samtools.github.io/hts-specs/VCFv4.2.pdf. To sum up, homref (homozygous wild type), hetero (heterozygous type) and homalt (homozygous mutant type) are three types of GT. The calculation methods of the GL of these three GTs are as follows:

GL(f, homref) = max(GLPowerLaw(f, homref), GLBinomial(f, homref));

GL(f, hetero) = max(GLPowerLaw(f, hetero), GLBinomial(f, hetero));

GL(f, homalt) = max(GLPowerLaw(f, homalt), GLBinomial(f, homalt)).

Among them, GLPowerLaw(f, homref) = 10×3×min(0, log _1o (0.02×(1-f)/f));

GLPowerLaw(f, hetero) = 10×3×min(0, log ₁₀ (min(f/(1-f), (1-f)/f)));

GLPowerLaw(f, homalt) = 10×3×min(0, log ₁₀ (0.02×f/(1-f)));

GLBinomial(f, homref) = 10/log(10) × d × KL Bernoulli(min(0.02, f), f);

GLBinomial(f, hetero)=10/log(10)×d×KL Bernoulli(0.5,f);

GLBinomial(f,homalt)=10/log(10)×d×KL Bernoulli(max(1-0.02,f),f);

KL Bernoulli(x,y)=y×log(y/x)+(1-y)×log((1-y)/(1-x));

If x=0, y=0, x is an infinitive or y is an infinitive, then KL Bernoulli(x,y)=0. For example, (0/0) is an infinitive.

And a refers to the number of reads supporting a mutant gene, d refers to the number of reads supporting any genotype, and f=a/d is the mutation frequency.

Then calculate the NLOD. NLOD=max(Qtumorhomref, Qnormalhomref).

where on the tumor sample, Qtumorhomref=-3-max(GL(f, homalt), GL(f, hetero)). Qtumorhomref is calculated using tumor sample data, so it has nothing to do with paired normal samples.

Wherein on the paired normal samples, Qnormalhomref=GL(f, homref)-max(GL(f, homalt), GL(f, hetero)). Qnormalhomref is calculated using paired normal sample data, so it has nothing to do with tumor samples.

(C5) The final mutation quality is: min(TLOD, NLOD+31). If the final mutation quality is at least 60, it is judged as positive and the mutant form is output; otherwise, it is judged as negative. A quality of 60 represents a false positive rate of one in a million. The mutation detection software MuTect (Cibulskis K, Lawrence M S, Carter S L, et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples[J]. Naturebiotechnology, 2013, 31(3): 213.) indirectly used The threshold of 60, and the maximum comparison quality that BWA can output is 60. Therefore a threshold of 60 was used. The mutation output format is VCF. The VCF format definition is described in detail in (Danecek P, Auton A, Abecasis G, et al. The variant call format and VCFtools[J]. Bioinformatics, 2011, 27(15): 2156-2158.).

In the second aspect, the present invention claims a gene variation detection system based on high-throughput sequencing.

The gene variation detection system based on high-throughput sequencing as claimed in the present invention includes device A, device B and device C;

The device A can preprocess the high-throughput sequencing data according to step (A) in the first aspect above, mainly including base quality correction and read deduplication;

The device B can calculate the preference generated by sequencing according to the step (B) in the first aspect above, and reduce the number of effective reads supporting the mutation through the preference to obtain the corrected number of reads supporting the mutation;

The device C can calculate the relevant quality of all individual samples according to the step (C) in the first aspect above by using the corrected number of reads supporting the mutant type, and then obtain the variation quality of a single sample, and judge by the variation quality.

In the third aspect, the present invention claims any of the following applications:

(1) Application of the method described in the first aspect above or the system described in the second aspect above in the preparation of products for early tumor screening, tumor prognosis, tumor classification and/or tumor drug guidance;

The product can detect gene variation from high-throughput sequencing data according to the steps of the method described in the first aspect above.

(II) Application of the method described in the first aspect above or the system described in the second aspect above in early tumor screening, tumor prognosis, tumor classification and/or tumor drug guidance.

Experiments prove that, for the WGS data of tissue samples, the present invention increases the true positive rate and reduces the false positive rate, and may reduce biological tests.

The present invention has the following advantages due to the adoption of the above technical scheme:

1. Almost no missed detection, that is to say, almost no false negative test results.

2. The number of false positive mutations produced is very small.

3. The software implemented using the algorithm described in the present invention runs faster than most software.

Description of drawings

Figure 1 is a performance comparison between the method of the present invention and similar software on data with molecular labels.

Figure 2 is a performance comparison between the method of the present invention and similar software on data with molecular labels.

Figure 3 is a performance comparison between the method of the present invention and similar software on data with molecular labels.

Fig. 4 is a performance comparison between the method of the present invention and similar software on data with molecular labels.

Fig. 5 is the performance comparison of the method of the present invention and similar software on tumor-normal tissue data.

Fig. 6 is the performance comparison done by the method of the present invention and similar software on tumor-normal tissue data.

In each figure, precision rate=number of true positives/(number of true positives+number of false positives). Recall = Number of True Positives/(Number of True Positives + Number of False Negatives).

Detailed ways

The experimental methods used in the following examples are conventional methods unless otherwise specified.

The materials and reagents used in the following examples can be obtained from commercial sources unless otherwise specified.

Embodiment 1, the establishment of the gene variation detection method based on high-throughput sequencing of the present invention

Taking SNV variation as an example, the high-throughput sequencing-based gene variation detection method established in the present invention is described in detail.

1. Perform high-throughput sequencing of the genomic DNA of the sample.

2. Preprocessing the data of high-throughput sequencing, mainly including: base quality correction, read deduplication; specific methods include:

The present invention adopts "Ewing B, Green P. Base-calling of automated sequencer tracing using phred.II. Error probabilities [J]. Genome research, 1998, 8 (3): 186-194." Quality ( quality), so any quality (for example: base quality and mutation quality) has a one-to-one correspondence with probability. So mass and probability are interchangeable.

The base quality correction is as follows :

a. If it is paired-end sequencing, merge before deduplication and filtering, merge the R1 and R2 ends, and then use the base quality of the merged reads.

For the overlapping region of R1 and R2, the calculation method of the combined base quality is as follows:

When the bases corresponding to p and q are the same: max(p, q);

When the bases corresponding to p and q are different: max(p, q)-min(p, q);

Among them, p is the base quality before R1 end correction, and q is the base quality before R2 end correction.

For the base quality in the non-overlapping region of R1 and R2, the combined base quality is equal to the base quality before correction.

b. If it is single-end sequencing, no merging is required (the base quality after correction is equal to the base quality before correction).

The reads are deduplicated, as follows :

The method includes: putting PCR or optically repeated reads into a family (cluster) to remove duplication, and the specific steps include the following:

a) For each locus (every locus on each chromosome, for example: chr1:1, chr1:2, ..., chrl:249250621, chr2:1, ...), find the The number of reads starting or ending (referred to as terminal depth).

b) Each locus x attracts other loci y, and calculates the attraction strength s of x to y. The attraction strength s is proportional to the depth of the end of x (the number of reads starting/terminating at this locus) (proportionality constant = 1), which is inversely proportional to the depth of the y end (proportionality constant = 1) and has an exponential decay relationship with the base number z between x and y (attenuation coefficient = 5).

If y is strongly attracted to x, then y is considered a PCR or optical duplicate of x.

The formula for the strength of attraction s is as follows:

s=x/y×5 ^-z ;

If s > 1 (y is strongly attracted to x), then y is considered to be a PCR or optical repeat of x, thus putting x and y into a family.

c) Filtering: calculate the threshold of base quality (threshold), the formula is as follows:

Threshold(P)=argmax _p∈P (p ^{-|{q∈P|q≤p}|} );

Among them, P is the set of sequencing error probabilities corresponding to all base qualities (one base quality corresponds to one sequencing error probability); the mathematical symbol ∈ represents the set belonging, so p∈P means that p is any element in P , and q∈P means that q is any element in P.

Filter the bases on the reads according to the base quality threshold: filter out the bases whose base quality after correction is lower than the base quality threshold. If a base is filtered out, it is considered that the base does not exist, and it can be understood that no base has been detected.

It should be noted that filtering bases and filtering reads are two concepts.

After filtering, for each family, the most common base (the base with the most occurrences) of the family is taken as the consensus family base after deduplication, that is, the base after deduplication. Here family can be either base-oriented or reads-oriented.

Then, integrate the repeated reads into a consensus read (consistent read). It should be mentioned here that no matter whether the deduplication is reads or bases, the final result is the same after deduplication. Therefore, it can be understood as deduplication of reads, or as filtering of bases.

3. Calculate the preference generated by sequencing, reduce the number of effective reads that support mutations through preference, and obtain the corrected number of reads that support mutations.

1. Calculate the bias generated during the sequencing process, including: deduplication bias, position bias, strand bias and mismatch bias.

Bias(a ₁ , a ₂ , b ₁ , b ₂ )=max(0, OR(LS(a ₁ , a ₂ , b ₁ , b ₂ ))-1)*(b ₂ ÷(b ₂ +b ₁ ))+1;

In the formula, OR stands for odds ratio; OR is calculated as follows:

OR(a ₁ , a ₂ , b ₁ , b ₂ )=(a ₁ ×b ₂ )÷(a ₂ ×b ₁ ).

In the formula, LS represents Laplace smoothing; LS formula is as follows: LS (a, b, c, d)=Laplace ₂ (Laplace ₁ (a, b, c, d)); Wherein, Laplace ₁ formula is as follows: Laplace ₁ ( a, b, c, d)=(a+p, b+p, c+p, d+p); Laplace ₂ formula is as follows: Laplace ₂ (a, b, c, d)=(a+p×α , b+p×α, c+p×β, d+p×β);

In the Laplace ₁ formula and Laplace ₂ formula, the default value of p is 0.5;

In the Laplace ₂ formula, the formulas of α and β are as follows: α=max(0, p×((a ₁ +a ₂ )/(b ₁ +b ₂ )-1)); β=max(0, p×( b ₁ +b ₂ )/(a ₁ +a ₂ )-1);

For deduplication bias, the meaning of each parameter is as follows:

a ₁ : the number of reads that support any alleles removed by deduplication;

a ₂ : After deduplication, the number of reads supporting any allele;

b ₁ : the number of reads supporting mutant alleles removed by deduplication;

b ₂ : After deduplication, the number of reads supporting the mutant allele.

If there are 7 reads before deduplication and 2 reads are left after deduplication, then 7-2=5 reads are removed after deduplication.

For position bias, the meaning of each parameter is as follows:

For each direction (a total of left and right directions), extend 1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66 bases according to the direction; then substitute the following four values:

a ₁ : Among the reads that have spanned the extended locus, the number of reads supporting any allele;

a ₂ : Among the reads that do not span the extended locus, the number of reads supporting any allele;

b ₁ : among the reads that have spanned the extended locus, the number of reads supporting the mutant allele;

b ₂ : Number of reads supporting the mutant allele among reads that do not span the extended locus.

Substituting 11 extended units in turn produces a sequence consisting of 11 values; taking every 3 adjacent values in the sequence as an array, thereby generating 9 arrays; taking the minimum value for each array, thereby generating 9 minimum values; Take the maximum value of the 9 minimum values as position bias.

For strand bias, substitute the following values:

a ₁ : In the forward direction of the chromosome strand S, the number of reads supporting any allele;

a ₂ : In the reverse direction of the chromosome strand S, the number of reads supporting any allele;

b ₁ : in the forward direction of the chromosome strand S, the number of reads supporting the mutant allele;

b ₂ : In the reverse direction of chromosome strand S, the number of reads supporting the mutant allele.

Simultaneously calculate base quality imbalance (abbreviated as BQI, represented by I in the formula, with 1 I) and strandposition imbalance (SPI for short, represented by I in the formula, with 2 Is), and calculate the values of 3 Is by the following method, Then take the minimum value; then divide the strand bias by this minimum value to normalize the strand bias; finally ensure that the strand bias is at least 1;

For BQI, c is the average base quality of all alleles of another chain, d is the average base quality of the mutant alleles of this chain, and I=10 ^dc ;

每一条read有两个末端：左边和右边。因此左末端和右末端分别产生两个I。For SPI, c is the average base distance to the end of the read for all alleles of the other strand, d is the average base distance to the end of the read for the mutant alleles of the suspected biased strand, and

Each read has two ends: left and right. So two I's are generated at the left and right ends respectively.

For the mismatch bias, let x=1, 2, 3, ..., 11; x is the mismatch number threshold (it can be inferred by the mismatch number<x). 1, 2, ..., 11 are thresholds for different attempts.

For each value of x, substitute the following four values:

a ₁ : the number of mismatches < x, the number of reads supporting any allele;

a ₂ : The number of mismatches ≥ x, the number of reads supporting any allele;

b ₁ : the number of mismatches < x, the number of reads supporting the mutant allele;

b ₂ : the number of mismatches ≥ x, the number of reads supporting the mutant allele;

Substituting 11 x values in turn produces a sequence consisting of 11 values; taking every 5 adjacent values in the sequence as an array, thus generating 7 arrays; taking the minimum value for each array, thereby generating 7 minimum values; Take the maximum value of the 7 minimum values as the mismatch bias.

Four preferences are calculated on each chromosome strand according to the above. Select the highest preference as the overall preference of the chromosome chain. Correct the number of reads that support the mutant type, and the number of reads after correction is equal to the number of reads before correction divided by the overall preference.

4. Use the corrected number of reads supporting the mutant type to calculate the relevant quality of all individual samples, and then obtain the variation quality of a single sample, and judge by the variation quality.

include:

1. For each locus whose quality needs to be calculated, use the following formula to calculate family quality, single-strand quality, double-strand quality, tumor-normal quality and variant quality.

Quality (p, B, b, r) = log _r (BIAS (B, p × B, B, b) × (r-1) + 1) × 10 × log ₁₀ (OR (B, p × B, B , b)).

in:

(1) For the reads in a family, let r=2 and use the following definition to calculate the family quality:

p: error probability corresponding to the threshold value of base quality calculated in step (a3) reads deduplication;

B: The number of reads that support any allele from the actual observation of the data;

b: Number of reads actually observed from the data supporting the mutant allele.

(2) If the high-throughput sequencing experiment uses the PCR amplification method to build the library, then r=1.5, and if the high-throughput sequencing experiment uses the capture or WGS method to build the library, then r=1.25. The families compared to a chromosome chain are all on a single-strand; for a single-strand, use the following definition to calculate the single-strand quality (the maximum family quality is as follows: for C: G > T: A mutant form is 44, For T: A>C:G mutant form is 48, for other SNV mutant forms is 52, for all InDel mutant forms is 60). If the single-strand quality is greater than the maximum family quality value, then let the single-strand quality be equal to the maximum family quality value;

p: The error probability corresponding to each possible family quality (ie try all family qualities);

B: The number of families actually observed to support any allele;

b: The number of families actually observed to support the mutant allele;

After trying all p, take the largest single-strand quality as the final single-strand quality.

(3) The forward or reverse direction of each single-strand on a double-strand; for a single-strand, use the following definition to calculate the double-strand quality:

v ₁ '=v ₁ +v ₂ ×min(1, (min(w ₁ +w ₂ ,w ₀ )-w ₁ )/w ₁ );

Among them, v ₁ ′ is the final mutation quality on a certain chain, v ₁ is the mutation quality of this chain (that is, the final single-strand quality), v ₂ is the mutation quality on another chain, and w ₁ is the mutation quality of this chain. The average value of family quality, w ₂ is the average value of the family quality of another chain, w ₀ is equal to the maximum family quality (the maximum family quality is as follows: for C: G>T: A mutant form is 44, for T: A>C: G mutant form is 48, for other SNV mutant forms is 52, for all InDel mutant forms is 60) plus 10;

Exchange v ₁ and v ₂ and exchange w ₁ and w ₂ to calculate v ₂ ′ according to the symmetry of the formula for calculating v ₁ ′;

v ₂ ′=v ₂ +v ₁ ×min(1, (min(w ₁ +w ₂ , w ₀ )−w ₂ )/w ₂ ).

2. Define double-strand quality as the variation quality of a single sample.

3. According to the above steps 1 and 2, calculate the single-sample mutation quality of the tumor sample and the normal control sample, respectively. Then calculate TLOD. TLOD represents the probability that the variation is neither derived from the body (soma) nor from the embryo (germline), as follows:

First, calculate the output value of the following function.

TNreward(a ₁ , a ₂ , b ₁ , b ₂ )=maxH(TNrewardBinomial(a ₁ , a ₂ , b ₁ , b ₂ ), PLQTN(b ₂ /b ₁ , a ₂ /a ₁ ));

Among them, TNrewardBinomial(a ₁ , a ₂ , b ₁ , b ₂ )=10/log(10)×(a1×KL Bernoulli(b ₂ /b ₁ , a ₂ /a ₁ ))

If max(a,-a)<max(b,-b), then maxH(a,b)=a, otherwise maxH(a,b)=b;

PLQTN(t, n)=3×10×log ₁₀ (min(t/n, 10 ^2.5/3 )).

TNpenal(a ₁ , a ₂ , b ₁ , b ₂ )=max(0, min(DSVQnormal, PLQ(a ₂ /a ₁ ))-12.5×(max(0, OR(a ₁ , a ₂ , b ₁ , b ₂ )-1)) ² ); where DSVQnormal refers to the single-sample mutation quality of the normal sample, that is, the double-strand quality of the normal sample.

a ₁ , a ₂ , b ₁ and b ₂ have the following meanings:

a ₁ : The sum of the family quality of the reads supporting any allele in the paired normal control sample;

a ₂ : The sum of the family quality of the reads supporting the mutant allele in the paired normal control sample;

b ₁ : the sum of family quality of reads supporting any allele in the tumor sample;

b ₂ : The sum of the family quality of the reads supporting the mutant allele in the tumor sample.

Then, calculate TLOD, the calculation formula is as follows:

TLOD(a ₁ ,a ₂ ,b ₁ ,b ₂ )=min(DSVQtumor,PLQ(b ₂ /b ₁ ))+TNreward(a ₁ ,a ₂ ,b ₁ ,b ₂ )-TNpenal(a ₁ ,a ₂ , b ₁ , b ₂ ).

Wherein, PLQ(f)=90+3×10×log ₁₀ (f). DSVQtumor refers to the single-sample mutation quality of the tumor sample, that is, the double-strand quality of the tumor sample.

(C4) Calculate NLOD. NLOD represents the probability that the variation originates from the embryo (germline). details as follows:

First, for each sample (there are two samples of tumor and paired normal in total), the genotype likelihood (GL) of each genotype (GT) is calculated. The definitions and descriptions of GT and GL are in https://samtools.github.io/hts-specs/VCFv4.2.pdf. To sum up, homref (homozygous wild type), hetero (heterozygous type) and homalt (homozygous mutant type) are three types of GT. The calculation methods of the GL of these three GTs are as follows:

GL(f, homref) = max(GLPowerLaw(f, homref), GLBinomial(f, homref));

GL(f, hetero) = max(GLPowerLaw(f, hetero), GLBinomial(f, hetero));

GL(f, homalt) = max(GLPowerLaw(f, homalt), GLBinomial(f, homalt)).

Among them, GLPowerLaw(f, homref) = 10×3×min(0, log ₁₀ (0.02×(1-f)/f));

GLPowerLaw(f, hetero) = 10×3×min(0, log ₁₀ (min(f/(1-f), (1-f)/f)));

GLPowerLaw(f, homalt) = 10×3×min(0, log ₁₀ (0.02×f/(1-f)));

GLBinomial(f, homref) = 10/log(10) × d × KL Bernoulli(min(0.02, f), f);

GLBinomial(f, hetero)=10/log(10)×d×KL Bernoulli(0.5,f);

GLBinomial(f,homalt)=10/log(10)×d×KL Bernoulli(max(1-0.02,f),f);

KL Bernoulli(x,y)=y×log(y/x)+(1-y)×log((1-y)/(1-x));

If x=0, y=0, x is an infinitive or y is an infinitive, then KL Bernoulli(x,y)=0. For example, (0/0) is an infinitive.

And a refers to the number of reads supporting a mutant gene, d refers to the number of reads supporting any genotype, and f=a/d is the mutation frequency.

Then calculate the NLOD. NLOD=max(Qtumorhomref, Qnormalhomref).

where on the tumor sample, Qtumorhomref=-3-max(GL(f, homalt), GL(f, hetero)). Qtumorhomref is calculated using tumor sample data, so it has nothing to do with paired normal samples.

Wherein on the paired normal samples, Qnormalhomref=GL(f, homref)-max(GL(f, homalt), GL(f, hetero)). Qnormalhomref is calculated using paired normal sample data, so it has nothing to do with tumor samples.

5. Calculate the final mutation quality

Calculate the final mutation quality according to the following formula: min(TLOD, NLOD+31).

If the final mutation quality is at least 60, it is judged as positive and the mutant form is output; otherwise, it is judged as negative. A quality of 60 represents a false positive rate of one in a million. The mutation detection software MuTect (Cibulskis K, Lawrence M S, Carter S L, et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples[J]. Nature biotechnology, 2013, 31(3): 213.) is used indirectly The threshold of 60 is set, and the maximum comparison quality that BWA can output is 60. Therefore a threshold of 60 was used. The mutation output format is VCF. The VCF format definition is described in detail in (Danecek P, Auton A, Abecasis G, et al. The variant callformat and VCFtools[J]. Bioinformatics, 2011, 27(15): 2156-2158.).

Detect InDel:

Compared with SNV variation, the step of "base quality estimation" is added before "base quality correction", and the rest are the same as SNV variation detection methods. The Base Quality Estimation method is as follows:

InDelQual(x, y, a, b, c) = min(x+10, y+10, STRQual(a, b))+10×log ₁₀ (c)

Among them, InDelQual is the InDel quality that has the same purpose as the base quality, x and y are the base quality of the previous base and the next base where InDel occurs, respectively, and a and b are the units of the short tandem repeat (STR) The size (repeat size) and the number of units (repeat number), c is the number of bases of InDel. The STRQual function is defined as follows:

STRQual(a, b)=44-10×log ₁₀ (1+8×softplus(a×b-8)/(θ×a ² ))

The softplus function is defined as follows:

softplus(v)=log(1+exp(v)), where log refers to the natural logarithm.

If InDel is a deletion of a STR unit, then θ=0.25, otherwise θ=1.

Example 2. Application examples of the gene variation detection method based on high-throughput sequencing of the present invention

One, the comparison of the method of the present invention and smCounter2 method

Adopt the method that the embodiment of the present invention 1 establishes and existing smCounter2 method (reference "Xu C, Gu X, Padmanabhan R, et al.smCounter2: an accurate low-frequency variant caller for targeted sequencing data with unique molecular identifiers [J] .Bioinformatics, 2019, 35(8):1299-1309.") Detected SNV variation on N13532 (SRR7526729) and N0261 (SRR7526728) samples, respectively.

N13532 (SRR7526729) sample and N0261 (SRR7526728) sample: data can be found at the following English websites. https://www.ncbi.nlm.nih.gov/sra/SRR7526728; https://www.ncbi.nlm.nih.gov/sra/SRR7526729. In addition, it is also described in the smCounter2 article. SRR7526728 product number: CDHS-13907Z-562; SRR7526729 product number: CDHS-13532Z-10181.

Figures 1 to 4 show the precision and recall rates of the present invention and smCounter2 on N13532 (SRR7526729) and N0261 (SRR7526728) samples (item number: SRP153933). As can be seen from the figure, the advantages of the present invention are higher F1 score and area under the curve.

smCounter2 needs to do machine learning training by repeatedly sequencing the same sample and obtaining some (but not all) known positive mutations of the sample in advance, but the present invention does not require these, so the present invention reduces the requirement for existing data . For the molecular label data of blood samples (Fig. 1 to Fig. 4), the present invention improves the true positive rate and reduces the false positive rate.

Two, the performance comparison of the inventive method, the Mutect2 method and the Strelka2 method on the standard mixed sample

The method established in Example 1 of the present invention and the existing Mutect2 method and Strelka2 method are used to detect SNV variation on standard mixed samples.

Strelka2 method: see "Kim S, Scheffler K, Halpern A L, et al. Strelka2: fast and accurate calling of germline and somatic variants [J]. Nature methods, 2018, 15(8): 591." Strelka2 source code is at https://github.com/Illumina/strelka.

Mutect2 method: see "Cibulskis K, Lawrence M S, Carter S L, et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancersamples[J]. Nature biotechnology, 2013, 31(3): 213." Mutect2 source code is at https://github.com/broadinstitute/gatk .

The standard mixed sample contains virtual tumor and normal WGS (whole genome) data: tumor is a mixture of two samples of NA12878/NA24385 with an average sequencing depth of 90×, and the normal is a mixture of NA12878/NA24385 samples with an average sequencing depth of 30× A single sample of pure NA24385. Standard mixture sample data are at: ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/use_cases/mixtures .

The germline variants of NA12878 are all in: ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/release/NA12878_HG001/NISTv3.3.2/GRCh37/HG001_GRCh37_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_ CHROM1-X_v .3.3.2_highconf_PGandRTGphasetransfer.vcf.gz.

The descriptive documentation for these variants is at:

ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/release/NA12878_HG001/NISTv3.3.2/README_NISTv3.3.2.txt.

Figure 5 and Figure 6 show the performance test results of the method of the present invention, the Mutect2 method and the Strelka2 method on standard mixed samples. As can be seen from the figure, the advantages of the present invention are higher F1 score and area under the curve.

Strelka2 needs to do machine learning training by sequencing other similar samples and obtaining all known positive mutations of other similar samples in advance, but the present invention does not need these. In summary, the present invention increases the true positive rate and reduces the false positive rate for WGS data of tissue samples, and potentially reduces biological testing.

The results of this example show that the method established in Example 1 of the present invention (abbreviated as UVC) is used, and compared and evaluated with other software. The evaluation showed that UVC performed very well: UVC both increased the true positive rate and decreased the false positive rate, regardless of whether the data were molecularly labeled or not.

Claims (8)

1. A gene variation detection method based on high-throughput sequencing comprises the following steps:

(A) Preprocessing the high-throughput sequencing data, mainly comprising base quality correction and read de-duplication;

(B) Calculating the preference generated by sequencing, and reducing the number of effective reads supporting mutation through the preference to obtain the corrected read number supporting mutation;

(C) Calculating the relative quality of all single samples by using the corrected read number of the support mutant, then obtaining the variation quality of the single sample, and judging through the variation quality;

the step (C) is carried out as follows:

(C1) For each site where quality needs to be calculated, calculating family quality, single-strand quality, double-strand quality by using the following formulas;

Quality(p，B，b，r)＝logr(BIAS(B，p×B，B，b)×(r-1)+1)×10×log ₁₀ (OR(B，p×B，B，b))；

wherein Bias (B, pxb, B) = max (0, or (LS (B, pxb, B)) -1) × (B ÷ (B + B)) +1;

OR represents an odds ratio; the OR is calculated as follows:

OR(B，p×B，B，b)＝(B×b)÷(p×B×B)；

LS stands for Laplace smoothing; LS formula is as follows:

LS(B，p×B，B，b)＝Laplace ₂ (Laplace ₁ (B，p×B，B，b))；

laplace1 is expressed as: laplace ₁ (B，p×B，B，b)＝(B+p′，p×B+p′，B+p′，b+p′)；

Laplace2 is expressed as: laplace ₂ (B，p×B，B，b)＝(B+p′×α，p×B+p′×α，B+p′×β，b+p′×β)；

In Laplace ₁ Formula and Laplace ₂ In the formula, the default value of p' is 0.5;

in Laplace ₂ In the formulas, α and β are as follows:

α＝max(0，p′×((B+p×B)/(B+b)-1))；

β＝max(0，p′×(B+b)/(B+p×B)-1)；

wherein:

(1) Let r =2 and calculate the family quality using the following definition for reads in a family:

p: error probability corresponding to the base quality threshold calculated in the step (A) reads deduplication;

b: the number of reads that actually observed to support any allele from the data;

b: actual observations from the data support the number of reads for the mutant allele;

(2) R =1.5 if the high throughput sequencing experiment is banked using PCR amplification, r =1.25 if the high throughput sequencing experiment is banked using capture or WGS; comparing that the families on one chromosome chain are on one single-strand; for a single-strand, the single-strand quality is calculated using the following definitions:

p: the error probability corresponding to each possible family quality;

b: actual observation of family numbers supporting any allele;

b: actual observation of family numbers supporting mutant alleles;

if the single-strand quality is greater than the maximum similarity quality value, then let the single-strand quality equal the maximum similarity quality value; the maximum family quality is as follows: for C: g > T: the a mutant form was 44, for T: a > C: the G mutant form was 48, 52 for the other SNV mutant forms, 60 for all InDel mutant forms;

after all p's are tried, the maximum single-strand quality is taken as the final single-strand quality;

(3) Forward or reverse direction of each single-strand on a double-strand; for a single-strand, double-strand quality was calculated using the following definitions:

v ₁ ′＝v ₁ +v ₂ ×min(1，(min(w ₁ +w ₂ ，w ₀ )-w ₁ )/w ₁ )；

wherein v is ₁ ' is the final variable mass on a strand, v ₁ Is the variable quality of this chain, i.e.the final single-strand quality, v ₂ Is the quality of variation on the other strand, w ₁ Is the average value of the family quality, w, of this chain ₂ Is the average value of the family quality of the other chain, w ₀ Equal to the maximum similarity quality, which is as follows, plus 10: for C: g > T: the a mutant form was 44, for T: a > C: the G mutant form was 48, 52 for the other SNV mutant forms, and 60 for all InDel mutant forms;

exchange v ₁ And v ₂ And exchange w ₁ And w ₂ According to the calculation of v ₁ ' symmetry calculation of formula v ₂ ′；

v ₂ ′＝v ₂ +v ₁ ×min(1，(min(w ₁ +w ₂ ，w ₀ )-w ₂ )/w ₂ )；

(C2) Defining double-strand quality as a single sample variation quality;

(C3) Calculating the single sample mutation quality of the tumor sample and the normal control sample according to the (C1) and the (C2) respectively; then calculating TLOD; TLOD represents the probability that the variation is neither somatic nor embryonic, as follows:

firstly, calculating the output values of the following functions;

TNreward(a ₁ ，a ₂ ，b ₁ ，b ₂ )＝maxH(TNrewardBinomial(a ₁ ，a ₂ ，b ₁ ，b ₂ )，PLQTN(b ₂ /b ₁ ，a ₂ /a ₁ ))；

wherein, TNrewardBonomial (a) ₁ ，a ₂ ，b ₁ ，b ₂ )＝10/log(10)×(a ₁ ×KLBernoulli(b ₂ /b ₁ ，a ₂ /a ₁ ) ); wherein, KLBernoulli (b) ₂ /b ₁ ，a ₂ /a ₁ )＝(a ₂ /a ₁ )×log((a ₂ /a ₁ )/(b ₂ /b ₁ ))+(1-(a ₂ /a ₁ ))×log((1-(a ₂ /a ₁ ))/(1-(b ₂ /b ₁ ) ); if b is ₂ /b ₁ ＝0、a ₂ /a ₁ ＝0、b ₂ /b ₁ Is an indefinite form or a ₂ /a ₁ Is an infinite form, then KLBernoulli (b) ₂ /b ₁ ，a ₂ /a ₁ ) =0; maxH (a, b) = a if max (a, -a) < max (b, -b), otherwise maxH (a, b) = b;

PLQTN(t，n)＝3×10×log ₁₀ (min(t/n，10 ^2.5/3 ))；

TNpenal(a ₁ ，a ₂ ，b ₁ ，b ₂ )＝max(0，min(DSVQnormal，PLQ(a ₂ /a ₁ ))-12.5×(max(0，OR(a ₁ ，a ₂ ，b ₁ ，b ₂ )-1)) ² ) (ii) a Wherein DSVQnormal refers to the single sample mutation quality of normal samples, namely double-strand quality of normal samples; OR represents the ratio of ratios, and the calculation formula of OR is as follows: OR (a) ₁ ，a ₂ ，b ₁ ，b ₂ )＝(a ₁ ×b ₂ )÷(a ₂ ×b ₁ )；a ₁ 、a ₂ 、b ₁ And b ₂ The meanings are as follows:

a ₁ : family q of reads supporting any allele in the paired normal control sample(ii) a total of the qualitys;

a ₂ : the sum of the family quality of reads supporting mutant alleles in the paired normal control sample;

b ₁ : family quality sum of reads supporting any allele in the tomor sample;

b ₂ : family quality sum of reads supporting mutant alleles in the tomor sample;

then, TLOD is calculated as follows:

TLOD(a ₁ ，a ₂ ，b ₁ ，b ₂ )＝min(DSVQtumor，PLQ(b ₂ /b ₁ ))+TNreward(a ₁ ，a ₂ ，b ₁ ，b ₂ )-TNpenal(a ₁ ，a ₂ ，b ₁ ，b ₂ )；

wherein PLQ (f) =90+3 × 10 × log ₁₀ (f) (ii) a DSVQtomor refers to the single-sample mutation quality of a tomor sample, i.e., the double-strand quality of a tomor sample;

(C4) Calculating NLOD; NLOD represents the probability of coming from an embryo for variation, as follows:

firstly, calculating the genotype likelihood of each genotype aiming at each sample; GT for genetic type, GL for genetic type likelihood for short; the GL of three GTs, homozygous wild-type, heterozygous and homozygous mutant, was calculated as follows:

GL formula for homozygous wild type: GL (f, homref) = max (GLPowerLaw (f, homref), GLBinomial (f, homref));

GL calculation formula for heterozygote type: GL (f, hetero) = max (GLPowerLaw (f, hetero), GLBinomial (f, hetero));

GL calculation formula for homozygous mutant: GL (f, homalt) = max (GLPowerLaw (f, homalt), GLBinomial (f, homalt));

where GLPowerLaw (f, homref) =10 × 3 × min (0,log) ₁₀ (0.02×(1-f)/f))；

GLPowerLaw(f,hetero)＝10×3×min(0,log ₁₀ (min(f/(1-f),(1-f)/f)))；

GLPowerLaw(f,homalt)＝10×3×min(0,log ₁₀ (0.02×f/(1-f)))；

GLBinomial(f,homref)＝10/log(10)×d×KLBernoulli(min(0.02,f),f)；

GLBinomial(f,hetero)＝10/log(10)×d×KLBernoulli(0.5,f)；

GLBinomial(f,homalt)＝10/log(10)×d×KLBernoulli(max(1–0.02,f),f)；

KLBernoulli (x, y) = y × log (y/x) + (1-y) × log ((1-y)/(1-x)); KLBernoulli (x, y) =0 if x =0, y =0, x is indeterminate or y is indeterminate;

and a refers to the number of reads supporting mutant genes, d refers to the number of reads supporting any genotype, f = a/d is the mutation frequency;

then, calculating NLOD; NLOD = max (qtumhorref, qnormalomref);

wherein, on the tomor sample, qtomormref = -3-max (GL (f, homalt), GL (f, hetero));

wherein Qnormalhomref = GL (f, homref) -max (GL (f, homalt), GL (f, hetero));

(C5) The final quality of variation was calculated according to the following formula: min (TLOD, NLOD + 31);

if the final quality of variation is at least 60, determining as positive and outputting a variant form; otherwise, the result is judged to be negative.

2. The method of claim 1, wherein: in the method, the genetic variation detected is an SNV variation or an InDel variation.

3. The method of claim 1, wherein: when the detected genetic variation is an InDel variation, the step (A) further comprises a step of estimating the base quality before the base quality correction; the base quality estimation is performed as follows:

InDelQual(x,y,a,b,c)＝min(x+10,y+10,STRQual(a,b))+10×log ₁₀ (c)；

wherein InDelQual is InDel mass having the same purpose as base mass, x and y are base mass of a previous base and a next base of the occurrence of InDel, a and b are unit size and unit number of the short tandem repeat sequence, and c is base number of InDel; the STRGUAL function is defined as follows:

STRQual(a,b)＝44-10×log ₁₀ (1+8×softplus(a×b-8)/(θ×a ² ))；

the softplus function is defined as follows:

softplus (v) = log (1 + exp (v)), where log refers to the natural logarithm;

if InDel is a deletion of one short tandem repeat unit, then θ =0.25, and in addition θ =1.

4. A method according to any one of claims 1-3, characterized in that: in the step (A), the base quality correction is carried out as follows:

if the high throughput sequencing is double-ended sequencing, merging before de-duplication and filtering, merging the R1 end and the R2 end, and then using the base quality of the merged reads;

for the overlapping region of the R1 end and the R2 end, the base quality after combination is calculated as follows:

when the bases corresponding to p and q are the same, the following formula is used: max (p, q);

when the bases corresponding to p and q are different, the following formula is used for calculation: max (p, q) -min (p, q);

wherein p is the base mass before the end R1 is corrected, and q is the base mass before the end R2 is corrected;

for base masses in non-overlapping regions of the R1 end and the R2 end, the combined base masses are equal to the base masses before correction;

if the single-ended sequencing is carried out, merging is not needed, and the base quality after correction is equal to the base quality before correction;

and/or

In the step (A), read deduplication is performed as follows:

placing PCR repeats or optically repeated reads into a family to remove repeats;

further, PCR repeats or optically repeated reads are placed in a family to remove duplicates according to a method comprising the following steps:

(a1) For each locus, the number of reads that start or stop at that locus, abbreviated as end depth, is found;

(a2) Each locus x attracts other loci y, the attraction strength s of x to y is calculated, the attraction strength s is in direct proportion to the depth of the x end, in inverse proportion to the depth of the y end and in exponential decay relation with the number z of bases between x and y;

the formula of the attractive force intensity s is as follows;

s＝x/y×5 ^-z ；

if s > 1, then y is considered a PCR or optical repeat of x, thereby placing x and y in a family;

(a3) And (3) filtering: the threshold for base mass was calculated as follows:

Threshold(P)＝argmax _p∈P (p ^{-l{q∈P|q≤p}|} )；

wherein P is a set of sequencing error probabilities corresponding to all base masses; the mathematical symbol e represents set attribution, so that the meaning of P e P means that P is any one element in P, and the meaning of q e P means that q is any one element in P;

bases on reads were filtered according to base quality thresholds: filtering out bases whose base quality is lower than a base quality threshold after correction;

after filtering, for each family, taking the base with the most frequent family as the base after de-weighting; the repeated reads are integrated into a consistent read.

5. A method according to any one of claims 1-3, characterized in that: in the step (B), four preferences of depletion bias, position bias, strand bias and mismatch bias are calculated in each chromosome chain; then selecting the highest preference as the integral preference of the chromosome chain; the number of reads supporting mutant is corrected, the corrected number of reads being equal to the number of reads before correction divided by the overall preference.

6. The method of claim 5, wherein: the four preferences of depletion bias, position bias, strand bias and mismatch bias are calculated as follows:

Bias(a ₁ ，a ₂ ，b ₁ ，b ₂ )＝max(0，OR(LS(a ₁ ，a ₂ ，b ₁ ，b ₂ ))-1)*(b ₂ ÷(b ₂ +b ₁ ))+1；

wherein OR represents an odds ratio; the OR is calculated as follows:

OR(a ₁ ，a ₂ ，b ₁ ，b ₂ )＝(a ₁ ×b ₂ )÷(a ₂ ×b ₁ )；

wherein LS represents Laplace smoothing; the LS formula is as follows: LS (a, b, c, d) = Laplace ₂ (Laplace ₁ (a, b, c, d)); wherein, laplace ₁ The formula is as follows: laplace ₁ (a, b, c, d) = (a + p, b + p, c + p, d + p); laplace2 has the following formula: laplace ₂ (a，b，c，d)＝(a+p×α，b+p×α，c+p×β，d+p×β)；

In Laplace ₁ Formula and Laplace ₂ In the formula, the default value of p is 0.5;

in Laplace ₂ In the formulas, α and β are as follows: α = max (0,p × ((a) ₁ +a ₂ )/(b ₁ +b ₂ )-1))；β＝max(0，p×(b ₁ +b ₂ )/(a ₁ +a ₂ )-1)；

For the reduction bias, the meaning of the parameters is as follows:

a ₁ : the number of reads supporting any allele that were deduplicated;

a ₂ : after deduplication, the number of reads supporting any allele;

b ₁ : the number of reads that support the mutant allele that are deduplicated;

b ₂ : number of reads supporting mutant alleles after deduplication;

for position bias, the meaning of each parameter is as follows:

for each direction, extending 1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66 bases in direction; the following four values are then substituted:

a ₁ : the number of reads that support any allele within reads that have spanned the extended locus;

a ₂ : the number of reads that support any allele in reads that do not span the extended locus;

b ₁ : in reads that have spanned the extended locus, the number of reads supporting the mutant allele;

b ₂ : number of reads that support mutant alleles in reads that do not span the extended locus;

sequentially substituting 11 extension units to generate a sequence consisting of 11 numerical values; taking every 3 adjacent values in the sequence as an array, thereby generating 9 arrays; taking the minimum value for each array, thereby producing 9 minimum values; taking the maximum value of the 9 minimum values as position bias;

for strand bias, the following values are substituted:

a ₁ : in the forward direction of the chromosome strand S, the number of reads supporting any allele;

a ₂ : in the reverse of chromosome chain S, the number of reads supporting any allele;

b ₁ : in the forward direction of chromosome strand S, the number of reads supporting the mutant allele;

b ₂ : in the reverse direction of the chromosome strand S, the number of reads supporting the mutant allele;

simultaneously calculating base quality inventory and strand position inventory, wherein the base quality inventory is BQI for short, the formula is represented by I, the formula comprises 1I, the strand position inventory is SPI for short, the formula is represented by I and comprises 2I, the numerical values of the 3I are calculated by the following method, and then the minimum numerical value is taken; then, dividing the strand bias by the minimum value to normalize the strand bias; finally ensuring that the minimum Strand bias is at least 1;

for BQI, c is the average base of all alleles of the other strandBase mass, d is the average base mass of the mutant allele of this strand, and I =10 ^d-c ；

For SPI, c is the average base distance to the end of read for all alleles of the other strand, d is the average base distance to the end of read for mutant alleles suspected of having a preferred strand, and

each read has two ends: left and right sides; the left and right termini thus yield two I, respectively;

let x =1,2,3, …,11 for mismatch bias; x is the mismatch number threshold, 1,2, …,11 are different attempted thresholds;

for each value of x, the following four values are substituted:

a ₁ : mismatch numbers < x, number of reads supporting any allele;

a ₂ : the mismatch number is more than or equal to x, and the number of reads of any allele is supported;

b ₁ : mismatch numbers < x, number of reads supporting mutant alleles;

b ₂ : the mismatch number is more than or equal to x, and supports the reads number of mutant alleles;

sequentially substituting 11 x values to generate a sequence consisting of 11 values; taking every 5 adjacent values in the sequence as an array, thereby generating 7 arrays; taking the minimum value for each array, thereby generating 7 minimum values; the maximum of the 7 minima was taken as the mismatch bias.

7. A gene variation detection system based on high-throughput sequencing comprises a device A, a device B and a device C;

the device A is capable of pre-processing high throughput sequencing data, mainly including base quality correction, read de-duplication, according to step (A) of any one of claims 1 to 6;

said device B is capable of calculating the preference generated by sequencing according to the step (B) of any one of claims 1 to 6, and obtaining the corrected read number of the support mutation type by reducing the effective read number of the support mutation type according to the preference;

the device C can calculate the correlation quality of all the individual samples by using the corrected read numbers of the supporting mutants according to the step (C) of any one of claims 1 to 6, and then obtain the variation quality of the individual samples, and perform the judgment by the variation quality.

8. Use of the method of any one of claims 1-6 or the system of claim 7 for the manufacture of a product for tumor pre-screening, tumor prognosis, tumor classification and/or tumor medication guidance.

Images