Long-read data robust variant calling method and system based on data quality improvement

Abstract

The application discloses a long-read data robust variant detection method and system based on data quality improvement, and belongs to the technical field of gene sequencing and biological information analysis. The method of the application is aimed at the characteristics of long-read data, and first, preprocessing is completed through sequence filtering and base correction, and a three-dimensional outlier feature system of sequencing support intensity, mapping quality reliability and strand-specific bias is constructed. Based on the median absolute deviation algorithm, the features are standardized and outlier scores are generated, a double-threshold dynamic optimization mechanism is constructed, and high-confidence data is screened. The variant candidate regions are identified through multi-tool joint comparison, the SNP and Indel are typed through a deep learning model, and the reliability of the results is ensured through triple cross-validation. The method does not need to manually set a fixed threshold, effectively suppresses interference factors, balances detection accuracy, robustness and scene universality, and provides high-quality variant data support for clinical diagnosis and genomics research.

Description

Technical Field

[0001] This invention belongs to the field of gene sequencing and bioinformatics analysis technology, specifically involving a robust variant detection method and system for long-read data based on data quality improvement. It is applicable to the accurate detection of single nucleotide polymorphisms (SNPs) and insertion / deletion variants (Indels) in long-read sequencing platforms such as PacBio and ONT, and can be widely used in clinical diagnosis, genomics research, species evolution analysis and other scenarios. Background Technology

[0002] Long-read sequencing technology, with its advantages of ultra-long sequencing fragments and high genome coverage, can effectively traverse repetitive and complex genomic regions, demonstrating significant advantages in the detection of SNPs, Indels, and other variant sites. It has become one of the core technologies for genomics research and clinical molecular diagnostics. Currently, mainstream long-read sequencing platforms include the PacBioHiFi / CLR platform and the Oxford Nanopore Technologies (ONT) platform, whose generated sequencing data provide crucial data support for complex genome analysis and rare variant screening.

[0003] The core of variant detection is to accurately identify sites that differ from the reference genome from massive sequencing data. However, long-read sequencing data is susceptible to interference factors such as base identification errors, read mapping biases, and strand distribution biases, leading to problems such as high false-positive variant rates and missed detections of true variant sites. To address these issues, existing technologies primarily improve detection accuracy by constructing quality control systems, optimizing alignment algorithms, and introducing statistical models for correction. Among these, strategies based on outlier screening and dynamic threshold adjustment have become key directions for data quality optimization.

[0004] Existing methods for detecting variants in long-read data often employ fixed thresholds or single-feature evaluation models for outlier screening, lacking a comprehensive consideration of the multi-dimensional characteristics of sequencing data. This makes them difficult to adapt to the complex scenarios of different sequencing platforms and genomes of different species. Dynamic threshold adjustment strategies often rely on single statistical indicators without combining data quality and stability verification, resulting in insufficient threshold adaptability and the potential for over-filtering of valid data or residual interference data. Furthermore, existing methods suffer from insufficient precision in resolving breakpoints and inadequate cross-validation mechanisms during variant candidate region identification and genotyping, further impacting the robustness and reliability of variant detection.

[0005] In practical applications such as clinical diagnosis and precision medicine, the accuracy of variant detection directly determines the reliability of diagnostic results. False positive variants may lead to misdiagnosis, while missed detections may delay treatment. In genomics research, low-quality variant data can interfere with the accuracy of gene function analysis and evolutionary relationship inference. Therefore, how to construct a multi-dimensional feature fusion outlier screening system, design an adaptive dynamic threshold adjustment mechanism, and combine precise alignment and multiple validation strategies to achieve efficient and accurate detection of variant sites in long-read data, while also considering the universality and robustness of the method, has become a key technical problem that urgently needs to be solved in the field of long-read sequencing data processing. This has significant theoretical value and practical application prospects. Summary of the Invention

[0006] This invention addresses the aforementioned problems by providing a robust mutation detection method and system for long-read data based on data quality improvement. This method constructs a three-dimensional outlier feature system and a robust statistically driven scoring mechanism, avoiding the misjudgment of outliers caused by single-feature evaluation in traditional methods, and eliminating the need for manually setting fixed thresholds. By combining dynamic threshold optimization and multiple cross-validation strategies, the method can accurately capture real mutation signals while effectively suppressing various interference factors, thus balancing the accuracy, robustness, and versatility of mutation detection.

[0007] To solve the above-mentioned technical problems, the present invention adopts the following technical solution:

[0008] This invention provides a robust mutation detection method for long-read data based on data quality improvement, comprising the following steps:

[0009] Step 1: Construct a long-read sequencing data system and perform basic preprocessing to determine core data elements and screening rules, and complete the initial data purification.

[0010] Step 2: Construct a three-dimensional outlier feature system, extract sequencing support strength features, mapping quality reliability features, and strand specificity bias features, and quantify them;

[0011] Step 3: Calculate outlier scores for reading segments based on robust statistical methods. Use the Median Absolute Deviation (MedAD) algorithm to achieve robust standardization of multi-dimensional features, avoiding the sensitivity of traditional standardization methods to extreme values. Integrate the standardized features of each dimension through an equal-weight fusion strategy to generate a comprehensive outlier score that can objectively reflect the degree of anomaly of the reading segment, providing a quantitative basis for outlier screening.

[0012] Step 4: Construct a dual-threshold dynamic optimization mechanism for sorting consistency verification. Design dual sorting indexes based on the original features and robust outlier scores respectively. Quantify the consistency of the two sortings using the Spearman rank correlation coefficient to determine the stability of data quality and the degree of abnormal interference. Adaptively adjust the hard threshold for outlier screening based on the verification results to form a dynamic threshold system that adapts to different data quality states.

[0013] Step 5: Based on the dual thresholds dynamically optimized in Step 4, perform high-quality data screening, using hard thresholds as the main approach and soft thresholds as a fallback to ensure the removal of extreme outliers. At the same time, combine the three core indicators of coverage depth, average base quality, and chain distribution balance to perform secondary verification of the screening results. The threshold parameters are dynamically corrected through a feedback adjustment mechanism to ensure that the screened data removes interference while retaining sufficient variant signals, ultimately obtaining a high-confidence sequencing dataset.

[0014] Step 6: A multi-tool joint alignment strategy is used to complete the accurate genome alignment. The alignment parameters are optimized for the characteristics of different long-read sequencing platforms. The accuracy of the alignment results is improved by combining secondary alignment and duplicate read filtering. The system identifies three types of candidate regions of variants: abnormal coverage depth, abnormal clustering of CIGAR strings, and split alignment. The credibility of the candidate regions is evaluated by weighted comprehensive scoring, and high-confidence candidate regions of variants are screened.

[0015] Step 7: Use a deep learning model to accurately detect and classify SNPs and Indels in high-confidence candidate regions. Introduce a multi-cross-validation mechanism including breakpoint back-substitution verification, homologous sequence alignment, and sample specificity verification to significantly reduce the false positive mutation rate. Finally, output the mutation data in the VCF standard format and generate a visualization report containing mutation statistics, chromosome distribution, and other information to support downstream analysis applications.

[0016] Furthermore, step 1 includes the following specific steps:

[0017] Step 1.1: Construct a long-read sequencing data system;

[0018] Obtain raw long-read sequencing data of the target species' genome. The data is from the PacBio HiFi / CLR or ONT sequencing platform, and the length of a single sequencing read is not less than 1000 bp to cover variant types such as SNPs and short fragment indels in the genome.

[0019] Step 1.2: Define core data elements and filtering rules;

[0020] ① Raw sequencing read set: contains nucleic acid sequence fragments generated by the sequencing platform. Each read is associated with four basic attributes: reference genome start mapping coordinates (in bp), reference genome end mapping coordinates (in bp), mapping confidence score, and CIGAR string.

[0021] ② Standard reference genome: Experimentally validated standardized genome sequences from NCBI RefSeq or Ensembl databases are used as benchmark templates for alignment and variant determination;

[0022] ③ Data filtering parameters: Set an effective sequence length threshold to ensure the integrity of sequence information for variant detection;

[0023] Step 1.3: Perform the basic preprocessing procedure;

[0024] Short reads and contaminating sequences below the effective sequence length threshold were removed using sequence filtering tools, and error correction algorithms were used to correct the sequencing data.

[0025] For ONT data: Use Guppy or Dorado neural network models to optimize base recognition accuracy;

[0026] For PacBio data: improve sequence consistency using the CCS iterative correction algorithm;

[0027] The final result is a preprocessed sequencing dataset.

[0028] Furthermore, step 2 includes the following specific steps:

[0029] Step 2.1: Extract sequencing support strength features;

[0030] For each preprocessed read, its effective alignment length within the target genome region is calculated using a formula to quantify the read's support for the target region. The calculation process is as follows:

[0031] (1)

[0032] in, The ability of the reading segment to support the target region, The starting coordinates of the read segment are mapped. The terminating coordinate is the coordinate of the map. When the terminating coordinate is less than the starting coordinate, the default effective support length is 1 to avoid calculation errors. The effective alignment lengths of all read segments are integrated to form a support strength set.

[0033] Step 2.2: Quantify the reliability characteristics of the mapping quality;

[0034] The mapping quality value of each read segment is extracted to construct a mapping quality set, which reflects the confidence level of the read segment alignment with the reference genome. At the same time, based on the base identification quality value of the read segment, the average confidence of a single base is calculated to construct a base quality set. Through the synergistic analysis of the mapping quality set and the base quality set, the reliability of the read segment mapping process is characterized.

[0035] Step 2.3, analyze chain-specific bias features;

[0036] By combining the relationship between the CIGAR string identifier and the read segment mapping coordinates, the DNA strand orientation of the read segment can be determined: when > When it is determined to be a positive chain, the eigenvalue is recorded as 1; when < When a read is identified as a reverse strand, its characteristic value is recorded as -1. The strand direction characteristic values ​​of all reads are integrated to form a strand bias characteristic set, which quantifies the strand distribution balance of the sequencing data.

[0037] Furthermore, step 3 includes the following specific steps:

[0038] Step 3.1, feature standardization based on median absolute deviation;

[0039] For the support strength set, mapping quality set, and chain bias feature set constructed in step 2, robust standardization is performed using the Median Absolute Deviation (MedAD) algorithm to avoid the sensitivity of the traditional Z-score method to extreme values.

[0040] For any feature set V, standardization is achieved through formulas (2)-(4):

[0041] (2)

[0042] (3)

[0043] (4)

[0044] Where Med(V) is the median of set V, MedAD(V) is the absolute deviation of the median of set V, and 0.6745 is the standardized correction coefficient under the normal distribution; standardized support strength sets are generated respectively. Standardized mapping quality set Standardized chain bias set ;

[0045] Step 3.2: Integrate multi-dimensional features to generate outlier scores;

[0046] Using an equal-weighted fusion strategy, the comprehensive outlier score for each reading segment is calculated, as shown in formula (5) below:

[0047] (5)

[0048] Among them, outlier in the formula s core i The outlier score for the i-th read segment. To standardize the i-th element in the set of support strengths, For the i-th element in the set of standardized mapping quality, Let be the i-th element in the standardized chain bias set; in the textual description, each absolute value term is the absolute value of the corresponding standardized value; integrate all outlier scores to form an outlier score set. ; Calculate the maximum value of the set Compared with the average As a global statistical indicator.

[0049] Furthermore, step 4 includes the following specific steps:

[0050] Step 4.1, calculate the initial double thresholds;

[0051] The outlier score set is sorted in ascending order to obtain the sorted outlier scores. The initial threshold is calculated using formulas (6) and (7).

[0052] (6)

[0053] (7)

[0054] Where S is the sorted set of outlier scores. Med(S) is the median of the sorted set of outlier scores, MedAD(S) is the absolute deviation of the median of the sorted set of outlier scores, and percentile(S, 75) is the 75th percentile of the sorted set of outlier scores.

[0055] Step 4.2: Perform a sorting consistency check;

[0056] Constructing a dual-sorting index: the original feature sorting index ( That is, the sorted index of the original support strength set after traditional Z-score standardization) and the robust score sorted index ( (i.e., a sorted index based on outlier scores), the Spearman rank correlation coefficient is calculated using formula (8). :

[0057] (8)

[0058] in, Let be the rank difference between the i-th read segment and the two rankings, and n be the total number of read segments; the closer the Spearman rank correlation coefficient is to 1, the higher the ranking consistency.

[0059] Step 4.3, adaptively adjust the threshold;

[0060] According to Spearman's rank correlation coefficient The result was an adjustment of the hard threshold. :

[0061] If the Spearman rank correlation coefficient is <0.3, the hard threshold is adjusted as follows:

[0062] =Med(S)+4×MedAD(S)(9)

[0063] If 0.3 ≤ Spearman's rank correlation coefficient < 0.7, keep the initial hard threshold unchanged;

[0064] If the Spearman rank correlation coefficient is ≥0.7, the hard threshold is adjusted as follows:

[0065] =Med(S)+2.5×MedAD(S)(10).

[0066] Furthermore, step 5 includes the following specific steps:

[0067] Step 5.1: Perform double threshold filtering;

[0068] The optimized dual threshold is applied to filter the preprocessed sequencing dataset: reads that satisfy the outlier score of the i-th read ≤ hard threshold are retained; if the proportion of valid reads is less than 30%, a soft threshold is used as an alternative screening threshold.

[0069] Step 5.2: Verify and correct the filtering results;

[0070] Calculate the coverage depth distribution, average base quality value, and chain distribution balance of the filtered dataset;

[0071] If the coverage depth is less than 10 times, the adjusted hard threshold = current hard threshold - 0.5 × MedAD(S), and the filter is re-filtered according to the adjusted threshold;

[0072] If the average base quality is less than 30 times, the adjusted hard threshold = current hard threshold + 0.3 × MedAD(S), and the screening is re-selected according to the adjusted threshold;

[0073] This ultimately results in a high-quality sequencing dataset.

[0074] Furthermore, step 6 includes the following specific steps:

[0075] Step 6.1: Perform a multi-tool joint comparison;

[0076] The reference genome was aligned using the mapping alignment tool (minimap2): PacBio data was aligned using the -ax map-pb parameter, and ONT data was aligned using the -ax map-ont parameter. The BWA-MEM algorithm was used for secondary alignment. The results were integrated, duplicate reads were marked, and low-quality alignments were filtered out using the sequence alignment mapping tool (SAMtools) to generate a binary sequence alignment mapping format file (BAM format file).

[0077] Step 6.2: Identify candidate mutation regions;

[0078] Based on the binary sequence alignment mapping format file (BAM file), three types of candidate regions are identified: coverage depth abnormal regions, alignment operation string (CIGAR string) abnormal clustering regions, and split alignment regions. The starting coordinates of the candidate region, the ending coordinates of the candidate region, the number of supported read segments, and the mutation type label of each region are recorded to construct a set of mutation candidate regions.

[0079] Step 6.3: Evaluate the quality of candidate regions;

[0080] Calculate the overall quality score for each candidate region. The score consists of the average mapping quality, average base quality, and read consistency of the reads within the region. Set a threshold and retain regions with scores ≥ the threshold to obtain a set of high-confidence variant candidates.

[0081] Furthermore, step 7 includes the following specific steps:

[0082] Step 7.1, Mutation type detection

[0083] A deep learning model was used to genotype the high-confidence candidate variant set: for SNP sites, the genotype was accurately determined by analyzing the base distribution pattern of all supporting reads at the candidate site and combining the degree of conservation of the site in homologous sequences; for Indel variants, the specific sequence fragments and lengths of insertions or deletions were reconstructed by accurately analyzing the alignment breakpoints and characteristics of reads with the reference genome, thus completing the variant genotyping.

[0084] Step 7.2, cross-validate the mutation results

[0085] The localization error was ensured to be ≤ 2bp through back-substitution verification at breakpoints; the rationality of the variation was evaluated by comparison with the genome aggregation database, Mendelian inheritance laws were verified through family data, and tumor-normal tissue pairing verification was performed through tumor data.

[0086] Step 7.3, Output the standardized results

[0087] Output the final variant detection results according to the variant call format standard (VCF v4.3 format). The results should include the full quality control labels generated by this method, including key information such as outlier filtering status of reads, comprehensive quality score of candidate regions, and variant validation pass rate. Simultaneously generate a visual analysis report in portable document format (PDF) or hypertext markup language format (HTML). The report should cover a statistical summary of variant types (SNP, Indel), a chromosome-level variant density distribution map, a detailed list of high-confidence variant sites, and a summary of core quality indicators, providing a clear basis for downstream functional annotation and experimental validation.

[0088] In another aspect, the present invention provides a robust variant detection system for long-read data based on data quality improvement, for performing the aforementioned method. The system specifically includes a computer processor and memory, a long-read sequencing data acquisition unit, an offline variant detection analysis and processing unit, and an online variant result output matching unit.

[0089] The long-read sequencing data acquisition unit uses a computer processor to perform length screening, contaminated sequence removal, and base correction preprocessing on the long-read raw sequencing data (from PacBio HiFi / CLR or Oxford Nanopore (ONT) sequencing platforms) collected during the genome sequencing of the target species, thereby constructing a long-read sequencing training dataset and loading it into computer memory.

[0090] The offline analysis and processing unit for mutation detection, based on the long-read sequencing training dataset generated by the long-read sequencing data acquisition unit, divides the three-dimensional outlier feature vectors in the training dataset—sequencing support strength, mapping quality reliability, and strand specificity bias—into several candidate mutation regions according to the credibility of the mutation signal and the degree of data anomaly.

[0091] The online mutation result output matching unit accurately matches the real-time sequencing data of the target species' genome with the variant candidate region categories obtained by the offline mutation detection analysis and processing unit. Combined with deep learning genotyping and multiple cross-validation results, the final SNP and Indel mutation detection results are output to the genome analysis management platform, which can provide effective evidence and support for various scenarios such as genome functional annotation, genetic variation association analysis, disease molecular mechanism research, and clinical diagnostic support.

[0092] Compared with the prior art, the present invention has the following advantages:

[0093] This invention addresses the core pain points of variant detection in long-read sequencing data by constructing a three-dimensional outlier feature system to replace the traditional single-feature evaluation model. This system comprehensively covers factors influencing data interference, avoiding misjudgment of outliers caused by single features, and is flexibly adaptable to complex scenarios involving different sequencing platforms such as PacBio and ONT, as well as different species genomes. This method employs data stability verification to achieve dynamic threshold adaptive optimization, eliminating the need for manually setting fixed thresholds and effectively solving the problems of insufficient adaptability and excessive data filtering or residual interference associated with traditional dynamic thresholds. Furthermore, this invention integrates multi-tool joint alignment and multiple cross-validation strategies, accurately capturing true variant signals while significantly suppressing interference factors such as base identification errors and mapping biases. This invention improves variant detection performance across the entire chain from data preprocessing to result output, balancing detection accuracy, robustness, and scenario versatility, providing high-quality variant data support for clinical diagnosis and genomics research. Attached Figure Description

[0094] To more clearly illustrate the technical solutions in the embodiments of this application, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0095] Figure 1 This is a schematic diagram of the method flow of the present invention;

[0096] Figure 2 This is a system structure diagram of the present invention. Detailed Implementation

[0097] To gain a deeper understanding of this invention, we will provide a comprehensive and detailed description. However, this invention has various implementations and is not limited to the specific examples listed herein. These examples are presented to enhance a full understanding of the disclosure of this invention.

[0098] This embodiment uses human genome (GRCh38 version) sequencing data generated by the PacBio HiFi sequencing platform as an example. Specifically, the dataset is long-read sequencing data with a sequencing coverage of 30-fold, a single read length range of 1000-20000 bp, and an average read length of 15000 bp. The reference genome is the human genome GRCh38 version indexed in the NCBI RefSeq database, containing 24 chromosomes (chromosomes 1-22, X chromosome, and Y chromosome) and mitochondrial genome sequences. This embodiment focuses on the detection of SNPs and Indel variations on human chromosome 1 to verify the effectiveness and advancement of the method of this invention.

[0099] like Figure 1 and Figure 2 As shown, a robust mutation detection method for long-read data based on data quality improvement includes the following steps:

[0100] Step 1.1, Constructing a long-read sequencing data system

[0101] Nucleic acid sequence fragments and associated attribute information are extracted from the raw sequencing data to construct a raw sequencing read set. Each read is associated with four basic attributes: reference genome start coordinates (pos), reference genome end coordinates (...), and reference genome termination coordinates (...). ), mapping confidence score ( The alignment operation string (CIGAR string) is used. The reference genome GRCh38 chromosome 1 starts at coordinate 0 and ends at coordinate 248956422.

[0102] Step 1.2: Define core data elements and filtering rules

[0103] ① Data filtering parameters: Set the minimum sequence length ( ≥1000bp, mapping confidence score ( A threshold of ≥20 is used as the basic screening threshold to remove short reads and low-confidence alignment reads, ensuring the integrity of sequence information and the reliability of alignment.

[0104] ② Standard reference genome: The GRCh38 version of the human genome sequence included in the NCBI RefSeq database was used, with the chromosome 1 sequence as the core alignment template. The gene annotation file (VCF format) of this version was also downloaded for subsequent functional annotation and verification of variant sites.

[0105] ③ Base quality standards: A single base identification quality value (Q value) ≥20 is set as a qualified base, and the proportion of qualified bases in a single read segment is ≥80%. Read segments below this standard are marked as low-quality read segments and are preferentially removed in subsequent preprocessing.

[0106] Step 1.3, execute the basic preprocessing procedure.

[0107] ① Sequence filtering: The raw sequencing data was filtered using the FASTP tool with the parameter set to "-- 1000 -- -- 20”, remove reads shorter than 1000bp, with single base quality values ​​below 20, and contaminated reads, and retain valid reads after filtering;

[0108] ② Base Correction: Based on the characteristics of PacBio HiFi data, the Circular Consistent Sequence (CCS) iterative correction algorithm is used for base correction, with the number of iterations set to 3. After correction, the average single base accuracy of the read segment is improved from 99.2% to 99.9%. At the same time, the Guppy neural network model (version 6.4.2) is used to optimize the base recognition results and further reduce the base recognition error.

[0109] ③ Data standardization: The filtered and corrected reads are stored in BAM format to construct a preprocessed reads set. The mapping coordinates, CIGAR strings and other attributes of the reads are uniformly formatted to ensure that the data format meets the requirements of subsequent analysis.

[0110] Step 2: Construct and quantify a three-dimensional outlier feature system. This step involves constructing a three-dimensional outlier feature system covering data support capability, alignment reliability, and strand distribution balance for the preprocessed sequencing data. Sequencing support strength features, mapping quality reliability features, and strand specificity bias features are extracted and quantified. Specifically, this includes the following steps:

[0111] Step 2.1, Extract sequencing support strength features

[0112] For each preprocessed read, its effective alignment length in the corresponding region of chromosome 1 is calculated using formula (1), quantifying the read's coverage support capability for the target region:

[0113] (1)

[0114] in, The starting coordinates of the read segment are mapped. The endpoint mapping coordinates of the read segment are used. When the endpoint mapping coordinates of the read segment are greater than the starting mapping coordinates, the effective support length is the difference between the two. When the endpoint mapping coordinates of the read segment are less than the starting mapping coordinates (in the reverse alignment scenario), the default effective support length is 1 to avoid calculation errors. In this embodiment, the effective support length of the read segment ranges from 1 to 18600 bp, with an average effective support length of 12800 bp. The effective alignment lengths of all read segments are integrated to form a support strength set (…). ).

[0115] Step 2.2, Quantify the mapping quality reliability characteristics

[0116] ① Extract the mapping confidence score for each read segment ( The value ranges from 0 to 60. A higher value indicates a higher confidence level in the alignment of the read with the reference genome. This is used to construct a mapping quality set. In this embodiment, The average value was 38.6 and the median was 40, indicating that the confidence level of most reads was good.

[0117] ② Based on the single-base identification quality value (Q value) of the read segment, calculate the average confidence level of each read segment using the formula "average confidence level of a single base = (sum of Q values ​​of all bases) / read segment length", and construct a base quality set. In this embodiment, The average value is 32.1, which meets the preset base quality standard;

[0118] ③ Through and The synergistic analysis of the two was conducted, and the Pearson correlation coefficient was used to measure the correlation between them. In this embodiment, the correlation coefficient was 0.72, which indicates that the read segment alignment confidence and the base recognition quality have a strong positive correlation and can be used together as a quantitative indicator of the reliability of the mapping quality.

[0119] Step 2.3, analyze chain-specific bias features

[0120] By combining the CIGAR string identifier with the read segment mapping coordinates, the DNA strand orientation of the read segment is jointly determined: when the termination mapping coordinate of the read segment is greater than the start mapping coordinate, it is identified as the positive strand (+ strand), and the feature value is recorded as 1; when the termination mapping coordinate of the read segment is less than the start mapping coordinate, it is identified as the negative strand (- strand), and the feature value is recorded as -1. The strand orientation feature values ​​of all read segments are integrated to form a strand bias feature set (…). The distribution balance of the sequencing data was quantified. In this embodiment, the proportion of sense reads was 51.2%, and the proportion of antisense reads was 48.8%, indicating a basically balanced distribution with no obvious chain bias.

[0121] Step 3: Calculate outlier scores based on robust statistical methods. This step uses the Median Absolute Deviation (MedAD) algorithm to robustly standardize multi-dimensional features, and then fuses the standardized features to generate a comprehensive outlier score, providing a quantitative basis for outlier screening. Specifically, this includes the following steps:

[0122] Step 3.1, Feature standardization based on median absolute deviation

[0123] For the construction in step 2 , and The MedAD algorithm is used for robust standardization to avoid the sensitivity of the traditional Z-score method to extreme values. For any feature set V, standardization is achieved through formulas (2)-(4):

[0124] (2)

[0125] (3)

[0126] (4)

[0127] Where Med(V) is the median of set V, MedAD(V) is the absolute deviation of the median of set V, and 0.6745 is the standardized correction coefficient under normal distribution. Let be the i-th element in set V. In this embodiment, the standardization results of each feature set are as follows: The standardized mean is 0.02 and the standard deviation is 0.98. The standardized mean is -0.01 and the standard deviation is 1.03. After standardization, the mean is 0.05 and the standard deviation is 0.97. The distributions of the standardized features tend to be consistent, allowing for fusion analysis. Standardized support strength sets are generated separately. Standardized mapping quality set Combined, standardized chain bias set .

[0128] Step 3.2: Integrate multi-dimensional features to generate outlier scores.

[0129] An equal-weighted fusion strategy is adopted to balance the influence of features of each dimension on outlier determination, and the comprehensive outlier score of each reading segment is calculated as follows: Formula (5) is as follows:

[0130] (5)

[0131] Among them, outlier s core i Let | be the outlier score of the i-th read segment. |、| |、| | represents the absolute value of the corresponding standardized numerical value; a higher score indicates a higher degree of anomaly in the read segment. Outlier scores from all read segments are integrated to form an outlier score set ( ), calculate the maximum value of the set ( The value is 4.82, and the average value is ( The average value was 0.86, and the median was 0.73, providing a reference for subsequent threshold setting.

[0132] Step 4: Construct a dual-threshold dynamic optimization mechanism. This step designs a dual-sorting index, verifies data stability using the Spearman rank correlation coefficient, adaptively adjusts the hard threshold for outlier filtering, and forms a dynamic threshold system that adapts to the current data quality status. Specifically, this includes the following steps:

[0133] Step 4.1, Calculate the initial double threshold.

[0134] right Sort in ascending order to obtain the outlier scores. The initial threshold is calculated using formulas (6) and (7):

[0135] (6)

[0136] (7)

[0137] Where S is the sorted set of outlier scores ( Med(S) is The median (0.73 in this example) is MedAD(S). The median absolute deviation (0.31 in this example) is percentile(S, 75). The 75th percentile (1.25 in this example). Substituting this into the calculation yields the initial hard threshold ( ) = 0.73 + 3 × 0.31 = 1.66, soft threshold ( =1.25.

[0138] Step 4.2, perform sorting consistency check.

[0139] Constructing a dual-sorting index: ① Original feature sorting index ( ): Set the original support strength ( ① After traditional Z-score standardization, the index is sorted in descending order of scores; ② Robust score-sorted index ( Based on outlier scores ( The scores are sorted in descending order. The Spearman rank correlation coefficient is calculated using formula (8). ), quantifying the consistency between the two types of sorting:

[0140] (8)

[0141] in, Let be the rank difference between the i-th read segment and the two sortings, and n be the total number of read segments (823,000 in this example). Calculations show that in this example... =0.78, which is close to 1, indicating that the two sorting methods have high consistency, the current data quality is stable, and the degree of abnormal interference is low.

[0142] Step 4.3, Adaptively adjust the threshold

[0143] according to The hard threshold is adjusted accordingly, following these rules:

[0144] like If the value is less than 0.3, the hard threshold is adjusted as follows:

[0145] =Med(S)+4×MedAD(S)(9)

[0146] If 0.3≤ If the threshold is less than 0.7, keep the initial hard threshold unchanged.

[0147] like ≥ 0.7, hard threshold adjusted to:

[0148] =Med(S)+2.5×MedAD(S)(10)

[0149] In this embodiment =0.78≥0.7, therefore the adjusted hard threshold =0.73+2.5×0.31=1.505, rounded to two decimal places as 1.51; the soft threshold remains unchanged at 1.25, thus forming a dynamic threshold system: hard threshold 1.51, soft threshold 1.25.

[0150] Step 5: Perform high-quality data screening based on dynamic thresholds. This step uses the adjusted dynamic thresholds as the core, performs double threshold filtering, and combines three core indicators to perform secondary verification and feedback correction on the screening results to obtain a high-confidence sequencing dataset. Specifically, it includes the following steps:

[0151] Step 5.1, perform double threshold filtering

[0152] The optimized dual threshold was applied to filter the preprocessed reads set: retaining those that met the criteria. Read segments with a hard threshold (1.51) or lower are selected; simultaneously, the percentage of valid read segments is calculated. If the percentage is less than 30%, a soft threshold (1.25) is used as an alternative screening threshold. In this embodiment, The number of read segments with a value ≤1.51 is 768,000, and the effective read segment ratio is 93.3%, which is higher than 30%. Therefore, there is no need to replace the threshold, and 768,000 read segments are obtained after filtering.

[0153] Step 5.2: Verify and correct the filtering results.

[0154] The three core metrics of the filtered dataset are calculated: coverage depth distribution, average base quality value, and chain distribution balance. The threshold parameters are dynamically adjusted based on the metric results, and the specific rules are as follows:

[0155] ① Coverage depth verification: Calculate the average coverage depth of each region of chromosome 1. In this embodiment, the average coverage depth is 27.8×, which is higher than the standard of 10×, so no threshold adjustment is required; if the coverage depth is lower than 10×, the current hard threshold is reduced by 0.5×MedAD(S) (0.155 in this embodiment), that is, the adjusted hard threshold = 1.51-0.155=1.355, and the screening is re-executed;

[0156] ② Base quality verification: Calculate the average base quality value of the reads after screening. In this example, it is 33.5, which is higher than the standard of 30×, so no threshold adjustment is required. If the average base quality is lower than 30×, the current hard threshold is increased by 0.3×MedAD(S) (0.093 in this example), that is, the adjusted hard threshold = 1.51 + 0.093 = 1.603, and the screening is performed again.

[0157] ③ Chain distribution balance verification: After screening, the proportion of positive chain reads is 51.5%, and the proportion of negative chain reads is 48.5%, which is a balanced distribution with no obvious bias.

[0158] After verification, the screening results of this embodiment meet all the indicator requirements, and no threshold adjustment is needed. The final result is a high-quality sequencing dataset containing 768,000 reads, which provides reliable data for subsequent alignment and variant detection.

[0159] Step 6: Multi-tool joint comparison and variant candidate region identification

[0160] This step employs a multi-tool joint alignment strategy to complete accurate genome alignment, systematically identifies three types of candidate variant regions and performs quality assessment, and screens high-confidence candidate variant regions. Specifically, it includes the following steps:

[0161] Step 6.1: Perform a multi-tool joint comparison.

[0162] ① Primary alignment: The minimap2 tool was used to align the high-quality sequencing dataset with the GRCh38 version chromosome 1 sequence. Based on the characteristics of PacBio HiFi data, the -ax map-pb parameter was selected, the number of alignment threads was set to 16, and the alignment results in SAM format were generated.

[0163] ② Secondary alignment: After converting the SAM format results to BAM format, the BWA-MEM algorithm is used for secondary alignment with the parameters set to "-t 16 -M" to further optimize the alignment position of the read segments and improve the alignment accuracy.

[0164] ③ Result optimization: The secondary alignment results were processed using SAMtools (version 1.17), including sorting, deduplication, and index building. Low-quality alignment reads with a mapping confidence score below 20 were filtered out, and finally, an optimized BAM format file was generated, with an alignment accuracy of 99.5%.

[0165] Step 6.2, Identify candidate mutation regions

[0166] Based on the optimized BAM file, the system identifies three types of candidate mutation regions, records the core information of each region, and constructs a set of candidate mutation regions. ):

[0167] ① Coverage depth abnormal areas: The coverage depth threshold is set to "average coverage depth ± 2 × standard deviation" (27.8 ± 5.6 in this embodiment, i.e. 22.2 × ~ 33.4 ×). Areas with coverage depth lower than 22.2 × or higher than 33.4 × are identified as abnormal areas. A total of 128 such areas were identified.

[0168] ② Abnormal clustering regions of CIGAR strings: For read segments of CIGAR strings containing "insertion (I)" and "deletion (D)" operations and with a length ≥5bp, their clustering regions were counted, and a total of 356 such regions were identified, which are speculated to be potential Indel mutation regions;

[0169] ③ Segmentation of alignment regions: Identify alignment regions where a single read segment is segmented to different positions on chromosome 1. A total of 89 such regions were identified, which are speculated to be regions associated with potential structural variations.

[0170] For each candidate region, record its starting coordinates ( ), Termination coordinates ( Supported number of segments to read () ) and variant type tags ( ), build It contains a total of 573 candidate regions for mutation.

[0171] Step 6.3, Evaluate the quality of candidate regions

[0172] The comprehensive quality score for each candidate region was calculated using a weighted summation method, with the following weights: average reading mapping quality (weight 0.4), average base quality (weight 0.3), and reading consistency (weight 0.3). The score range was 0-100 points, with a threshold of 25 points. Regions with scores ≥25 points were retained as the high-confidence variant candidate set. ).

[0173] In this embodiment, a total of 542 regions with a comprehensive quality score ≥25 were identified, including 119 regions with abnormal coverage depth, 338 regions with abnormal CIGAR string aggregation, and 85 regions with splitting and comparison features. This provides accurate candidate regions for subsequent variant detection and typing.

[0174] Step 7: Precise mutation detection, cross-validation, and result output

[0175] This step uses a deep learning model to accurately classify high-confidence candidate regions using SNP and Indel, introduces a multiple cross-validation mechanism to reduce the false positive rate, and outputs standardized variant data and a visualization report. Specifically, it includes the following steps:

[0176] Step 7.1, Detection of variant types

[0177] The DeepVariant model (version 1.5.0) was used for... Perform mutation typing, and set the model parameters to "-- PACBIO -- 16”, employing a differential analysis strategy for different types of variation:

[0178] ① SNP genotyping: The base distribution pattern of all supporting reads at candidate sites was analyzed, and the sequence conservation characteristics of chromosome 1 in GRCh38 version (obtained through UCSC Genome Browser) were combined to accurately determine the genotype (homozygous or heterozygous). A total of 2863 SNP sites were detected.

[0179] ② Indel typing: Accurately analyzes the alignment breakpoints and characteristics of reads with the reference genome, reconstructs the specific sequence fragments and lengths of insertions or deletions, and detects a total of 789 Indel variants, including 412 insertion variants and 377 deletion variants, with lengths ranging from 1 to 50 bp.

[0180] Step 7.2, cross-validate the mutation results

[0181] A triple cross-validation mechanism is introduced to verify the reliability of mutation results at each level, thereby reducing the false positive mutation rate.

[0182] ① Breakpoint back-substitution verification: The detected variant sites and breakpoint information are back-substituted into the original alignment data to re-verify the alignment consistency and ensure that the positioning error is ≤2bp; In this embodiment, the positioning error of 99.2% of the variant sites is ≤2bp, which meets the standard;

[0183] ② Database comparison and verification: The mutation results were compared with the human chromosome 1 mutation data in the gnoMedAD database (version 3.1) to evaluate the rationality of the mutations; in this embodiment, 2789 SNP sites and 763 Indel mutations could be matched in the database, with a matching rate of 94.3%, indicating that the mutation results are highly reliable;

[0184] ③ Sample specificity verification: The peripheral blood sample used in this embodiment has no known genetic diseases, and the verification results show no pathogenic variant sites, which is consistent with the sample characteristics; if it is a family sample, it is necessary to verify whether the variant sites conform to Mendel's laws of inheritance; if it is a tumor sample, it is necessary to screen for tumor tissue-specific somatic variant sites.

[0185] Step 7.3, Output the standardized results

[0186] ① Standardized data output: The final variant detection results are output in VCF v4.3 format, which includes complete information on all 3652 validated variant sites (2863 SNPs and 789 Indels), with quality control labels including outlier filtering status of reads, comprehensive quality score of candidate regions, and variant validation pass rate.

[0187] ② Visualization Report Generation: Generates a PDF visualization report covering: statistical charts of variant types (SNPs account for 78.4%, Indels account for 21.6%), a distribution map of chromosome 1 variant density (variables are mainly concentrated in gene-rich regions), a detailed list of high-confidence variant sites (including coordinates, genotype, and validation status), and a summary of core quality indicators (effective read rate 93.3%, variant validation pass rate 94.3%, and localization error ≤2bp proportion 99.2%), providing clear evidence for downstream gene function annotation and clinical analysis.

[0188] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. For example, this embodiment uses PacBio HiFi data, but it can also be applied to ONT sequencing data.

[0189] A robust variant detection system for long-read data based on data quality improvement includes a computer processor and memory, a long-read sequencing data acquisition unit, an offline variant detection analysis and processing unit, and an online variant result output and matching unit.

[0190] The long-read sequencing data acquisition unit uses a computer processor to perform length screening, contaminated sequence removal, and base correction preprocessing on the long-read raw sequencing data (from PacBio HiFi / CLR or Oxford Nanopore (ONT) sequencing platforms) collected during the genome sequencing of the target species. This preprocessing constructs a long-read sequencing training dataset, which is then loaded into computer memory. The offline variant detection and analysis unit, based on the long-read sequencing training dataset generated by the long-read sequencing data acquisition unit, categorizes the three-dimensional outlier feature vectors (sequencing support strength, mapping quality reliability, and strand specificity bias) in the training dataset into several candidate variant regions according to the confidence level of the variant signal and the degree of data anomaly. The online variant result output matching unit accurately matches the real-time genome sequencing data of the target species with the candidate variant regions obtained by the offline variant detection and analysis unit. Combining deep learning genotyping and multiple cross-validation results, the final SNP and Indel variant detection results are output to the genome analysis management platform. This provides effective evidence and support for various scenarios, including genome functional annotation, genetic variant association analysis, disease molecular mechanism research, and clinical diagnostic support.

[0191] This embodiment focuses on human chromosome 1, but can also be extended to other species' genomes and whole-genome variation detection scenarios. Therefore, the embodiment should be considered exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the claims are intended to be included within the invention.

[0192] Contents not described in detail in this specification are prior art known to those skilled in the art. Although illustrative specific embodiments of the invention have been described above to facilitate understanding by those skilled in the art, it should be understood that the invention is not limited to the scope of the specific embodiments. Various modifications are readily apparent to those skilled in the art as long as they fall within the spirit and scope of the invention as defined and determined by the appended claims, and all inventions utilizing the concept of this invention are protected.

Claims

1. A robust mutation detection method for long-read data based on data quality improvement, characterized in that, Includes the following steps: Step 1: Construct a long-read sequencing data system and perform basic preprocessing to determine core data elements and screening rules, and complete the initial data purification. Step 2: Construct a three-dimensional outlier feature system, extract sequencing support strength features, mapping quality reliability features, and strand specificity bias features, and quantify them; Step 3: Calculate outlier scores for reading segments based on robust statistical methods. Use the median absolute deviation algorithm to achieve robust standardization of multi-dimensional features, avoiding the sensitivity of traditional standardization methods to extreme values. Integrate the standardized features of each dimension through an equal-weight fusion strategy to generate a comprehensive outlier score that can objectively reflect the degree of anomaly of the reading segment, providing a quantitative basis for outlier screening. Step 4: Construct a dual-threshold dynamic optimization mechanism for sorting consistency verification. Design dual sorting indexes based on the original features and robust outlier scores respectively. Quantify the consistency of the two sortings using the Spearman rank correlation coefficient to determine the stability of data quality and the degree of abnormal interference. Adaptively adjust the hard threshold for outlier screening based on the verification results to form a dynamic threshold system that adapts to different data quality states. Step 5: Based on the dual thresholds dynamically optimized in Step 4, perform high-quality data screening, using hard thresholds as the main approach and soft thresholds as a fallback to ensure the removal of extreme outliers. At the same time, combine the three core indicators of coverage depth, average base quality, and chain distribution balance to perform secondary verification of the screening results. The threshold parameters are dynamically corrected through a feedback adjustment mechanism to ensure that the screened data removes interference while retaining sufficient variant signals, ultimately obtaining a high-confidence sequencing dataset. Step 6: A multi-tool joint alignment strategy is used to complete the accurate genome alignment. The alignment parameters are optimized for the characteristics of different long-read sequencing platforms. The accuracy of the alignment results is improved by combining secondary alignment and duplicate read filtering. The system identifies three types of candidate regions of variants: abnormal coverage depth, abnormal clustering of CIGAR strings, and split alignment. The credibility of the candidate regions is evaluated by weighted comprehensive scoring, and high-confidence candidate regions of variants are screened. Step 7: Use a deep learning model to accurately detect and classify SNPs and Indels in high-confidence candidate regions. Introduce a multi-cross-validation mechanism including breakpoint back-substitution verification, homologous sequence alignment, and sample specificity verification to significantly reduce the false positive mutation rate. Finally, output the mutation data in the VCF standard format and generate a visualization report containing mutation statistics, chromosome distribution, and other information to support downstream analysis applications.

2. The robust mutation detection method for long-read data based on data quality improvement according to claim 1, characterized in that, Step 1 includes the following specific steps: Step 1.1: Construct a long-read sequencing data system; Obtain raw long-read sequencing data of the target species' genome. The data is from the PacBio HiFi / CLR or ONT sequencing platform, and the length of a single sequencing read is not less than 1000 bp. Step 1.2: Define core data elements and filtering rules; ① Raw sequencing read set: contains nucleic acid sequence fragments generated by the sequencing platform. Each read is associated with four basic attributes: reference genome start mapping coordinates, reference genome end mapping coordinates, mapping confidence score, and CIGAR string. ② Standard reference genome: Experimentally validated standardized genome sequences from NCBI RefSeq or Ensembl databases are used as benchmark templates for alignment and variant determination; ③ Data filtering parameters: Set an effective sequence length threshold to ensure the integrity of sequence information for variant detection; Step 1.3: Perform the basic preprocessing procedure; Short reads and contaminating sequences below the effective sequence length threshold were removed using sequence filtering tools, and error correction algorithms were used to correct the sequencing data. For ONT data: Use Guppy or Dorado neural network models to optimize base recognition accuracy; For PacBio data: improve sequence consistency using the CCS iterative correction algorithm; The final result is a preprocessed sequencing dataset.

3. The robust mutation detection method for long-read data based on data quality improvement according to claim 2, characterized in that, Step 2 includes the following specific steps: Step 2.1: Extract sequencing support strength features; For each preprocessed read, its effective alignment length within the target genome region is calculated using a formula to quantify the read's support for the target region. The calculation process is as follows: (1) in, The ability of the reading segment to support the target region, The starting coordinates of the read segment are mapped. The terminating coordinate is the coordinate of the map. When the terminating coordinate is less than the starting coordinate, the default effective support length is 1 to avoid calculation errors. The effective alignment lengths of all read segments are integrated to form a support strength set. Step 2.2: Quantify the reliability characteristics of the mapping quality; Extract the mapping quality value of each read segment and construct a mapping quality set. Simultaneously, based on the base identification quality value of the read segment, calculate the average confidence of a single base and construct a base quality set. Through the collaborative analysis of the mapping quality set and the base quality set, characterize the reliability of the read segment mapping process. Step 2.3, analyze chain-specific bias features; By combining the relationship between the CIGAR string identifier and the read segment mapping coordinates, the DNA strand orientation of the read segment can be determined: when > When it is determined to be a positive chain, the eigenvalue is recorded as 1; when < When a read is identified as a reverse strand, its characteristic value is recorded as -1. The strand direction characteristic values ​​of all reads are integrated to form a strand bias characteristic set, which quantifies the strand distribution balance of the sequencing data.

4. The robust mutation detection method for long-read data based on data quality improvement according to claim 3, characterized in that, Step 3 includes the following specific steps: Step 3.1, feature standardization based on median absolute deviation; For the support strength set, mapping quality set, and chain bias feature set constructed in step 2, robust standardization is performed using the median absolute deviation algorithm to avoid the sensitivity of the traditional Z-score method to extreme values. For any feature set V, standardization is achieved through formulas (2)-(4): (2) (3) (4) Where Med(V) is the median of set V, MedAD(V) is the absolute deviation of the median of set V, and 0.6745 is the standardized correction coefficient under the normal distribution; standardized support strength sets are generated respectively. Standardized mapping quality set Combined, standardized chain bias set ; Step 3.2: Integrate multi-dimensional features to generate outlier scores; Using an equal-weighted fusion strategy, the comprehensive outlier score for each reading segment is calculated, as shown in formula (5) below: (5) Among them, outlier in the formula s core i The outlier score for the i-th read segment. To standardize the i-th element in the set of support strengths, For the i-th element in the set of standardized mapping quality, Let be the i-th element in the standardized chain bias set; in the textual description, each absolute value term is the absolute value of the corresponding standardized value; integrate all outlier scores to form an outlier score set. ; Calculate the maximum value of the set Compared with the average As a global statistical indicator.

5. A robust mutation detection method for long-read data based on data quality improvement according to claim 4, characterized in that, Step 4 includes the following specific steps: Step 4.1, calculate the initial double thresholds; The outlier score set is sorted in ascending order to obtain the sorted outlier scores. The initial threshold is calculated using formulas (6) and (7). (6) (7) Where S is the sorted set of outlier scores. Med(S) is the median of the sorted set of outlier scores, MedAD(S) is the absolute deviation of the median of the sorted set of outlier scores, and percentile(S, 75) is the 75th percentile of the sorted set of outlier scores. Step 4.2: Perform a sorting consistency check; Construct a dual-ranking index: the original feature-ranking index and the robust score-ranking index. Calculate the Spearman rank correlation coefficient using formula (8). : (8) in, Let be the rank difference between the i-th read segment and the two rankings, and n be the total number of read segments; the closer the Spearman rank correlation coefficient is to 1, the higher the ranking consistency. Step 4.3, adaptively adjust the threshold; According to Spearman's rank correlation coefficient The result was an adjustment of the hard threshold. : If the Spearman rank correlation coefficient is <0.3, the hard threshold is adjusted as follows: =With(S)+4×WithAD(S)(9) If 0.3 ≤ Spearman's rank correlation coefficient < 0.7, keep the initial hard threshold unchanged; If the Spearman rank correlation coefficient is ≥0.7, the hard threshold is adjusted as follows: =Med(S)+2.5×MedAD(S)(10).

6. The robust mutation detection method for long-read data based on data quality improvement according to claim 5, characterized in that, Step 5 includes the following specific steps: Step 5.1: Perform double threshold filtering; The optimized dual threshold is applied to filter the preprocessed sequencing dataset: reads that satisfy the outlier score of the i-th read ≤ the hard threshold are retained; If the percentage of valid read segments is less than 30%, a soft threshold will be used as an alternative screening threshold. Step 5.2: Verify and correct the filtering results; Calculate the coverage depth distribution, average base quality value, and chain distribution balance of the filtered dataset; If the coverage depth is less than 10 times, the adjusted hard threshold = current hard threshold - 0.5 × MedAD(S), and the filter is re-filtered according to the adjusted threshold; If the average base quality is less than 30 times, the adjusted hard threshold = current hard threshold + 0.3 × MedAD(S), and the screening is re-selected according to the adjusted threshold; This ultimately results in a high-quality sequencing dataset.

7. A robust mutation detection method for long-read data based on data quality improvement according to claim 6, characterized in that, Step 6 includes the following specific steps: Step 6.1: Perform a multi-tool joint comparison; Reference genome alignment was performed using a mapping alignment tool: PacBio data was aligned with the -ax map-pb parameter and ONT data was aligned with the -ax map-ont parameter. The BWA-MEM algorithm was used for secondary alignment. The results were integrated, duplicate reads were marked, and low-quality alignments were filtered out by the sequence alignment mapping tool to generate a binary sequence alignment mapping format file. Step 6.2: Identify candidate mutation regions; Based on binary sequence alignment mapping format files, three types of candidate regions are identified: regions with abnormal coverage depth, regions with abnormal clustering of alignment operation strings, and regions with split alignment. The starting coordinates of the candidate region, the ending coordinates of the candidate region, the number of supported reading segments, and the mutation type label of each region are recorded to construct a set of mutation candidate regions. Step 6.3: Evaluate the quality of candidate regions; Calculate the overall quality score for each candidate region. The score consists of the average mapping quality, average base quality, and read consistency of the reads within the region. Set a threshold and retain regions with scores ≥ the threshold to obtain a set of high-confidence variant candidates.

8. A robust mutation detection method for long-read data based on data quality improvement according to claim 7, characterized in that, Step 7 includes the following specific steps: Step 7.1, mutation type detection; A deep learning model was used to genotype the high-confidence candidate variant set: for SNP sites, the genotype was accurately determined by analyzing the base distribution pattern of all supporting reads at the candidate site and combining the degree of conservation of the site in homologous sequences; for Indel variants, the specific sequence fragments and lengths of insertions or deletions were reconstructed by accurately analyzing the alignment breakpoints and characteristics of reads with the reference genome, thus completing the variant genotyping. Step 7.2, cross-validate the mutation results; The localization error was ensured to be ≤ 2bp through back-substitution verification at breakpoints; the rationality of the variation was evaluated by comparison with the genome aggregation database, Mendelian inheritance laws were verified through family data, and tumor-normal tissue pairing verification was performed through tumor data. Step 7.3, output the standardized results; The final mutation detection results should be output according to the standard mutation call format. The results should include the full quality control labels generated by this method, including key information such as the outlier filtering status of the read segment, the comprehensive quality score of the candidate region, and the mutation verification pass rate. Simultaneously generate a visual analysis report in portable document format or hypertext markup language format. The report should include a statistical summary of variant types, a distribution map of chromosome-level variant density, a detailed list of high-confidence variant sites, and a summary of core quality indicators, providing a clear basis for downstream functional annotation and experimental verification.

9. A robust mutation detection system for long-read data based on data quality improvement, characterized in that: The system is used to perform the method according to any one of claims 1 to 8, and the system specifically includes a computer processor and memory, a long read sequencing data acquisition unit, an offline analysis and processing unit for variant detection, and an online output and matching unit for variant results. The long read sequencing data acquisition unit uses a computer processor to perform length screening, contaminated sequence removal, and base correction preprocessing on the long read raw sequencing data collected during the genome sequencing of the target species, thereby constructing a long read sequencing training dataset and loading it into computer memory. The offline analysis and processing unit for mutation detection divides the three-dimensional outlier feature vectors in the training dataset into several candidate regions based on the confidence level of the mutation signal and the degree of data anomaly, according to the long-read sequencing training dataset generated by the long-read sequencing data acquisition unit. The online mutation result output matching unit accurately matches the real-time sequencing data of the target species' genome with the variant candidate region categories obtained by the offline mutation detection analysis and processing unit. Combined with deep learning genotyping and multiple cross-validation results, the final SNP and Indel mutation detection results are output to the genome analysis management platform, which can provide effective evidence and support for various scenarios such as genome functional annotation, genetic variation association analysis, disease molecular mechanism research, and clinical diagnostic assistance.

Images