DNA sequence local alignment method based on k-bounded matching statistics

Abstract

The invention discloses a DNA sequence local alignment method based on k-bounded matching statistics, and aims to solve the problems of slow alignment, limitation on downstream application and more resource occupation caused by internal contradictions and limitations in speed, sensitivity, coordinate information retention and resource consumption when large-scale genome data is processed by the conventional DNA sequence alignment method. The method comprises the following steps of: constructing a composite index for a reference DNA sequence set by utilizing spectrogram Burrows-Wheeler conversion, scanning and querying basic groups of DNA sequences one by one by utilizing the composite index, generating a k-bounded matching statistical array, generating a de-randomized matching statistical array by utilizing a random matching model, defining a grammar rule and a translation rule, and carrying out translation on the de-randomized matching statistical array. And processing the de-randomized matching statistical array by using a grammar rule and a translation rule to obtain a character type comparison event of the query DNA sequence, and finding a comparison result corresponding to the query DNA sequence in the reference DNA sequence set based on the character type comparison event.

Description

Technical Field

[0001] This invention relates to the field of genome sequence alignment, specifically to a method for local DNA sequence alignment based on k-bounded matching statistics combined with an efficient k-mer index data structure. Background Technology

[0002] DNA sequence alignment is the cornerstone of modern computational biology and genomics. Its core task is to find high-quality local similarity regions for a given query sequence (such as a genome assembly fragment or sequencing read) within a large-scale reference sequence database. With the widespread adoption of Next-Generation Sequencing (NGS) technology, genomic data is growing at an unprecedented rate, placing extremely high demands on the computational efficiency and scalability of sequence alignment algorithms. Most current mainstream alignment tools, such as BLAST and BWA, employ a "seed-and-extend" strategy. This strategy first uses rapid "seed" matching to screen for candidate alignment regions, and then uses a more time-consuming "extend" algorithm (such as Smith-Waterman dynamic programming) to obtain the final accurate alignment. However, this classic paradigm faces inherent challenges when dealing with massive amounts of data:

[0003] (1) Trade-off between speed and sensitivity: Heuristic tools (such as BLAST) use longer k-mers (shorter words) as seeds in pursuit of speed, but this leads to a decrease in sensitivity when processing sequences with many variations or errors, which may miss important biological signals. However, using shorter k-mers to improve sensitivity leads to an exponential increase in the number of seed matches, which overwhelms subsequent processing and significantly increases runtime. This inherent contradiction forces existing tools to compromise between speed and discovery capability, making it difficult to simultaneously meet the requirements of high efficiency and high sensitivity in high-throughput analysis.

[0004] (2) Loss of coordinate information and limitations of alignment-free methods: In recent years, to address the challenges of large data scale, many alignment-free methods based on k-mer counting have emerged. These methods perform rapid similarity assessment by comparing the k-mer sets between sequences, avoiding time-consuming base-by-base alignment, and perform well in scenarios such as metagenomic classification and sequence clustering. However, their main drawback is that they completely discard the position (coordinate) information of k-mers in the sequence when processing them, thus failing to provide the precise positioning required for alignment and limiting their application in downstream analyses that require precise coordinates (such as variant detection and gene annotation).

[0005] (3) Computational resource bottleneck: Traditional alignment tools heavily rely on data transfer between the CPU and main memory. Due to the large size of genome sequences, the index structure often cannot be fully loaded into the cache, leading to frequent memory accesses and forming the so-called "memory wall" bottleneck. In addition, many tools require a large amount of temporary disk space during operation, which poses an obstacle in high-performance computing clusters and resource-constrained environments (such as personal computers or lightweight cloud instances), limiting their deployment flexibility and scalability. Extensive research has been conducted both domestically and internationally to address this challenge. Internationally, research directions include developing new index data structures (such as FM-index, B-tree, etc.) to optimize query speed, and designing dedicated hardware accelerators (such as FPGA, ASIC). However, most of these methods are still limited to optimizing a certain link in the existing "seed-extend" framework. Domestic research also focuses on algorithm optimization and parallel computing, but there is still room for exploration in proposing new algorithm theories that can fundamentally change the alignment paradigm. Summary of the Invention

[0006] To address the inherent contradictions and limitations in speed, sensitivity, coordinate information preservation, and resource consumption of existing DNA sequence alignment methods when processing large-scale genomic data, which lead to slow alignment, limited downstream applications, and high resource consumption, this invention proposes a local DNA sequence alignment method based on k-bounded matching statistics.

[0007] The technical solution adopted in this invention is:

[0008] It includes the following steps:

[0009] S1. Construct a composite index using the Burrows-Wheeler transformation of the spectral data to a set of reference DNA sequences;

[0010] S2. Use the composite index to scan the bases of the query DNA sequence one by one to generate a k-bounded matching statistics array. The k-bounded matching statistics array includes the matching length between the query DNA sequence and the reference DNA sequence set.

[0011] S3. Establish a random matching model, and based on the k-bounded matching statistics array, use the random matching model to generate a derandomized matching statistics array.

[0012] S4. Define syntax rules and translation rules, process the derandomized matching statistics array using syntax rules and translation rules to obtain character-type alignment events of the query DNA sequence, and find the alignment result corresponding to the query DNA sequence in the reference DNA sequence set based on the character-type alignment events.

[0013] The beneficial effects of this invention are as follows:

[0014] This invention employs a method combining k-bounded matching statistics with an efficient k-mer index data structure. The core idea is to transform discrete k-mer matching events into a continuous statistical signal rich in location information. Local sequence alignment of the query sequence is then performed using signal processing techniques such as random matching models, defined syntax rules, and translation rules. This approach retains the extremely high query speed of the k-mer method while effectively utilizing the sequence's location information. Therefore, it proposes an efficient, fast, and accurate local alignment method for large-scale sequence data, achieving high sensitivity and high accuracy in local alignment while maintaining high speed, high sensitivity, and low resource consumption.

[0015] This invention fundamentally solves the problem of speed versus sensitivity in existing technologies by transforming the discrete sequence matching problem into a continuous signal recognition problem. On the one hand, the core computation of this invention is a one-time linear scan of the query sequence, the time complexity of which is proportional to the length of the query sequence, ensuring stable and predictable high-speed performance. On the other hand, by repairing signal interruptions caused by sequence mutations through a derandomization algorithm, this method has natural robustness to errors and mutations in the sequence, and can ensure high sensitivity without sacrificing speed, thereby discovering alignments that traditional heuristic tools may miss.

[0016] This invention innovatively combines the computational efficiency of alignment-free methods with the coordinate-level accuracy of traditional alignment methods by intrinsically preserving the precise location information of each matching event in a k-bounded matching statistics array. Each element of the MS array strictly corresponds to a base position in the query sequence, ensuring lossless preservation of coordinate information in the initial computation stage. Subsequent derandomization and translation steps can accurately parse the boundary coordinates of local alignments from the signal, effectively solving the problems of existing alignment-free techniques failing to provide precise location information and traditional alignment techniques incurring excessive computational overhead.

[0017] This invention significantly improves deployment flexibility. The core algorithm of this invention executes entirely in main memory, without relying on any temporary disk files, completely eliminating disk I / O bottlenecks. Simultaneously, the SBWT index it relies on is a space-efficient compressed data structure, effectively reducing memory consumption. Due to its lightweight and disk-independent characteristics, the method of this invention can be compiled into WebAssembly and run directly in modern web browsers, greatly lowering the technical threshold for bioinformatics analysis and providing an ultimate solution for data privacy protection—a unique advantage unmatched by existing desktop or server-side comparison tools.

[0018] This invention represents a methodological shift from traditional combinatorial pattern matching to a modern signal processing paradigm, offering broad versatility and scalability. It redefines the sequence alignment problem as the task of detecting and interpreting non-random signals within a data stream filled with random noise. This pipeline framework of "computing the original signal -> denoising and enhancement -> signal translation and interpretation" is a highly general and modular design, and its core ideas can be directly applied to pattern discovery in other types of sequence data, such as protein sequence analysis. Attached Figure Description

[0019] Figure 1 This is an overall architecture diagram of the present invention;

[0020] Figure 2 This is a schematic diagram illustrating the translation results applied to the derandomized matching statistics array;

[0021] Figure 3 This is a flowchart of the search, mutation detection, and sequence alignment processes in this invention; Detailed Implementation

[0022] Specific implementation method one: Combining Figures 1-3 This embodiment describes a local DNA sequence alignment method based on k-bounded matching statistics, which includes the following steps:

[0023] S1. Construct a composite index using the Burrows-Wheeler transformation of the spectral data as a reference DNA sequence set. The specific process is as follows:

[0024] S11. Obtain a set of reference DNA sequences and set the length of each k-mer in the DNA sequence to k, where k is any positive integer, and the total number of k-mers is n. For bacterial genomes, k can be set to 31. For more complex mammalian genomes, a larger k value can be chosen. The choice of k value is a trade-off between the specificity and sensitivity of the alignment; specificity avoids random matching, while sensitivity allows for differences in sequences.

[0025] S12. Based on all k-mers in the reference DNA sequence set, construct an initial k-spectrum set R, which is the set of all different k-mers in the reference sequence. The SBWT data structure is theoretically a compact encoding of the de Bruijn graph, where nodes typically represent... The substring of length k is represented by an edge, and the edges represent k-mer nodes. To ensure the path integrity of this graph structure (i.e., to ensure that every non-virtual k-mer node in the graph has a backtracking path of at least length k), the initial k-spectrum set R needs to be filled to obtain a complete filled k-spectrum. The fill operation is as follows:

[0026] A subset is obtained by selecting from the initial k-spectral set R. subset All k-mers Length prefix It is not a suffix of any other k-mer in the initial k-spectral set R. In the de Bruijn graph, this corresponds to a node lacking an "ingress path," so it needs to be extended to the subset. Add several (Z) virtual k-mers prefixed with the special placeholder $ (a character whose lexicographical order is less than any base) to the subset, and then... Add to filter subset From the initial k-spectrum set R, a complete k-spectrum is obtained. .

[0027] S13, based on k-spectrum An SBWT index is obtained using the Burrows-Wheeler transform of the spectrogram. The SBWT index compactly stores the results of all k-mers sorted according to colexicographical order. Colexicographical order is a dictionary-based sorting method that compares characters from right to left, ensuring that all k-mers sharing the same suffix (e.g., all ending in GCA) are clustered in contiguous intervals after sorting. The SBWT index is a space-efficient compressed data structure; its core function is to... (Extending to the right) and (Left shortening) The operation simulates the extension and shortening of the sequence, achieving efficient sequence traversal.

[0028] Next, to enhance the SBWT index, an additional Longest Common Suffix (LCS) array corresponding to the k-mer co-order sorting in the SBWT index is calculated and stored. This array records the length of the longest common suffix between each k-mer and its preceding k-mer under the co-order sorting. The LCS array provides crucial information for the subsequent fast calculation of the MS array. This invention uses the SBWT index and the longest common suffix array together to form a composite index for subsequent fast queries. The composite index is a composite data structure that supports efficient k-mer queries.

[0029] S2. Using a composite index, scan the DNA sequence base by base to generate a k-bounded matching statistics array. The specific process is as follows:

[0030] S21. Obtain a query DNA sequence Q and create an integer array MS with the same length as the query DNA sequence Q. The array MS is used to store the matching length results.

[0031] S22. Initialize the matching length of the query DNA sequence. And an initial matching range representing all k-mers in the composite index. ,in, To locate a position on a DNA sequence, The initial value is the entire interval of all k-mer values. Let be the starting point for all k-mer intervals. This represents the end of all k-mer intervals.

[0032] S23. Using the composite index from S13, scan the bases of the query DNA sequence one by one to generate a k-bounded matching statistical array. This calculation process is a linear, pipelined operation, which can be viewed as sliding a window across the query DNA sequence and searching in real time for the longest match between the sequence within the window and the reference DNA sequence set. The specific process is as follows:

[0033] Perform an "expand-contract" matching search in each scan:

[0034] (a) Attempt to match extensions:

[0035] (1)

[0036] in, To query the first DNA sequence The bases at each position, To query the first DNA sequence The matching range corresponding to each position To query the first DNA sequence The purpose of this operation is to query the matching range corresponding to each position. The ending, length is Does the suffix exist in the reference DNA sequence?

[0037] (ii) Successful extension: If Returns a non-empty range, indicating that the match continues, i.e., with... The ending, length is The suffix exists in the reference DNA sequence, thus determining the current position. Corresponding matching range At this moment, the current position The match length is increased, and it is updated to However, it cannot exceed the preset k-mer upper limit k, i.e. .

[0038] (iii) Extension Failure: If Returns an empty range, representing the current match length. Unable to Continuing to extend, that is, at a matching length of Append to the matching suffix The corresponding term could not be found in the reference DNA sequence. The ending, length is The suffix does not exist in the reference DNA sequence, and the match is within this length. The above interrupt occurs. At this point, the algorithm does not start from zero, but enters a shrinking loop, which involves repeatedly calling... The operation intelligently shortens the left end of the previous matching range, searching for a shorter one that can still be matched. The successfully extended suffix indicates that this process is repeated cyclically. and decrease After each contraction, repeat the matching extension in step (a) until the extension attempt is successful or... This returns the current matching range and matching length.

[0039] Based on the above scanning and matching process, the obtained positions Match length Store in array The entire process ensures the integrity of the integer array. The record is based on The longest possible match at the end, i.e., the longest length value that can completely match a substring in the reference DNA sequence set, has an upper limit of 1. This step calculates the match length between each base position in the query DNA sequence and the reference DNA sequence in linear time complexity. After all bases in the query DNA sequence have been scanned and matched, the match lengths for all positions in the query DNA sequence are stored in the array MS, which is the complete k-bounded matching statistics array.

[0040] S3. Establish a random matching model to distinguish between random noise and true homologous signals. The random matching model is based on the following: under the null hypothesis (i.e., the query DNA sequence and the reference DNA sequence are unrelated), the length of a single random match follows a geometric distribution, while the longest random match length in the entire reference DNA sequence set follows an extreme value distribution. A dynamic significance threshold for the match length is calculated based on the random matching model. The establishment of the random matching model considers the probability of k-mers appearing randomly in a large genome, ensuring the scientific validity and adaptability of the threshold. The specific process is as follows:

[0041] S31. Due to random collisions, the MS array contains a large amount of noise signals generated by accidental matches. These signals are graphically represented as low, cluttered "troughs." This step aims to filter out the noise and highlight the true alignment signals, which are graphically represented as towering, continuous "mountains." Therefore, this invention proposes allowing users to define the statistical false positive rate tolerance. (For example (), used to control the balance between sensitivity and specificity of the present invention. Based on the total number of k-mers. and statistical false positive rate tolerance The significance threshold for matching length is determined by solving the following inequality. :

[0042] (2)

[0043] in, Represents the maximum matching length that is generated entirely randomly. Exceeding the set threshold The probability, It is a random variable representing the maximum randomly generated match length. Its cumulative distribution function, based on the extreme value distribution, is:

[0044] (3)

[0045] in, The probability of two nucleotides randomly matching (e.g., 0.25) reflects the background nucleotide composition of the DNA sequence.

[0046] Solving according to equations (2)-(3), we can obtain:

[0047] (4)

[0048] Significance threshold for match length This ensures that only matches that occur by chance with extremely low probability are considered potential signals.

[0049] S32. The MS array is derandomized and augmented using a significance threshold to eliminate noise and generate a smooth derandomized matching statistics (DMS) array. The specific process is as follows:

[0050] Create a derandomized matching statistics array DMS of the same length as the array MS, and initialize the last value of the array DMS. , Given the sequence length, a significance threshold based on the matching length is then applied. From the second to last position of the queried DNA sequence Iterate backwards from position 1 to the corresponding array MS:

[0051] (1) Signal recognition: If ,and This indicates an upward trend, that is ,express If it is itself a strong and continuous signal from the same source, then it is considered... It is reliable, so retain it and directly copy its value to obtain the array DMS. .

[0052] (2) Noise filling: If ,but If it is judged as noise, ,but If it is determined to be a matching plateau period, ,but It was determined to be a match breakpoint. When If it is determined to be noise, a matching plateau, or a matching breakpoint, then it is considered... Unreliable. In this case, the invention utilizes the inherent mathematical properties of the array MS ( ),from The next repaired location Linear extrapolation is performed to fill the current position, thereby ensuring the smoothness of the array DMS. .

[0053] This process employs a reverse iterative algorithm to smooth and correct the MS array generated by S2, effectively filtering out statistically insignificant matches below the stated significance threshold and repairing spurious signal breaks caused by single base differences. The core idea behind derandomizing the MS array is to preserve the ascending slopes of the MS array "mountain view" and outwardly fill in the "valleys" and "flat peaks," much like backfilling from a mountain peak to a valley. This repairs the signal "cliffs" caused by single variations or errors into smooth slopes, thereby connecting the alignment regions that should be continuous and outputting a DMS array with a higher signal-to-noise ratio that better reflects the sequence collinearity trend. This process is similar to filtering and smoothing in signal processing, aiming to restore the true matching signal interrupted by small variations or errors in the DNA sequence.

[0054] S4. The defined syntax rule decodes a continuous numerical array (DMS) into discrete character-based comparison events. This rule is a one-way traversal algorithm based on local context. arrive Iterate through the array DMS, at each position The algorithm will check a by , and The local numerical pattern is composed of [various elements]. Based on the predefined translation rules matched by this pattern, the position [is determined]. This will be decoded into corresponding character-type alignment events. Character-type alignment events are basic alignment events with clear biological significance, such as matches (M), mismatches or single-base insertions (X), single-base deletions or structural breakpoints (R), and regions specific to the query DNA sequence (I). By combining character-type alignment events at different positions on the query DNA sequence, a complete local alignment relationship between the query DNA sequence and the reference DNA sequence can be reconstructed. For example, a continuous increase in the value in the DMS array corresponds to a continuous base match, while a sudden drop or plateau in the value in the DMS array indicates the presence of a variation or alignment breakpoint. This conversion from numerical signals to symbolic sequences is a key bridge connecting statistics and biological meaning.

[0055] The translation rules are defined as a mapping table that associates local numerical patterns of different shapes in the DMS array with specific biological events. The mapping table must contain at least the following rules:

[0056] 1) For the "V"-shaped valley pattern of mismatch or single base insertion:

[0057] When the array DMS is querying DNA sequences The position presents a pattern such as positive numbers, zero or negative numbers, or 1. This pattern is interpreted as a single-point mutation causing a matching interruption and rapid recovery. Therefore, the position... Translated as "Event X".

[0058] 2) "Cliff" mode for single base deletions or structural breakpoints:

[0059] When the array DMS is querying DNA sequences The position shows values ​​that are much higher than the significance threshold. A sudden drop below the significance threshold The pattern, interpreted as a discontinuous break in a significantly homologous matching region, corresponds to structural differences (such as missing, long insertions, or rearrangements) between the query and reference sequences. Therefore, the position... and Translated as the left and right boundaries of the "R" event.

[0060] 3) "Plains" mode for querying DNA sequence-specific regions:

[0061] when If it is not part of the "V"-shaped valley pattern, it is interpreted as not having a corresponding match in the reference DNA sequence, and the position is... Translated as "I" event.

[0062] 4) For the matched "platform / ramp" pattern:

[0063] when Furthermore, if the above three patterns are not satisfied, it is interpreted as a continuous homologous match; therefore, the position is... Translated as the "M" event.

[0064] like Figure 2 As shown, the derandomized matching statistics array is processed sequentially using syntax rules and translation rules to obtain character-type alignment events in the query DNA sequence. Based on the character-type alignment events, the alignment result corresponding to the query DNA sequence is found in the reference DNA sequence set.

[0065] Example 1

[0066] Building upon Specific Implementation Method 1, this paper further elaborates on how to accurately detect single nucleotide variants (SNVs) and short insertions or deletions (Indels) in the query DNA sequence Q relative to the reference DNA sequence R using the core concept of this invention. This method introduces a bidirectional query and anchor point localization strategy to achieve higher resolution accuracy. The specific process is as follows:

[0067] Step 1: Building a bidirectional index:

[0068] To achieve bidirectional comparison, a composite index is constructed according to S1 in Specific Implementation Method 1 for the query DNA sequence Q and the reference DNA sequence R, denoted as... and This step forms the basis for subsequent two-way verification and precise localization, allowing the algorithm to not only query the match of Q in R, but also the match of R in Q.

[0069] Step 2: Calculate the forward MS array:

[0070] According to the description in S2 of the first specific implementation method, using Generate an array MS of the query DNA sequence Q relative to the reference DNA sequence R, denoted as... This array reflects the matching status of each position on the queried DNA sequence Q in the reference DNA sequence R, and is the main basis for identifying mutation events.

[0071] Step 3: Identify candidate regions of variation (significant descent points);

[0072] Traversal Actively seek out "significant decline points". If the location... satisfy (Indicating that the previous position is a significant match) and (Indicates that the current position matching was suddenly interrupted), position This was defined as a significant drop point. This abrupt signal mutation strongly suggests the location of the queried DNA sequence. Initially, a mutation event occurs relative to the reference DNA sequence, then the location... It is considered a candidate position for the left boundary of the variation.

[0073] Step 4: Locate the unique anchor point to determine the right boundary:

[0074] For each identified significant descent point The algorithm at the significant descent point In the query DNA sequence on the right, find the nearest anchor point that uniquely matches a certain k-mer in the reference DNA sequence R. A unique match means that the k-mer appears only once in the entire reference DNA sequence R, which can be determined by querying... This unique anchor point acts like a precise GPS coordinate, reliably re-locking the right boundary of the variant region to a specific location in the reference DNA sequence. This effectively avoids the positioning ambiguity caused by repetitive genomic sequences and is crucial for ensuring the accuracy of variant detection.

[0075] Step 5: Bidirectional Precise Analysis of the Variation:

[0076] With the left boundary (significant descent point) ) and right boundary (unique anchor point) The mutation event is confined to a specific interval. Next, variant sequences are extracted from the query DNA sequence. These variant sequences refer to the subsequences of the query DNA sequence Q between its left and right boundaries, also known as the mutated query alleles. Variant sequences are extracted from the reference DNA sequence. To determine the original sequence (also known as the reference allele) corresponding to the variant in the reference DNA sequence R, it is necessary to use... Reverse calculation from anchor point Initially, the match between the reference DNA sequence R and the query DNA sequence Q is determined. This reverse calculation allows for the precise location of the mutation's initiation position within the reference DNA sequence R. and end position The reference allele is... Finally, through comparison and By knowing the sequence content and length, the type of mutation and its precise coordinates on the reference DNA sequence R can be accurately defined.

[0077] like and , indicating a single nucleotide variant (SNV). If Empty or length much smaller , indicates a short insertion. If Empty or length much smaller , indicates short deletion (Indel).

[0078] Example 2

[0079] To verify the effectiveness of the method of this invention, three benchmark tests covering typical bacterial genomics analysis tasks were designed and implemented, and real data were used for evaluation. The experiments were conducted on a portable computing device equipped with an Intel i7-1165G7 quad-core processor and 16GB of memory. All tools were run at default parameters to ensure the objectivity and reproducibility of the comparisons.

[0080] Variance Detection Dataset: This invention utilizes seven closely related Staphylococcus aureus genome assemblies from a study published in Access Microbiology by Wick et al. (2025). These assemblies originate from the same strain but were generated using different sequencing technologies (including short and long reads) and assembly software, thus representing the data diversity that may be encountered in practical applications. The ground-truth set is based on a high-quality, mixed-assembly 2.9 Mb Staphylococcus aureus reference genome, obtained using high-precision short-read alignment methods (BWA-MEM and freebayes), which is considered the "gold standard" under current technological conditions.

[0081] Gene search dataset: This invention utilizes a large-scale strain library from the Norwegian Resistant Microbes Monitoring Project (NORM), from which 1999 whole chromosome sequence assemblies of *E. coli* (average length 4.63 Mb) were selected. The query sequences were 19 gene sequences (approximately 50,000 bases in total length) extracted from the pks pathogenicity island of the *E. coli* standard strain IHE3034. The actual gene locations were generated using bakta, a recognized rapid annotation tool in the field, as a benchmark for evaluating the accuracy of the gene search in this invention.

[0082] Reference sequence alignment dataset: This invention uses 112 whole-chromosome assembled sequences of ST131 Escherichia coli from the NORM project, a sequence type with significant clinical importance globally. Actual alignment results were obtained through the following procedure:

[0083] First, a high-quality hybrid assembly was randomly selected as a reference sequence. Then, the raw paired-end short-read sequencing data of the remaining 111 strains were aligned to this reference sequence using BWA-MEM2. Finally, the samtoolspileup function was used to generate a consensus sequence for each strain relative to the reference. This workflow simulates the standard analysis process based on raw sequencing data, and its results can serve as the gold standard for evaluating assembly-based alignment methods.

[0084] Benchmark comparison method:

[0085] To comprehensively evaluate the performance of this invention, several advanced prior art technologies representative of various task domains were selected as benchmarks for comparison:

[0086] Variation detection task: Performance comparison with MUMmer4 (based on whole-genome alignment), Shred (simulating short-read alignment process), and SKA (based on split k-mer alignment) – an alignment-free tool.

[0087] Gene search task: Speed ​​and accuracy compared to blast+, the most widely used tool in the field of biological sequence search.

[0088] Reference sequence alignment task: Performance comparison with SKA and Snippy, which are also fast alignment tools oriented towards assemblies.

[0089] Performance evaluation metrics:

[0090] This invention employs a series of quantitative indicators to comprehensively evaluate the performance of the invention and benchmark methods from multiple dimensions, including accuracy, sensitivity, and resource consumption.

[0091] Accuracy and sensitivity metrics:

[0092] Precision: Precision measures the proportion of true variants among the variants reported by the method. This metric directly reflects the reliability of the variant detection results.

[0093]

[0094] in, This represents the number of true variants that were correctly detected. This indicates the number of non-variant sites reported in error reports.

[0095] Recall / Sensitivity: Recall represents the proportion of all real variants that are correctly detected by the model. This metric assesses the model's ability to find real variants.

[0096]

[0097] in, This indicates the true number of variants that the model failed to detect.

[0098] Alignment error rate: In the reference sequence alignment task, the total number of inconsistent bases (including A, C, G, T, and empty or ambiguous bases) is calculated by comparing the shared sequence generated by this invention with the gold standard shared sequence position by position. The error rate is reported in two forms: the absolute number of erroneous bases and the relative error rate relative to the total length of the reference genome.

[0099] Calculate resource consumption indicators:

[0100] Running time: Records the total wall-clock time required from the start of the program to the completion of the task.

[0101] Peak memory usage: Records the maximum amount of resident memory (RSS) used during program execution.

[0102] Temporary disk space: Records the maximum disk space occupied by temporary files or directories used by the program during its operation.

[0103] Experimental Results and Analysis

[0104] The performance of the method proposed in this invention was evaluated in three typical application scenarios and compared with existing technologies. The experimental results are shown in the table below.

[0105] Performance of mutation detection task:

[0106] In the task of detecting variations in the Staphylococcus aureus genome, the method of the present invention has demonstrated excellent performance in the detection of both single nucleotide variants (SNVs) and insertions / deletions (Indels).

[0107] Table 1. Comparison of detection performance for single nucleotide variants (SNVs)

[0108]

[0109] Table 2. Performance comparison of insertion / missing information detection (SKA does not support insertion detection).

[0110]

[0111] Experimental results show that the method of this invention achieves 100% sensitivity, similar to MUMmer4 and Shred, successfully detecting all real SNVs and Indels. However, in terms of precision, this invention significantly outperforms the benchmark methods. Its precision for SNV detection (0.32) is more than twice that of MUMmer4 (0.14) and Shred (0.16), and its precision for Indel detection (0.26) is nearly three times that of the latter two. This demonstrates the significant advantage of this invention in effectively controlling false positives. The SKA method performs poorly on this dataset, exhibiting low sensitivity and precision.

[0112] Performance of gene search task:

[0113] In a task involving the search for PKS pathogenic island genes in 1999 *E. coli* genomes, the method of this invention demonstrated excellent performance in both speed and accuracy. This method used four parallel processes to process all samples, with a total processing time of 509 seconds (average 0.25 seconds per sample) and a peak memory usage of 195 MB. In comparison, BLAST+ took 364 seconds (average 0.18 seconds per sample) and had a peak memory usage of 189 MB. Regarding accuracy, in 934 samples containing PKS gene islands, this method found gene regions that were completely consistent with Bakta annotations in 894 (95.7%). In the remaining 39 samples, the average difference between the aligned regions reported by this invention and the actual annotations was only 5 bases.

[0114] Reference sequence alignment task performance:

[0115] In the task of aligning reference sequences of 112 E. coli genomes, the method of this invention outperforms existing technologies in terms of accuracy.

[0116] Table 3. Number of absolute error bases in the alignment of 112 E. coli genomes.

[0117]

[0118] Table 4. Relative error rate of 112 E. coli genome alignments

[0119]

[0120] Experimental results show that the method of this invention generates 73,814 median erroneous bases in the shared sequence (relative error rate of 1.5%), significantly lower than SKA (104,990, 2.1%) and Snippy (149,894, 3.0%), demonstrating that this invention has higher accuracy in generating high-quality alignments based on assemblies. In terms of computational resource consumption, the total time of this invention is 624 seconds (average 5.7 seconds per sample), comparable to SKA (669 seconds, 6.0 seconds), and much faster than Snippy (2361 seconds, 21.1 seconds).

[0121] In summary, the method of the present invention demonstrates superior performance compared to existing advanced technologies in terms of accuracy, sensitivity, or computational efficiency in multiple key application scenarios, proving the effectiveness of the present invention.

[0122] The above examples of the present invention are merely illustrative of the computational model and process of the present invention, and are not intended to limit the implementation of the present invention. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is impossible to exhaustively list all possible implementations here. Any obvious variations or modifications derived from the technical solutions of the present invention are still within the scope of protection of the present invention.

Claims

1. A method for local alignment of DNA sequences based on k-bounded matching statistics, characterized in that: It includes the following steps: S1. Construct a composite index using the Burrows-Wheeler transformation of the spectral data to a set of reference DNA sequences; S2. Use the composite index to scan the bases of the query DNA sequence one by one to generate a k-bounded matching statistics array. The k-bounded matching statistics array includes the matching length between the query DNA sequence and the reference DNA sequence set. S3. Establish a random matching model, and based on the k-bounded matching statistics array, use the random matching model to generate a derandomized matching statistics array. S4. Define syntax rules and translation rules, process the derandomized matching statistics array using syntax rules and translation rules to obtain character-type alignment events of the query DNA sequence, and find the alignment result corresponding to the query DNA sequence in the reference DNA sequence set based on the character-type alignment events.

2. The method for local DNA sequence alignment based on k-bounded matching statistics according to claim 1, characterized in that: The specific process of S1 is as follows: S11. Obtain a set of reference DNA sequences and set the length of the k-mer to k and the total number to n; S12. Construct an initial k-spectrum set based on the reference DNA sequence set. Select a subset from the initial k-spectrum set, and then select all k-mers from the subset. Length prefix If a subset is not a suffix of any other k-mer in the initial k-spectrum set, add Z virtual k-mers prefixed with the special placeholder $ to the subset. Then merge the subset with the virtual k-mers with the initial k-spectrum set after filtering the subset to obtain a k-spectrum. S13. Based on the k-spectrum, the SBWT index is obtained by using the Burrows-Wheeler transform of the spectrogram. The SBWT index stores the results of all k-mers sorted by co-word order. In addition, calculate and store an array of the longest common suffixes corresponding to the k-mer co-word order in the SBWT index. The array of the longest common suffixes includes the length of the longest common suffix of each k-mer and the previous k-mer under the k-mer co-word order. A composite index is formed by combining the SBWT index and the longest common suffix array.

3. The method for local DNA sequence alignment based on k-bounded matching statistics according to claim 2, characterized in that: The specific process of S2 is as follows: S21. Obtain a query DNA sequence and create an integer array MS with the same length as the query DNA sequence; S22. Initialize the matching length of the query DNA sequence. And an initial matching range representing all k-mers in the composite index. ,in, To locate a position on a DNA sequence, The initial value is the entire interval of all k-mers; S23. Using the composite index from S13, scan the bases of the DNA sequence one by one to generate a k-bounded matching statistics array.

4. The method for local DNA sequence alignment based on k-bounded matching statistics according to claim 3, characterized in that: The specific process of S23 is as follows: Perform an extended-contraction matching search in each scan: (a) Attempt to match extensions: （1） in, To query the first DNA sequence The bases at each position, To query the first DNA sequence The matching range corresponding to each position To query the first DNA sequence The matching range corresponding to each position Indicates extension to the right; (ii) If Returns a non-empty range, indicating that the matching continues, and retrieves the current position. Corresponding matching range Current location The matching length is updated to and ; (iii) If Returns an empty range, indicating a match within the specified length. Interrupt, execute in a loop and decrease ,in, This indicates that the left side is shortened. After each shortening, step (1) is repeated until the attempt to extend is successful. This gives the current matching range and matching length; Based on the above scanning and matching process, the obtained positions Match length Store in array Similarly, after the DNA sequence has been scanned and matched for all bases, the matching lengths of all positions in the DNA sequence are stored in the array MS. The resulting array MS is the k-bounded matching statistics array.

5. The DNA sequence local alignment method based on k-bounded matching statistics according to claim 4, characterized in that: The random matching model in S3 is based on the following: When the query DNA sequence is not associated with the reference DNA sequence, the length of a single random match follows a geometric distribution, and the longest random match length in the entire set of reference DNA sequences follows an extreme value distribution.

6. The method for local DNA sequence alignment based on k-bounded matching statistics according to claim 5, characterized in that: The specific process of S3 is as follows: S31. Customize the statistical false positive rate tolerance, and determine the significance threshold of the match length based on the total number of k-mers and the statistical false positive rate tolerance; S32. Create a derandomized matching statistics array DMS of the same length as the k-bounded matching statistics array, and initialize the last value of the derandomized matching statistics array DMS. , Given the sequence length, a significance threshold based on the matching length is then applied. Traverse the k-bounded matching statistics array in reverse order from the second-to-last position to the first position of the query DNA sequence; (1) If ,and ,think reliable, ; (2) If ,but For noise; like ,but To match the plateau period; like ,but To match the breakpoint; when When considering noise, matching plateau period, or matching interruption point, it is assumed that... Unreliable ; The process continues until the traversal is complete, at which point the derandomized matching statistics array DMS is obtained.

7. The method for local DNA sequence alignment based on k-bounded matching statistics according to claim 6, characterized in that: The specific process of S31 is as follows: （2） in, Represents the maximum matching length that is generated entirely randomly. Significance threshold exceeding match length The probability, For custom statistical false positive rate tolerance; The cumulative distribution function based on the extreme value distribution is: （3） in, The probability of two nucleotides randomly matching; Solving according to equations (2)-(3), we get: （4）。 8. The method for local DNA sequence alignment based on k-bounded matching statistics according to claim 7, characterized in that: The specific process of S4 is as follows: S41. Define the syntax rules, which are to decode the continuous array DMS into discrete character-type comparison events; S42. Define translation rules, which are mapping tables that associate local numerical patterns of different shapes in the DMS array with specific biological events; S43. Process the derandomized matching statistics array using syntax rules and translation rules to obtain character-type alignment events in the query DNA sequence. Based on the character-type alignment events, find the alignment result corresponding to the query DNA sequence in the reference DNA sequence set.

9. A method for local DNA sequence alignment based on k-bounded matching statistics according to claim 8, characterized in that: The specific process of S41 is as follows: Iterate through the array DMS, at each position According to , and The local numerical patterns are composed of preset corresponding character-type alignment events, including mismatches or single base insertions, single base deletions or structural breakpoints, matching and querying DNA sequence-specific regions.

10. A method for local DNA sequence alignment based on k-bounded matching statistics according to claim 9, characterized in that: The specific process of S42 is as follows: 1) For V-shaped valley patterns caused by mismatches or single base insertions: When the array DMS is querying DNA sequences If the position exhibits a pattern of positive numbers, zero or negative numbers, or 1, and this pattern is interpreted as a single-point mutation causing a matching interruption and rapid recovery, then the position... Translated as Event X; 2) Cliff mode for single base deletions or structural breakpoints: When the array DMS is querying DNA sequences The position shows values ​​ranging from above the significance threshold to below the significance threshold. The pattern, which is interpreted as a discontinuous break in a significantly homologous matching region, corresponding to a structural difference between the query sequence and the reference sequence, will then be used to determine the position. and location Translated as the left and right boundaries of the R event; 3) Plain mode for querying DNA sequence-specific regions: when If it is not part of the V-shaped valley pattern, it is interpreted as not having a corresponding match in the reference DNA sequence, and the position is... Translated as Event I; 4) For matching platform or ramp modes: when If none of the above three patterns are satisfied, and the result is interpreted as a continuous homologous match, then the position will be... Translated as the M event.

Images