A method for reviewing high-throughput sequencing gene variation detection results

Abstract

The present invention discloses a method for reviewing the results of high-throughput sequencing gene variation detection. The method comprises the following steps: constructing a training set, including sequencing data of a number of positive variation sites and negative non-variation sites; extracting and vectorizing the variation site features; constructing a model using the random forest method, and then using the model to determine whether the variation at the site to be tested is a variation site. The method of the present invention can quickly and accurately complete the process of manually reviewing mutations, greatly saving labor costs and improving the overall variation detection efficiency; achieving or even exceeding the accuracy of manual review to a certain extent, and greatly reducing the subjectivity of manual inspection; using machine learning algorithms, it has certain innovations in this field; developing a model based on real clinical data, which is closer to real cases, and performs the same as manual results on the gold standard data set, and can be adapted to various different cell variation detection software, with good practicality.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics, and in particular to a method for reviewing high-throughput sequencing gene variation detection results. Background Art

[0002] Currently, second-generation sequencing (NGS) technology has been widely used in somatic variant detection, providing an essential foundation for clinical medication guidance and prognosis monitoring. Therefore, accurate variant detection has become a prerequisite for precision medicine. However, accurately detecting somatic variants (especially low-frequency variants) faces significant challenges, including sample preparation damage, sequencing errors, and misalignment with the reference genome. Bioinformatics software for somatic variant detection is a downstream and crucial component. Numerous excellent variant detection software (such as Mutect, Strelka, and Varscan) are widely used in the somatic variant detection market. However, when faced with complex variants, commonly used variant detection software often produces a certain percentage of errors (ranging from 15% to 30%). Therefore, further manual review is required using visual variant visualization software (e.g., the most commonly used software, IGV (Integrative Genome Viewer) [http: / / www.igv.org / ]) to eliminate false positives and present the final true variant results. However, this manual review process has three limitations: First, using visualization software to review variants and generate accurate results requires extensive bioinformatics expertise and extensive training, and is also subject to significant subjective judgment. Second, as sequencing data volumes grow, the efficiency and labor costs of manual review hinder the timeliness and cost-effectiveness of variant detection. Third, the subjectivity of manual review and the lack of uniformity in review standards also pose a threat to variant detection accuracy. Due to the significant subjectivity of IGV's review rules, Erica K. et al. systematically investigated and tested the standardization of its criteria ("Standard operating procedure for somatic variant refinement of sequencing data with paired tumor and normal samples" DOI:10.1038 / s41436-018-0278-z). This study established and validated a set of criteria for manual variant review using a set of gold-standard datasets. By learning this set of criteria, researchers improved the accuracy of somatic variant calls from 77.4% to 94.1%. Although the IGV review rules have been professionally established and verified, the problem of low efficiency of manual review has not been effectively solved. Summary of the Invention

[0003] In response to the above problems, the purpose of the present invention is to use machine learning algorithms and manually annotated gold standard data sets to provide a method for replacing or automating manual review of somatic variants, that is, a method for reviewing gene variation detection results.

[0004] In a first aspect, the present invention claims protection for a method for reviewing the results of high-throughput sequencing gene variation detection.

[0005] The method for reviewing high-throughput sequencing gene variation detection results claimed in the present invention may include the following steps:

[0006] (A) A training set was constructed, including sequencing data of several positive variant sites and negative non-variant sites.

[0007] (B) extracting and vectorizing features of the variant site from the training set; the features of the variant site include any 6 or more of the following, such as any 7, 8, 9, 10, or 11:

[0008] Mutation support number: the number of reads in tumor tissue that support the mutation site to be detected;

[0009] Mutation frequency: the frequency of the mutation site to be detected in tumor tissue;

[0010] Base quality: the average quality of the bases at the mutation site to be detected in the reads supporting the mutation site to be detected in the tumor tissue;

[0011] Misalignment rate: The average misalignment rate of the 50 bp upstream and downstream sequences of the mutation site to be detected in the reads supporting the mutation site to be detected in tumor tissue;

[0012] HDR value: the score of suspected homology alignment errors in reads supporting the mutation site to be detected in tumor tissue;

[0013] SideB value: edge preference score of reads supporting the mutation site to be detected in tumor tissue;

[0014] StrandB value: the strand preference score of reads supporting the mutation site to be detected in tumor tissue;

[0015] Base ratio: the ratio of the number of reads supporting the mutation site to be detected in normal control tissue and tumor tissue;

[0016] Other variant type ratio: If there are other variant types in the tumor tissue besides a certain variant type, the ratio of the number of supporting reads of the other variant types to the certain variant type;

[0017] Alignment length: the average length of reads supporting the mutation site to be detected in tumor tissue;

[0018] InDel length: the length of the deletion or insertion of Indel;

[0019] (C) Using the feature results obtained in step (B), a random forest method is used to construct a model, and then the model is used to determine whether the site to be tested is a variant site.

[0020] In the above method, the characteristics of the variant site may further include any one or more of the following, for example, any two, three, four, five or six:

[0021] Alignment quality: the average alignment quality of reads supporting the mutation site to be detected in tumor tissue;

[0022] Environmental quality: The average base quality of the 50 bp sequences upstream and downstream of the mutation site to be detected in the reads supporting the mutation site to be detected in the tumor tissue;

[0023] Insert size: the mean insert size of reads supporting the mutation site to be detected in tumor tissue;

[0024] Genome complexity score: The complexity score of the reference genome sequence 20 bp upstream and downstream of the mutation site to be detected;

[0025] Normal coverage depth: the coverage depth of the mutation site to be detected in normal control tissue;

[0026] Normal InDel Presence Rate: Whether an InDel is present 50 bp upstream and downstream of the mutation site in the tumor tissue and / or normal tissue. If so, the product of the InDel's variation frequency and length is calculated, or the products are further added together. If not, the value is 0. For example, if an InDel is present 50 bp upstream and downstream of the mutation site in the tumor tissue or normal tissue, the product of the InDel's variation frequency and length is calculated. If an InDel is present 50 bp upstream and downstream of the mutation site in both the tumor tissue and normal tissue, the product of the InDel's variation frequency and length is calculated and then added together.

[0027] In the above method, the type of genetic variation is base substitution or InDel; the model is constructed based on the variation type, using sequencing data of variant sites of the same variation type to construct a training set, extracting and vectorizing features, and constructing the model. For example, a training set is constructed using sequencing data of variant sites with base substitution variation type, and after extracting and vectorizing features, a model is constructed.

[0028] In a specific embodiment of the present invention, when the type of genetic variation is base substitution, the characteristics of the variation site preferably include: mutation support number, mutation frequency, base quality, mismatch rate, HDR value, SideB value, StrandB value and base ratio; more preferably, include mutation support number, mutation frequency, base quality, mismatch rate, HDR value, SideB value, StrandB value, base ratio and alignment quality.

[0029] In a specific embodiment of the present invention, when the type of the genetic variation is InDel, the characteristics of the variation site preferably include mutation support number, mutation frequency, false alignment rate, HDR value, SideB value, StrandB value, base ratio, other variation type ratio, InDel length and alignment quality.

[0030] The base substitution may be a single base (SNV), a double base (DNV) and / or a triple base (TNV) substitution.

[0031] In the above method, the HDR value is calculated as follows:

[0032] When the number of reads supporting the mutation site to be detected is less than 5, or when the number of reads supporting the mutation site to be detected is greater than or equal to 5 and the absolute value of the difference in mutation frequency between a certain variant upstream and downstream of the reads supporting the mutation site to be detected in normal control tissue and tumor tissue is less than 0.7, the HDR value is 0;

[0033] When the number of reads supporting the mutation site to be detected is greater than or equal to 5, and the absolute value of the difference between the mutation frequencies of a certain variant upstream and downstream of the reads supporting the mutation site to be detected in normal control tissue and tumor tissue is greater than or equal to 0.7, the HDR value is calculated as follows:

[0034]

[0035] Among them, I i Equal to 1; n is the number of sites that meet the HDR conditions.

[0036] The smaller the HDR value, the higher the probability of alignment error of the suspected homology.

[0037] In the above method, the SideB value is calculated as follows:

[0038] When the number of reads supporting the mutation site to be detected is greater than or equal to 10, the SideB value is calculated according to the following formula:

[0039]

[0040] Where n = 50, Dti is the depth of supporting reads for the i-th bp upstream of the mutation site; Dn i is the depth of supporting reads at the ith bp downstream of the mutation site. abs represents the absolute value.

[0041] When the number of reads supporting the mutation site to be detected is less than 10, the SideB value is 0.

[0042] The lower the SideB score, the more marginal preference the site has.

[0043] In the above method, the StrandB value is calculated as follows:

[0044] When the number of reads supporting the mutation site to be detected is greater than or equal to 10, the StrandB value is calculated according to the following formula:

[0045]

[0046] Where R is the ratio of positive strands in reads supporting the mutation site to be detected. abs represents the absolute value.

[0047] When the number of reads supporting the mutation site to be detected is less than 10, the StrandB value is 0.

[0048] The higher the StrandB value, the less strand preference the mutation has.

[0049] Furthermore, the complexity score is calculated as follows: if 6 or more consecutive single-base repeat sequences or tandem repeat sequences (for example, 'AAAAAA' or 'ATCATCATCATCATCATCATC') appear within 20 bp upstream and downstream of the mutation site to be detected on the reference genome sequence, the feature value is 1, otherwise it is 0; the consecutive single-base repeat sequence means that the same base appears repeatedly continuously, for example, 'AAAAAA'; in the consecutive tandem repeat sequence, the tandem sequence refers to a sequence of two or more bases in series, for example, 'AT' or 'ATC', and the consecutive tandem repeat sequence is, for example, 'ATCATCATCATCATCATC'.

[0050] Before extracting and vectorizing features of variant sites from the training set in step (B), the method may further include filtering reads in the training set that support the mutation site to be detected as follows: filtering out reads with at least one of the following three conditions: read length less than 55, read alignment quality less than 10, and soft clip length greater than 20 bp.

[0051] In a specific embodiment of the present invention, pysam (v0.11) of python 2.7.8 is used to filter and exclude reads in paired samples of candidate variant sites.

[0052] Through the above steps, unreliable reads can be filtered out to prevent the introduction of unnecessary noise into the input training set.

[0053] Furthermore, the model can be constructed using the randomForest package of R language (such as the randomForest_4.16-14 package in R3.5.1).

[0054] Furthermore, the ten-fold cross-validation method is used to optimize the main parameters of the model (such as the number of trees (ntree) and the number of randomly extracted features (mtry) to construct the gradient).

[0055] More specifically, for base substitution mutations, the optimized main parameters are as follows: ntree is set to 900; mtry is set to 5; for InDel mutations, the optimized parameters are as follows: ntree is set to 2000; mtry is set to 4.

[0056] In step (A) of the method, the positive variant sites and negative non-variable sites refer to positive variant sites or negative non-variable sites that have been detected according to existing software, annotated, and manually reviewed.

[0057] The detection includes detecting sequencing data using existing variation detection software, for example, using Mutect1 (version 3.1) in GATK3.1 for base substitution variation, and using Strelka (v1.0.14) for InDel variation (insertion and deletion) (the above two software are used here for variation detection, but the algorithm is not limited to the detection results of the above two software in terms of applicability).

[0058] After variant detection, the variant results are annotated using, for example, VEP (v83).

[0059] After annotating the variant results, the following filtering conditions can be used to filter candidate variants: filter variants outside the WES (whole exome) region; filter out synonymous mutations; filter out variants with a mutation frequency less than 5% against non-hotspot variants in the COSMIC database (https: / / cancer.sanger.ac.uk / cosmic / ); for hotspot variants included in the COSMIC database, filter out variants with a mutation frequency less than 1%; filter out variants with a site mutation frequency greater than 1% in the commonly used 1000 gene database (suspected common variants); filter out low-coverage sites (specifically, site coverage less than 40X in tumor, or less than 10X in normal); filter out variants with an absolute support number less than or equal to 6 reads.

[0060] Among them, the variant sites with a site variation frequency greater than 1% in the commonly used 1000-gene database include:

[0061] Variant sites with a population frequency greater than 1% in the NHLBI-ESP database (https: / / evs.gs.washington.edu / EVS / );

[0062] Variant sites with a population frequency greater than 1% in the 1000 Genomes Project database (http: / / phase3browser.1000genomes.org / index.html);

[0063] Variants with a frequency greater than 1% in East Asian populations from the 1000 Genomes Project data (http: / / phase3browser.1000genomes.org / index.html).

[0064] After the above filtering steps, a large number of low-frequency variants that are difficult to identify and variants with unclear clinical significance will be filtered out.

[0065] Finally, true and false variants are obtained through manual review as positive variant sites and negative non-variant sites; for example, IGV is used to evaluate the true and false labels of candidate variants.

[0066] The gene variation refers to the gene variation occurring in human tumor tissue compared with normal control tissue, which may be a base substitution or InDel.

[0067] In the present invention, the base substitution includes SNV variation, DNP variation and TNP variation, that is, single base variation, double base variation and / or triple base variation.

[0068] In a second aspect, the present invention claims protection for a system for reviewing the results of high-throughput sequencing gene variation detection.

[0069] The system for reviewing gene variation detection results claimed in the present invention comprises device A, device B, and device C;

[0070] The device A is capable of constructing a training set of sequencing data including a plurality of positive variant sites and negative non-variant sites according to step (A) of the method described above;

[0071] The device B can extract and quantify the features of the variant site from the training set according to step (B) of the method described above;

[0072] The device C can extract and vectorize the features of the variant site using the device B according to step (C) of the method described above, construct a model using the random forest method, and then use the model to determine whether the site to be tested is a variant site.

[0073] In a third aspect, the present invention claims protection for any of the following applications:

[0074] (I) Use of the method according to the first aspect or the system according to the second aspect in the preparation of a product for early tumor screening, tumor prognosis, tumor classification and / or tumor medication guidance;

[0075] The product can review the results of high-throughput sequencing gene variation detection according to the steps of the method described in the first aspect above.

[0076] (II) Use of the method described in the first aspect or the system described in the second aspect in early tumor screening, tumor prognosis, tumor classification, and / or tumor medication guidance;

[0077] (III) Use of the method described in the first aspect or the system described in the second aspect in detecting gene mutations.

[0078] The present invention has the following advantages due to the adoption of the above technical solution:

[0079] (1) The algorithm can quickly and accurately complete the process that originally required manual mutation inspection, which can greatly save labor costs and improve the efficiency of the overall mutation detection process;

[0080] (2) The algorithm surpasses the accuracy of manual review to a certain extent and can greatly reduce the subjectivity of manual inspection;

[0081] (3) The algorithm uses machine learning algorithms, which is innovative in this field;

[0082] (4) The algorithm is based on real clinical data to develop a model that is closer to real cases. Its performance on the gold standard dataset is comparable to that of manual results. It can also be adapted to various cell variation detection software and has good practicality. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] Figure 1 This is an algorithm flow chart of the method for reviewing gene variation detection results of the present invention.

[0084] Figure 2 The error rate curve calculated by cross-validation on the training set when screening features.

[0085] Figure 3 It is the AUC indicator of base substitution and InDel models.

[0086] Figure 4 Prediction of base substitution and InDel models on an independent validation set. DETAILED DESCRIPTION

[0087] Unless otherwise specified, the experimental methods used in the following examples are conventional methods.

[0088] Unless otherwise specified, the materials and reagents used in the following examples can be obtained from commercial sources.

[0089] Example 1: Establishment and application of the method for reviewing gene mutation detection results of the present invention

[0090] According to Figure 1 As shown in the flowchart, the model building method includes the following steps:

[0091] 1. Training Data Preparation

[0092] In order to automate the process of manual review of somatic variants (tumor tissue and normal control tissue samples) using a supervised learning algorithm (random forest), 94 cancer whole-exome sequencing samples were selected and sequenced by Illumina with an average depth of approximately Tumor 200X and Normal 100X.

[0093] The sequencing fastq files were aligned using the alignment software bwa-0.7.10, and then the output bam data was tested for variation. For example, Mutect1 (version 3.1) in GATK3.1 was used to detect SNPs, DNPs, and TNPs (single-base, double-base, and triple-base variations), and Strelka (v1.0.14) was used to detect InDels (insertions and deletions). (The above two softwares were used for variation detection here, but the algorithm is not limited to the detection results of the above two softwares.) VEP (v83) was then used to annotate the variation results. Before the final manual review of the annotation results, the candidate variations were filtered according to the clinical significance of the mutation using the following filtering conditions:

[0094] 1. Filter variants outside the WES (whole exon) region;

[0095] 2. Filter out synonymous mutations;

[0096] 3. Compare the non-hotspot variants in the COSMIC database (https: / / cancer.sanger.ac.uk / cosmic / ) and filter out variants with a frequency of less than 5%;

[0097] 4. For hotspot variants included in the COSMIC database, variants with a frequency of less than 1% were filtered out;

[0098] 5. Filter out variants with a site variation frequency greater than 1% in the commonly used 1000-gene database (suspected common variants).

[0099] Among them, the variant sites with a site variation frequency greater than 1% in the commonly used 1000-gene database include:

[0100] Variant sites with a population frequency greater than 1% in the NHLBI-ESP database (https: / / evs.gs.washington.edu / EVS / );

[0101] Variant sites with a population frequency greater than 1% in the 1000 Genomes Project database (http: / / phase3browser.1000genomes.org / index.html);

[0102] Variants with a frequency greater than 1% in East Asian populations from the 1000 Genomes Project data (http: / / phase3browser.1000genomes.org / index.html).

[0103] 6. Filter out low coverage sites, specifically sites with coverage less than 40X in tumor or less than 10X in normal;

[0104] 7. Filter out variants with absolute support numbers less than or equal to 6 reads.

[0105] After the above filtering steps, a large number of low-frequency variants that are difficult to identify and variants with unclear clinical significance will be filtered out.

[0106] The true labels of the candidate variants were determined by two researchers with two or more years of relevant bioinformatics experience (those skilled in the art). They standardized and standardized the IGV evaluation rules by studying the IGV SOP [Standard operating procedure for somatic variant refinement of sequencing data with paired tumor and normal samples DOI:10.1038 / s41436-018-0278-z]. The candidate variants were then evaluated for true and false labels using IGV. After the evaluation, sites where the two researchers disagreed were removed, and the remaining sites were used as gold standard sites for training the supervised learning model. After this step, a total of 11,223 variants were identified, including 9,181 SNPs and 2,042 indels.

[0107] True somatic variants should have sufficient molecular support (reads) (referring to molecules with unique identifiers after amplification, alignment, and deduplication), generally requiring ≥ 7 reads. Common erroneous variants may occur in the following situations:

[0108] (1) Suspected germline mutation. When the same mutation type is detected in normal control cells (usually white blood cells) and tumor tissues at the same time, and the possibility of tumor cells contaminating normal tissue cells can be ruled out with a high probability (the mutation frequency of the mutation in normal tissue cells is also high), it is highly likely that the mutation is a germline mutation.

[0109] (2) No read coverage or low coverage depth in normal tissue cells. When a candidate variant site has no sequencing coverage or a coverage depth of less than 10 layers in normal cells, it is impossible to determine whether the variant is supported in normal tissue cells. Therefore, the site in this case will generally be classified as a false positive.

[0110] (3) When there is no read coverage or the coverage depth is low (less than 40 layers) in cancer tissue cells, the site is generally classified as a false positive.

[0111] (4) Reads have a high degree of strand bias (e.g., >90% of reads are positive or negative strands); or (and) a certain percentage of reads contain misaligned bases. In this case, the site label is more difficult to determine.

[0112] (5) There are many variant bases on the supported reads, but the variants do not exist in the corresponding normal tissue cells. This situation is often caused by alignment errors or errors in the reference genome sequence.

[0113] (6) Multiple variant types appear simultaneously and each accounts for a certain proportion.

[0114] (7) Mutation frequency is too low. Generally, mutation frequencies below 5% are questionable, but this threshold is highly correlated with sequencing depth and varies from experiment to experiment. The average sequencing depth of tumors in this study was approximately 200X, and was uniformly filtered according to the criteria of 5% and mutation hotspots of 1%.

[0115] (8) There are inconsistent and sporadic base variations in each supported read.

[0116] (9) The alignment quality is too low. For example, more than 60% of the reads have an alignment quality less than 30.

[0117] (10) The base quality is too low. For example, the average base quality of supporting reads is less than 20.

[0118] (11) There are indels in the adjacent positions (upstream and downstream of 50 bp) of the supported reads. This situation increases the probability of alignment errors.

[0119] (12) Edge preference: that is, in all supporting reads, the variant sites are located at the same end of the reads.

[0120] (13) Low complexity of adjacent sequences.

[0121] 2. Filtering reads

[0122] This filtering step uses pysam (v0.11) of Python 2.7.8 to filter out reads in paired samples of candidate variant sites. The filtering conditions include the following:

[0123] (1) Read length less than 55;

[0124] (2) Read alignment quality is less than 10;

[0125] (3) Soft truncation length is greater than 20 bp.

[0126] Through the above steps, unreliable reads can be filtered out to prevent the introduction of unnecessary noise into the input training set.

[0127] 3. Extracting and Vectorizing the Features of Variant Sites

[0128] The following relevant features were extracted and vectorized for each variant site. Specifically, based on the reads after the above filtering steps and combined with the site interpretation considerations commonly used in manual review, the relevant features of a single site and the 50bp upstream and downstream sites were calculated. The relevant features and calculation methods are as follows:

[0129] Mutation support count (mut_Count): The number of reads supporting the mutation site in tumor tissue.

[0130] Mutation frequency (mut_Freq): The frequency of the mutation site to be detected in the tumor tissue (the number of mutation reads supported divided by the total depth of the site).

[0131] Mapping quality (mapping_Quality): The mean mapping quality of mutation-supporting reads in tumor tissue.

[0132] Base quality (base_Quality): The average quality of the bases at the mutation sites in the supporting mutation reads in tumor tissue.

[0133] Environmental quality (BQ_average): The average base quality of the 50 bp upstream and downstream sequences of the mutation site to be detected in the reads supporting the mutation site to be detected in the tumor tissue.

[0134] Mismatch_Average: The average mismatch rate of the 50 bp upstream and downstream sequences of the mutation site in the supporting mutation reads in tumor tissue.

[0135] HDR value (HDR_score): The score of suspected homology alignment errors in supporting mutation reads in tumor tissue.

[0136] Homology alignment errors often result in multiple variants with the same or similar mutation frequencies at other base positions upstream and downstream of the read supporting the mutation, which can easily lead to false positives at that site. This value is calculated only when the number of reads supporting the mutation is 5 or more; otherwise, the HDR value is 0.

[0137] Based on this error characteristic, we can first count the number of upstream and downstream position mutations, and then calculate the corresponding HDR value. When the absolute value of the difference between the mutation frequency of a certain upstream and downstream mutation in normal control tissue and tumor tissue is greater than or equal to 0.7, the mutation position is considered to meet the HDR condition. The HDR calculation method is:

[0138]

[0139] Among them, I i is equal to 1, and n is the number of sites that meet the HDR conditions.

[0140] When the absolute value of the difference between the mutation frequencies of a certain upstream and downstream mutation in normal control tissue and tumor tissue is less than 0.7, the mutation position is considered to be not in compliance with HDR and the HDR is 0.

[0141] The smaller the HDR value, the higher the probability of alignment errors for suspected homology.

[0142] SideB value (sideBias): The edge bias score of supporting mutation reads in tumor tissue, calculated as:

[0143] When the number of reads supporting the mutation is greater than or equal to 10:

[0144]

[0145] Where n = 50, Dt i is the depth of supporting reads for the i-th bp upstream of the mutation site; Dn i is the depth of supporting reads at the ith bp downstream of the mutation site. abs represents the absolute value.

[0146] When the number of reads supporting the mutation is less than 10, SideB is 0.

[0147] The lower the SideB score, the more marginal preference the site has.

[0148] Insert size (insertSize): The mean insert size of reads supporting the mutation site to be detected in tumor tissue;

[0149] StrandB value (strandBias): The strand bias score of mutation-supporting reads in tumor tissue is calculated by counting the proportion R of positive strands in these reads and calculating the strand bias score.

[0150] For variant sites with 10 or more reads supporting the mutation:

[0151]

[0152] Where R is the ratio of positive strands in reads supporting the mutation. abs represents the absolute value.

[0153] When the number of reads supporting the mutation is less than 10, StrandB is 0.

[0154] The higher the StrandB score, the less strand-biased the mutation is.

[0155] Genome complexity score (repeative_Flag): The complexity score of the reference genome sequence 20 bp upstream and downstream of the mutation position.

[0156] If 6 or more consecutive single-base repeat sequences or tandem repeat sequences (for example, 'AAAAAA' or 'ATCATCATCATCATCATCATC') appear within 20 bp upstream and downstream of the mutation site to be detected on the reference genome sequence, the feature value is 1, otherwise it is 0; the continuous single-base repeat sequence means that the same base appears repeatedly continuously, for example, 'AAAAAA'; in the continuous tandem repeat sequence, the tandem sequence refers to a sequence of two or more bases in series, for example, 'AT' or 'ATC', and the continuous tandem repeat sequence is, for example, 'ATCATCATCATCATCATC'.

[0157] Base ratio (ntRatio): The ratio of reads supporting mutations in normal control tissue and tumor tissue.

[0158] Other variant type ratio (var_TypeRatio): If there are other variant types in the tumor tissue besides a certain variant type, the ratio of the number of supporting reads of the other variant type to the number of supporting reads of the certain variant type.

[0159] Normal coverage depth (normal_Coverage): the coverage depth of the mutation site in normal control tissue.

[0160] Alignment length (query_Length): The average length of reads supporting mutations in tumor tissues.

[0161] Normal Indel Presence Rate (normal_Indels): Whether an InDel is present in the 50 bp upstream and downstream of the mutation site in tumor and / or normal tissues. If so, the product of the InDel's variation frequency and length is calculated, or these products are added together. If not, the value is 0. If an InDel is present in the 50 bp upstream and downstream of the mutation site in either tumor or normal tissue, the product of the InDel's variation frequency and length is calculated. If an InDel is present in the 50 bp upstream and downstream of the mutation site in both tumor and normal tissues, the product of the InDel's variation frequency and length is calculated and added together.

[0162] InDel length (indel_Lengths): The length of the InDel deletion or insertion.

[0163] 4. Model Training, Evaluation, and Prediction

[0164] Because base substitutions and InDel mutation characteristics have significant differences, such as the length and complexity of the InDel sequence itself having a significant impact on the authenticity of the site itself, and the InDel error rate is higher in tandem repeat regions, base substitutions and InDels have different reference standards and thresholds for measuring site authenticity. Therefore, we divided the base substitution and InDel training sets into two independent data sets and trained the models separately.

[0165] 1. Model building process

[0166] The data from the relevant features extracted in step 3 were divided into training and validation sets. Specifically, 6181 cases and 1430 variant sites were randomly selected from base substitutions (including SNVs, DNVs, and TNVs) and indels, respectively, as the training set, and the remaining sites were used as the test set, as shown in Table 1.

[0167] Table 1. Partitioning of training and test sets for model construction

[0168] Model training set Test set base substitution 6181 3000 Indel 1430 612

[0169] Next, base substitution and indel model construction and optimization were performed on the 6181 and 1430 training sets, respectively, and model performance was evaluated on the test set. Model construction was performed using the randomForest_4.16-14 package in R3.5.1.

[0170] (1) The features used in the construction of the base substitution model are shown in Table 2. The features in Table 2 are incremented in number from top to bottom (for example, when the model takes 3 features, the selected features are "mismatch_Average", "ntRatio", and "mutCount"), and the error rate of the model under each gradient is evaluated using the 5-fold cross-validation method on the training set (using the function "randomForest::rfcv"). The results are as follows Figure 2 As shown in the middle left figure, the horizontal axis is the number of features, and the vertical axis is the cross-validation error rate. It can be seen that for the base substitution model, when the number of features increases to 8, the error rate trend curve tends to be flat.

[0171] The models generated under the above different numbers of features were used to predict the test set, evaluate the area under the curve (AUC), and calculate the accuracy of the model. The number of gradient 1 and 2 features was too small and was not included in the analysis. The results are shown in Table 3.

[0172] Table 2 Characteristics of base substitution model construction

[0173] Serial number Feature Name 1 misMatch_Average 2 ntRatio 3 mut_Count 4 HDR_score 5 Base_Quality 6 strandBias 7 sideBias 8 Mut_Freq 9 Mapping_Quality 10 BQ_average 11 insertSize 12 normal_Coverage 13 normal_Indels 14 repeative_Flag

[0174] Table 3 Model effects with different numbers of features

[0175] Number of features Test set AUC Test set accuracy 3 0.981±0.005 0.939 4 0.982±0.004 0.939 5 0.989±0.003 0.955 6 0.992±0.003 0.96 7 0.992±0.003 0.96 8 0.992±0.003 0.961 <![CDATA[ 9 ]]> <![CDATA[ 0.993±0.003 ]]> <![CDATA[ 0.963 ]]> 10 0.994±0.003 0.964 11 0.992±0.003 0.965 12 0.992±0.003 0.967 13 0.995±0.002 0.968 14 0.995±0.002 0.967

[0176] (2) The features used in the InDel model construction are shown in Table 4. The features in Table 4 are quantitatively increased from top to bottom, and the error rate of the model under each gradient is evaluated using the 5-fold cross validation method on the training set (using the function "randomForest::rfcv"). The results are as follows Figure 2 As shown in the middle right figure, the horizontal axis is the number of features, and the vertical axis is the cross-validation error rate. It can be seen that for the InDel model, when the number of features increases to 10, the error rate trend curve tends to be flat.

[0177] The models generated under the above different numbers of features were used to predict the test set, evaluate the area under the curve (AUC), and calculate the accuracy of the model. The number of gradient 1 and 2 features was too small and was not included in the analysis. The results are shown in Table 5.

[0178] Table 4 Characteristics of InDel model construction

[0179]

[0180]

[0181] Table 5 Model effects with different numbers of features

[0182] Number of features Test set AUC Test set accuracy 3 0.986±0.009 0.943 4 0.991±0.007 0.948 5 0.991±0.007 0.951 6 0.992±0.007 0.948 7 0.992±0.007 0.953 8 0.992±0.007 0.956 9 0.992±0.007 0.951 <![CDATA[ 10 ]]> <![CDATA[ 0.994±0.006 ]]> <![CDATA[ 0.962 ]]> 11 0.993±0.007 0.964 12 0.993±0.006 0.966 13 0.994±0.006 0.971 14 0.994±0.006 0.969 15 0.994±0.006 0.969 16 0.994±0.006 0.972

[0183] The key model parameters, such as the number of trees (ntree) and the number of randomly extracted features (mtry), were optimized using a 10-fold cross-validation method on the training set. The optimal parameters after optimization are detailed in Table 6.

[0184] Table 6 Optimal parameters for SNV and Indel models

[0185] Model Best ntree Best mtry base substitution 900 5 Indel 2000 4

[0186] 2. Model performance on the test set

[0187] The model trained in the above steps was used to make predictions and evaluate the model on the independent test set.

[0188] The site features to be reviewed are extracted and vectorized according to step 3, and then predicted using the prediction model. Negative or positive is judged based on the threshold. For example, when the threshold is 0.5, a value below the threshold (0.5) is negative, and a value above or equal to the threshold (0.5) is positive.

[0189] in accordance with Figure 2Based on the cross-validation error rate curve in

[15] , the AUC metric on the test set, and the significance of the reference features themselves, a model with a base substitution feature count of 9 and an indel feature count gradient of 10 was selected. As shown in Tables 3 and 5, the AUCs were 0.993 ± 0.003 and 0.994 ± 0.006, respectively. The accuracy of the model on the independent test set at each gradient is shown in Tables 4 and 6. Accuracy = (number of predicted true positives + number of predicted true negatives) / total number of samples.

[0190] When the threshold is 0.5, the sensitivity and specificity are shown in Table 7. The AUC values ​​of the base substitution and InDel models are 0.993±0.002 and 0.994±0.006, respectively. Figure 3 The box plots of the test set prediction scores are shown in Figure 4 . Figure 4 As shown in the figure, the two models have good effects on distinguishing negative and positive sites in the independent test set.

[0191] Table 7 Performance indicators of base substitution and Indel models on the test set

[0192]

[0193] Note: The true positives and true negatives marked in the table are evaluated using the IGV method described previously.

[0194] The above embodiments are only used to illustrate the present invention, wherein the model feature calculation and algorithm tuning are subject to change. Any equivalent transformations and improvements based on the technical solution of the present invention should not be excluded from the scope of protection of the present invention.

Claims

1. A method for reviewing high-throughput sequencing gene variation detection results, comprising the following steps: (A) Construct a training set, including sequencing data of several positive variant sites and negative non-variant sites; (B) Extracting and vectorizing features of the variant site from the training set; the features of the variant site include: Mutation support number: the number of reads in tumor tissue that support the mutation site to be detected; Mutation frequency: the frequency of the mutation site to be detected in tumor tissue; Base quality: the average quality of the bases at the mutation site to be detected in the reads supporting the mutation site to be detected in the tumor tissue; Misalignment rate: The average misalignment rate of the 50 bp upstream and downstream sequences of the mutation site to be detected in the reads supporting the mutation site to be detected in tumor tissue; HDR value: the score of suspected homology alignment errors in reads supporting the mutation site to be detected in tumor tissue; SideB value: edge preference score of reads supporting the mutation site to be detected in tumor tissue; StrandB value: the strand preference score of reads supporting the mutation site to be detected in tumor tissue; Base ratio: the ratio of the number of reads supporting the mutation site to be detected in normal control tissue and tumor tissue; Other variant type ratio: If there are other variant types in the tumor tissue besides a certain variant type, the ratio of the number of supporting reads of the other variant types to the certain variant type; Alignment length: the average length of reads supporting the mutation site to be detected in tumor tissue; InDel length: the length of the deletion or insertion of Indel; (C) Using the feature results obtained in step (B), a random forest method is used to construct a model, and then the model is used to determine whether the variant at the test site is a variant site; The calculation method of the HDR value is: When the number of reads supporting the mutation site to be detected is less than 5, or when the number of reads supporting the mutation site to be detected is greater than or equal to 5 and the absolute value of the difference in mutation frequency between a certain variant upstream and downstream of the reads supporting the mutation site to be detected in normal control tissue and tumor tissue is less than 0.7, the HDR value is 0; When the number of reads supporting the mutation site to be detected is greater than or equal to 5, and the absolute value of the difference between the mutation frequencies of a certain variant upstream and downstream of the reads supporting the mutation site to be detected in normal control tissue and tumor tissue is greater than or equal to 0.7, the HDR value is calculated as follows: in, is equal to 1; n is the number of sites that meet the HDR conditions; The SideB value is calculated as follows: When the number of reads supporting the mutation site to be detected is greater than or equal to 10, the SideB value is calculated according to the following formula: ； Where n=50, Dt i is the depth of supporting reads for the i-th bp upstream of the mutation site; Dn i is the depth of supporting reads at the i-th bp downstream of the mutation site; When the number of reads supporting the mutation site to be detected is less than 10, the SideB value is 0; The StrandB value is calculated as follows: When the number of reads supporting the mutation site to be detected is greater than or equal to 10, the StrandB value is calculated according to the following formula: ； Where R is the ratio of the positive strand in the reads supporting the mutation site to be detected; When the number of reads supporting the mutation site to be detected is less than 10, the StrandB value is 0; The method is a non-disease diagnosis and treatment method.

2. The method according to claim 1, characterized in that The characteristics of the variant site also include any one or more of the following: Alignment quality: the average alignment quality of reads supporting the mutation site to be detected in tumor tissue; Environmental quality: The average base quality of the 50 bp sequences upstream and downstream of the mutation site to be detected in the reads supporting the mutation site to be detected in the tumor tissue; Insert size: the mean insert size of reads supporting the mutation site to be detected in tumor tissue; Genome complexity score: The complexity score of the reference genome sequence 20 bp upstream and downstream of the mutation site to be detected; Normal coverage depth: the coverage depth of the mutation site to be detected in normal control tissue; Normal InDel Presence Rate: Whether there are InDels 50 bp upstream and downstream of the mutation site to be detected in tumor tissues and / or normal tissues; if so, the product of the InDel variation frequency and length is calculated, or the products are further added; if not, the value is 0.

3. The method according to claim 1, wherein: The type of gene variation is base substitution or InDel; the model is constructed according to the variation type, and a training set is constructed using sequencing data of variation sites of the same variation type, and features are extracted and vectorized to construct a model.

4. The method according to claim 2, wherein: The genomic complexity score is calculated as follows: if 6 or more consecutive single-base repeat sequences or tandem repeat sequences appear within 20 bp upstream and downstream of the mutation site to be detected on the reference genomic sequence, the genomic complexity score is 1, otherwise it is 0.

5. The method according to any one of claims 1 to 4, characterized in that: Before extracting and vectorizing features of variant sites from the training set in step (B), the method further includes filtering reads in the training set that support the mutation sites to be detected as follows: filtering out reads with at least one of the following three conditions: read length less than 55, read alignment quality less than 10, and soft truncation length greater than 20 bp.

6. A system for reviewing gene variation detection results, comprising device A, device B, and device C; The device A is capable of constructing a training set of sequencing data including a plurality of positive variant sites and negative non-variant sites according to step (A) of any one of claims 1 to 5; The device B is capable of extracting and quantifying the features of the variant site from the training set according to step (B) of any one of claims 1 to 5; The device C can extract and vectorize the variant site features using the device B according to step (C) in any one of claims 1-5, use the random forest method to build a model, and then use the model to determine whether the variant site to be tested is a variant site.

7. Any of the following applications: (I) Use of the system according to claim 6 in the preparation of products for early tumor screening, tumor prognosis, tumor classification and / or tumor medication guidance; (II) Use of the system of claim 6 in detecting gene mutations; said application is a non-disease diagnosis and treatment application.

Images