Parallel precise local sequence alignment method based on layered filtering and adaptive parameter optimization

Abstract

The invention discloses a parallel precise local sequence alignment method based on layered filtering and adaptive parameter optimization, and belongs to the technical field of biology. The method aims to solve the problem that speed and accuracy cannot be both considered in the existing local comparison field technology. The method comprises the following steps: establishing a probability model containing a false positive rate and a false negative rate, obtaining the total cost, and determining an optimal k-mer length and a counting threshold which enable the total cost to be minimum based on comparison requirements of a user; according to the method, producer-consumer parallel computing is adopted, a producer thread performs preliminary screening on a query sequence by using a global and coarse-grained IBF index, and tasks with a potential matching relationship are packaged into a cart; a consumer thread obtains a'cart 'from the queue, a finer index is constructed in a local range for accurate filtering, and finally an accurate banded Waterman-Eggert comparison algorithm is executed only on a minimum candidate area; and finally, performing duplicate removal on the result to obtain full-amount accurate local comparison.

Description

Technical Field

[0001] This invention belongs to the field of biological technology, specifically relating to a method for local sequence alignment between genome, transcriptome, and metagenomic sequence databases. Background Technology

[0002] Accurate and rapid searching for locally similar sequences (i.e., local alignment) in massive genomic datasets is a fundamental and crucial computational task in modern computational biology and genomics research. The pressures of evolutionary selection have resulted in high conservation of functional DNA fragments (such as gene coding regions and regulatory elements) across species at the sequence level. Therefore, large-scale local alignment can systematically identify these conserved regions across the entire genome, providing crucial clues for gene function prediction, evolutionary relationship reconstruction, and disease association studies. However, with the maturity of third-generation sequencing technologies (such as PacBio and Oxford Nanopore), which can generate millions of reads exceeding 10kb in a single sequencing run, and milestone events such as the T2T consortium's completion of the complete human genome map, the scale and complexity of data have increased exponentially, posing a severe performance challenge to existing local alignment techniques.

[0003] Current mainstream local alignment methods have inherent technical trade-offs in the three core indicators of "accuracy", "speed" and "resource consumption", which are specifically reflected in:

[0004] Limitations of Accuracy and Sensitivity of Heuristic Methods: Heuristic methods, represented by the BLAST family of tools, are currently the most widely used strategy. Their core idea is to locate candidate regions through fast "seed" matching (usually short, perfectly matching k-mers), and then execute time-consuming dynamic programming algorithms only within these candidate regions. To ensure speed, BLAST uses a longer k-mer as the seed by default (e.g., 28-mer). However, this results in insufficient sensitivity to homologous regions containing multiple mismatches, insertions, or deletions (i.e., large edit distances). Experiments show that BLAST with default parameters misses a significant proportion of biologically significant local alignments. If users shorten the seed length to improve sensitivity, the number of seed matches increases exponentially, making subsequent extension and verification steps computationally overloaded. This not only significantly increases runtime but also generates a massive number of statistically meaningless short matches, interfering with the recognition of real biological signals.

[0005] Computational Performance and Memory Bottlenecks of Precise Alignment Methods: In contrast to heuristic methods, precise alignment methods guarantee the integrity of results. The theoretical basis of these methods is the use of a "lossless filtering" strategy to ensure that all target alignments meeting specific scoring criteria are found. For example, the early Stellar tool used a filter called SWIFT, combined with the k-mer counting lemma, to achieve precise results. However, this comes at the cost of extremely high computational costs. Because their algorithms are single-threaded, runtimes for processing modern mammalian genomes are typically measured in hours or even days. More seriously, the memory model of these algorithms makes effective parallel scaling difficult. Even with simple parallelization schemes (e.g., processing the genome in chunks independently), each parallel task requires loading a huge index or constructing complex intermediate data structures. The total memory requirement increases linearly with the number of cores, easily leading to Out of Memory errors, which severely limits the practical application of precise alignment methods on large-scale datasets.

[0006] The complexity of parameter selection and the challenges of handling repetitive sequences: The performance of almost all alignment tools heavily depends on the selection of a series of key parameters, such as the length and shape of the k-mer (whether it includes spacers), the scoring matrix, and the penalty. The optimal values ​​of these parameters vary with the evolutionary distance of the species, the error rate of the sequence, and the user's research objectives. Their selection often requires users to possess deep expertise and conduct extensive trial-and-error experiments, posing a significant technical barrier for non-expert users. Furthermore, all eukaryotic genomes contain a large number of repetitive sequences (such as SINEs and LINEs). These regions generate massive amounts of non-specific seed matching, leading to a large amount of redundant computation in these regions and significantly reducing overall efficiency. The traditional approach is "hard masking," which involves pre-deleting or replacing these regions with 'N' using tools such as RepeatMasker. However, increasing research shows that many repetitive sequences themselves carry important biological functions, and brute-force masking can lead to the permanent loss of this information, resulting in systematic bias in the analysis results.

[0007] In summary, existing technologies fail to provide an effective solution that simultaneously meets the five core requirements of "completely accurate results," "heuristic-level speed," "controllable parallel memory," "automatic parameter optimization," and "intelligent handling of repetitive sequences." Therefore, a new technological paradigm is urgently needed in this field to address the stringent demands of large-scale sequence alignment in the post-genomic era. Summary of the Invention

[0008] This invention aims to address the problem that existing local comparison techniques cannot simultaneously achieve both speed and accuracy.

[0009] A parallel, precise local sequence alignment method based on hierarchical filtering and adaptive parameter optimization includes:

[0010] S100. Determine ε-match based on the comparison criteria set by the user and the length of the input sequence, and determine the combination cost based on the false positive rate (FPR) and false negative rate (FNR). By minimizing the combinatorial cost, a set of theoretically optimal k-mer filtering parameters for the current task is determined, including the k-mer length k, shape, and shared k-mer counting threshold. ;

[0011] S300, two-level hierarchical filtering, and parallel distribution of "producers-consumers":

[0012] Global coarse-grained filtering: executed in parallel by multiple producer threads; all reference fragments are logically divided into different "warehouses," each warehouse being a reference warehouse, and a globally shared IBF index is built for all warehouses; all parallel producer threads access the same IBF index to determine whether the kmer in the query sequence exists in a certain warehouse, based on a threshold. Filter out some combinations of "query fragment - reference repository";

[0013] Task Packaging and Distribution: The producer thread will package the "query fragment - reference warehouse" tasks through the first layer of filtering into a data structure called "cart". The cart will be placed into a global, first-in-first-out task queue.

[0014] Local fine-grained filtering and precise comparison: This is executed in parallel by multiple consumer threads. Each consumer thread retrieves a "cart" from the task queue and, for all query segments contained within the "cart," immediately builds a local full-text index in memory that records the precise position of the k-mer in the query segment. Using the full-text index, a SWIFT filter is executed. After filtering, the results are fed into the Waterman-Eggert dynamic programming algorithm to find all local comparisons between two sequences that exceed a threshold score.

[0015] S400: Collect the local alignment results output by all consumer threads, identify and merge duplicate alignments based on the aligned genomic coordinates, and finally output the results.

[0016] Furthermore, the comparison criterion is that the comparison length must be greater than or equal to the minimum comparison length. And the number of edit operations compared must be less than or equal to the maximum number of edit operations. .

[0017] Furthermore, the number of edit operands must be less than or equal to the maximum number of edit operands. The judgment method is as follows:

[0018] For a length of ,Include The comparison of each editing operation is performed by maximizing the normalization of editing operation conditions. To make a judgment.

[0019] Furthermore, the false negative rate (FNR) is obtained using a three-dimensional dynamic programming algorithm, including:

[0020] A DP table, which is a three-dimensional dynamic programming matrix, is pre-computed and populated to store the number of false negative error configurations. The value of FNR is determined by querying this pre-computed DP table during the process of optimal k-mer filtering parameters.

[0021] DP table status middle Represents the length of the sequence. Represents the total number of editing operations. Represents the length of k-mer, Represents the shared k-mer count threshold required for filtering; state as follows:

[0022]

[0023] The total number of error distributions is determined by the number of combinations. And thus obtain .

[0024] Furthermore, the false positive rate (FPR) is obtained through the following steps:

[0025] First, a sequence of length L is modeled as a set of k-mers contained in the sequence;

[0026] The expected number of occurrences of a specific “canonical k-mer” as follows:

[0027]

[0028] in: The effective sequence size correction coefficient is mentioned above; Represents the size of the alphabet. This represents the total number of k-mers in the sequence;

[0029] The probability of a single k-mer producing a "spurious match". :

[0030]

[0031] False positive rate due to hash collisions in IBF data structures and The total false positive probability of a single k-mer is obtained. Thus, the probability is obtained. :

[0032]

[0033] in, This represents the total number of k-mers in the pattern; This represents what is observed in the pattern. The cumulative probability of one or fewer k-mer hits; Representative observed The probability of hitting k-mer multiple times is the final estimate of the FPR.

[0034] Furthermore, during the task packaging and distribution process, dynamic priority scheduling of the repeating sequence is first performed, and then the producer packages the priority-scheduled candidate tasks into the cart; the process of dynamic priority scheduling of the repeating sequence includes:

[0035] When determining candidate tasks, if the number of hit warehouses for a query fragment Q exceeds half of the total number of warehouses, the system determines that the query fragment Q is highly repetitive. For such queries, the system prioritizes pairing the query fragment Q with the reference warehouse with the highest overall variability among all hit warehouses to form candidate tasks, thereby directing computing resources to the region most likely to generate specific alignments.

[0036] Furthermore, the overall variability calculation process of the warehouse includes:

[0037] While constructing the global IBF index, the variability score of the reference sequence fragments in each bin is calculated. It is defined as the number of unique k-mers within the reference sequence segment i. Total number of k-mers within the fragment The ratio; the overall variability score obtained by averaging the variability scores of all reference segments contained in a given warehouse.

[0038] Furthermore, during the global coarse-grained filtering process, multiple producer threads are started based on the number of available CPU cores; each producer thread is responsible for processing a portion of the query fragments; for each fragment, the producer extracts the k-mer and queries the global IBF index, accumulating the shared k-mer hit count between each "query fragment-reference warehouse" pair; if this number exceeds the threshold t determined by S100, then the "query fragment-reference warehouse" pair is considered a candidate task.

[0039] Furthermore, the construction method of full-text index is as follows:

[0040] Multiple consumer threads independently retrieve carts from the queue. For each task in a cart, the consumer first constructs a local and temporary full-text index I in its thread-private memory, which contains only the k-mer precise position of the query segment required by the current task. Local means that this index I is limited to the batch of query segments contained in the current cart. Index I is constructed in the private memory of the consumer thread in real time. When the task in the cart is finished, index I is discarded and the memory is released.

[0041] Furthermore, in the process of two-level hierarchical filtering and parallel distribution of "producers and consumers", a symmetrical k-mer pattern with added intervals is adopted.

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

[0043] 1. Achieving both accuracy and high computational speed: This invention employs a mathematically rigorous, lossless filtering process to ensure the identification of all accurate local alignments that meet user-defined conditions, achieving theoretical 100% sensitivity. Simultaneously, through innovative two-stage filtering and a highly parallelized architecture, its operating speed is up to 900 times faster than its predecessor (the single-threaded accurate alignment tool Stellar), reaching a speed comparable to mainstream heuristic tools like BLAST. This systematically solves the long-standing trade-off between speed and accuracy in the field of local alignment.

[0044] 2. Innovative Parallel Architecture Achieves Exceptional Space Efficiency: The "producer-consumer" model proposed in this invention, combined with a strategy of "global probabilistic indexing + local real-time precise indexing," is key to its space efficiency. The global IBF index itself is extremely compact, while the computationally and memory-intensive precise index is dynamically and temporarily created only when the consumer thread processes a small batch of data. This architecture avoids the problem of rapidly increasing memory consumption caused by each thread needing to load a huge index in traditional parallelization methods, making it possible to perform genome-wide precise alignments on standard multi-core servers.

[0045] 3. Automated and intelligent selection of filtering parameters: This invention introduces for the first time a mechanism for automatically deriving filtering parameters based on the FPR / FNR probability model in the field of precise local alignment. This mechanism transforms the complex, experience-dependent parameter selection problem into a calculable mathematical problem with a clear optimization objective, greatly improving the ease of use of the tool and its performance stability and robustness under different species and alignment requirements.

[0046] 4. This invention provides an efficient and lossless strategy for processing repetitive sequences: Compared to traditional hard masking methods that permanently lose sequence information, the dynamic priority scheduling strategy based on "variability" proposed in this invention is a non-destructive processing method. Without discarding any sequence information, it intelligently guides computing resources to regions with richer information, effectively avoiding time traps in repetitive sequences and achieving complete and efficient analysis of the entire genome (including repetitive regions). Attached Figure Description

[0047] Figure 1 The flowchart is for an example.

[0048] Figure 2 This is a processing logic diagram for an example.

[0049] Figure 3 This is a schematic diagram of an IBF cart queue with 4 bins. Detailed Implementation

[0050] This invention relates to the fields of large-scale bioinformatics, computational genomics, and high-performance parallel computing, aiming to address the enormous computational overhead and memory consumption caused by sequence length and repetitive sequences when performing precise local alignment on whole-genome-level data. The purpose of this invention is to overcome the contradiction between the slow speed of traditional precise alignment tools and the fact that heuristic tools (such as BLAST), while ensuring speed, may miss a large number of significant alignments or be overwhelmed by massive amounts of low-quality matches. Therefore, the core objective of this invention is to provide a novel precise local sequence alignment system and method that integrates parallel computing, efficient indexing, and intelligent parameter optimization, systematically addressing the multiple technical challenges faced by existing technologies in processing large-scale genomic data, including speed, memory usage, accuracy, and ease of use.

[0051] This invention relates to a parallel, spatially efficient, and completely accurate local sequence alignment method for ultra-large-scale genome, transcriptome, and metagenomic sequence databases. More specifically, this invention relates to an integrated system and method that incorporates a two-level hierarchical filtering strategy, a producer-consumer asynchronous parallel computing model, and an adaptive parameter derivation mechanism based on a rigorous probabilistic model. This aims to ensure lossless alignment results while achieving computational speeds comparable to, and even surpassing, mainstream heuristic tools in certain scenarios. This invention also covers specific applications of this system and method in advanced bioinformatics analysis, including but not limited to: deep comparative genomics between distantly related species, precise breakpoint resolution of complex structural variations (SVs), personalized genome assembly assistance based on long-read sequencing data, and whole-genome scanning of conserved non-coding regions (CNEs) and other functional elements. Specific implementation method one:

[0053] The parallel precise local sequence alignment method proposed in this embodiment, based on hierarchical filtering and adaptive parameter optimization, has a key innovation in its computational framework: a two-level hierarchical filtering based on a producer-consumer architecture, guided by an adaptive parameter optimization model. This framework decomposes a computationally intensive, large-scale alignment problem into a series of well-defined, computationally decoupled, and resource-controllable subtasks. The specific process includes:

[0054] Step 1: Adaptive Parameter Derivation: Before executing the comparison task, the system first runs an automated parameter derivation program. This program is based on the user-defined comparison criteria (i.e., the comparison length must be greater than or equal to the minimum comparison length). And the number of edit operations compared must be less than or equal to the maximum number of edit operations. A probabilistic model is constructed using the alignment criteria (referred to as ε-match) and the length of the input sequence, including the false positive rate (FPR) and false negative rate (FNR). FNR refers to the probability that a true ε-match is incorrectly discarded by the filter, while FPR refers to the probability that two sequence segments with no ε-match relationship share more than a threshold number of occurrences by chance. k-mer (i.e., length of) The probability that the sequence is incorrectly identified as a candidate region by the filter.

[0055] By minimizing the combinatorial cost function FPR+FNR, the program automatically derives and determines a set of theoretically optimal k-mer filtering parameters for the current task, including the length and shape of the k-mers and the shared k-mer count threshold required for subsequent filtering. This step automates the complex parameter tuning process, eliminating the need for manual user intervention.

[0056] Step Two: Two-Level Hierarchical Filtering and Parallel "Producer-Consumer" Distribution: This invention employs a highly scalable parallel computing architecture to decouple the complex comparison task into two main stages:

[0057] Global coarse-grained filtering: This is executed in parallel by multiple "producer" threads. First, the reference sequence is segmented, and all resulting segmented reference sequence fragments are placed into different bins, with the number of bins less than the number of reference sequence fragments. A globally shared, space-efficient Interleaved Bloom Filter (IBF) index is built for all bins. Then, the producer threads read the query sequence fragments, perform a fast k-mer membership lookup using the IBF index, and calculate the shared k-mer count threshold derived in step one. This quickly filters out the vast majority of "query fragment-reference bin" combinations that do not contain enough shared k-mers. This is a probabilistic but extremely fast global filtering process.

[0058] Task Packaging and Distribution: The producer thread packages the "query fragment-reference warehouse" tasks with potential homologous relationships, filtered through the first layer, into a data structure called a "cart". A "cart" does not contain a single "query fragment-reference warehouse" pair, but rather a set of query fragments and a reference warehouse. These query fragments were found to match the same reference warehouse in the IBF index filtering. The carts will then be placed into a global, first-in-first-out task queue.

[0059] Local fine-grained filtering and precise comparison: This is executed in parallel by multiple "consumer" threads. Each consumer thread retrieves a "cart" from the task queue and, for all query segments contained within that "cart," immediately builds a local index in memory that records the precise position of each k-mer within the query segment. This index is called a full-text index because, unlike the first-level global IBF index, it records the precise position of each k-mer within the query segment. This index is also local, as it is not built for all sequences but only for the current batch of sequences within the cart, and it is temporary; the index is cleaned up after subsequent tasks have used it. Using this full-text index, a second-level lossless SWIFT filter (an efficient lossless filtering algorithm based on the k-mer counting lemma) is executed for more rigorous collinearity and density checks. Only a very small number of candidate regions that pass through this second layer of fine filtering will be finally sent to a strip-shaped Waterman-Eggert dynamic programming algorithm module for the final, base-by-base precise local alignment calculation. The Waterman-Eggert dynamic programming algorithm is used to find all local alignments between two sequences that score above a threshold. Strip-shaped dynamic programming is a standard dynamic programming optimization technique that greatly accelerates the alignment calculation by calculating only a narrow "strip" region around the diagonal of the alignment matrix.

[0060] Step 3: Result Integration and Output: After all parallel consumer threads have completed their respective alignment tasks, the system enters the final integration phase. In this phase, the local alignment results generated by all consumer tasks are collected. Because the input sequence was segmented into partially overlapping fragments in the initial stage, the same alignment may be found repeatedly by different threads within the overlapping region. Therefore, a key step in this phase is to perform a merging process to eliminate duplication, thereby ensuring that the final alignment result output to the user is accurate and unique.

[0061] The efficient processing strategy for repetitive sequences described above quantifies the degree of repetition by calculating the "variability" (unique k-mer ratio) of each reference sequence segment. When a query segment is identified as having high repetition (i.e., it has potential matches with a large number of reference segments), the system intelligently prioritizes processing reference segments with higher variability (i.e., lower repetition), thereby effectively bypassing computationally intensive regions and avoiding unnecessary computational overhead.

[0062] Example

[0063] Combination Figure 1 and Figure 2 This embodiment describes the search for all genomes that meet specific criteria between two very large mammalian genomes: human GRCh38 and mouse GRCm39. (For example, minimum alignment length) =100 bp, maximum number of edit operations Taking a precise local alignment as an example, the steps include:

[0064] Step 1: Determine the intelligent adaptive filtering parameters:

[0065] Step 11: Construct a formal problem definition and cost function.

[0066] The goal of this invention is to find all local alignments that satisfy the user-defined "ε-match" criterion. This criterion is strictly defined as: a segment of length... ,Include A pair of edit operations is formed if and only if they simultaneously satisfy the minimum length condition. and maximum normalized editing operation conditions When this occurs, it is considered a valid "ε-match".

[0067] To find the optimal filtering parameters, the system first establishes a combined cost function to be optimized:

[0068]

[0069] This function aims to find an optimal balance point. It serves as a weighting factor to balance the two filtering errors: the false positive rate (FPR) and the false negative rate (FNR). The length of k-mer, The threshold for counting shared k-mers required for filtering, i.e., the two sequence fragments must share at least one k-mer. Only when a certain number of k-mers are detected will it be considered a possible match and retained by the filter for further analysis.

[0070] Steps 1 and 2: Combinatorial mathematical modeling and dynamic programming solution of FNR (false negative rate).

[0071] To accurately estimate FNR, a true ε-match is defined as one where the number of shared k-mers is below the filtering threshold. The probability of false rejection is calculated using a three-dimensional dynamic programming algorithm. This algorithm aims to provide a more accurate estimate, overcoming the limitation of the k-mer counting lemma, which only applies to the worst-case scenario where the error probability follows a uniform distribution. The core of this algorithm is the pre-computation and filling of a dynamic programming (DP) table, a three-dimensional matrix that stores the number of false negative configurations. This DP table only needs to be pre-computed once (taking approximately 1 second and consuming about 6MB of memory) and can be repeatedly consulted in all subsequent parameter derivation tasks. The state of the DP table... Strictly defined as: in a length of Randomly introduced in the sequence One edit operation caused the sequence to share fewer k-mers with the original sequence than a threshold. The total number of error distribution patterns.

[0072] state The symbols in the Chinese language are defined as follows: Represents the length of the sequence (i.e., the minimum length of the ε-match). ); Represents the total number of editing operations (i.e., errors); Represents the length of k-mer; This represents the threshold number of shared k-mers required for filtering. The three-dimensional dynamic programming algorithm efficiently populates this DP table using a rigorous recursive formula.

[0073]

[0074] Finally, for a given combination of parameters The FNR value can be quickly obtained by querying this pre-computed DP table and normalizing it (i.e., dividing by the total number of error distributions). The total number of error distributions is determined by the number of combinations. Given:

[0075]

[0076] Therefore, the final formula for calculating FNR is:

[0077]

[0078] Step 13: Probabilistic statistical modeling of FPR (false positive rate) and correction of effective sequence length.

[0079] To accurately estimate the FPR (First-Place Ratio), which is the ratio of the number of frames that two query sequences and reference sequence segments with no real alignment share more than a certain percentage of frames by chance. To address the probability of a reference fragment being incorrectly identified as a candidate region due to the presence of k-mers, this invention establishes a probabilistic statistical model. This model first models a sequence of length L as the set of k-mers contained within that sequence. In an ideal random DNA sequence, the expected number of occurrences of a specific "canonical k-mer" (i.e., the k-mer itself and the one with the smaller hash value in its inverse complementary strand) is calculated. ,as follows:

[0080]

[0081] in, Represents the total length of the sequence. Represents the length of k-mer, Represents the size of the alphabet (for DNA, ), This represents the total number of k-mers in the sequence.

[0082] This invention recognizes that the k-mer diversity of the real genome is significantly lower than that of a random sequence of the same length due to the presence of numerous repetitive sequences and low-complexity regions. To correct this bias, this invention introduces the key concept of "effective sequence size," using a correction factor... To compensate for this decrease in k-mer diversity, since dynamically calculating the number of k-mers for each dataset at runtime is considered "infeasible," this invention employs a preset constant correction factor. Corrected expected number of occurrences Calculated using the following formula:

[0083]

[0084] in: The effective sequence size correction coefficient is given.

[0085] Based on this correction The probability of a single k-mer match producing a "false match" Defined by the following formula:

[0086]

[0087] The false positive rate of the IBF data structure due to hash collisions is essentially set as an acceptable target for the user to calculate memory overhead. This is the false positive rate of the IBF data structure due to hash collisions when calculating memory overhead. The calculation is as follows:

[0088]

[0089] in, The number of hash functions, i.e., the number of hash functions required to insert or retrieve an element. A number of different hash values. : The number of elements to be inserted, i.e., the total number of k-mers expected to be stored in a single bin of IBF. : Total number of bits in the filter, i.e., the memory size occupied by a single bin of the IBF. The base of the natural logarithm (approximately 2.71828). The user sets an acceptable target. (For example, 0.00001, which represents a false positive rate of 0.001%). The system uses the above formula to calculate how much memory it needs to allocate to a single bin of IBF to achieve the user's desired performance. value.

[0090] and The combined total false positive probability of a single k-mer is obtained .

[0091] Finally, the system uses the cumulative probability function of the binomial distribution to accurately calculate a probability distribution containing... In k-mer query patterns, there is occasional sharing of more than The probability of k-mer is the final estimate of FPR. .Should The calculation is based on the following formula:

[0092]

[0093] in, Represents the total number of k-mers in this pattern (i.e. ). This represents the shared k-mer count threshold required for filtering. This represents what is observed in the pattern. The cumulative probability of one or fewer k-mer hits. Representative observed The probability of hitting k-mer multiple times is the final estimate of the FPR.

[0094] Step 14: Optimization Solution and Parameter Output; After determining the calculation models for FNR and FPR described in Steps 12 and 13, the system operates within a preset parameter space (e.g., From 7 to 23, Perform an efficient grid search from 1 to 25. At each grid point... The above uses the formulas for FNR and FPR to calculate the cost respectively. Finally, through global optimization, a method is found that minimizes the cost function. The parameter pair that obtains the minimum value These parameters will be established as the optimal execution parameters for this comparison task and passed to subsequent stages.

[0095] After determining the theoretically optimal k-mer length through grid search in step 1-4, the system instantiates this abstract length value into a concrete k-mer pattern for performing the alignment. The system has a predefined, highly sensitive pattern library containing a series of proven efficient spaced k-mer patterns (including symmetric spaced k-mer patterns). A spaced k-mer is a sequence seed whose shape is defined by a binary string of 1s (must match) and 0s (ignored / spaced). The effective length of the spaced k-mer is the number of 1s in the corresponding binary string. These patterns correspond one-to-one with their respective effective length values. The system selects a pattern from this library whose effective length matches the optimal length derived in step 1-4. The values ​​match exactly.

[0096] Step 2: Building the Parallel Computing Architecture and Data Preprocessing:

[0097] Step 21: Selection and configuration of core seed strategy.

[0098] The selection strategy in the morphology library of this invention preferably adopts a symmetrical spaced-interval k-mer morphology, meaning that the binary string corresponding to the selected spaced-interval k-mer is symmetrical. This strategy is to address the technical problems that arise when combined with the canonical k-mer strategy used in this invention. The canonical k-mer strategy is a standard method for efficiently processing DNA double strands, which selects only the k-mer with the smaller hash value from either the forward strand or its anticomplementary strand. When using an asymmetric seed, the canonical k-mer strategy can cause the number of physically overlapping bases between adjacent k-mers to become inconsistent and potentially increase significantly, thereby reducing the sensitivity of the filter. Therefore, this invention adopts a symmetrical spaced-interval k-mer morphology (i.e., its binary representation is centrosymmetric). This configuration ensures that regardless of whether the canonical k-mer is selected from the forward or anticomplementary strand, the number of physically overlapping bases with adjacent canonical k-mers remains constant and minimal. This configuration ensures that the k-mer filter has the highest and most predictable sensitivity and robustness when edit errors exist in the sequence.

[0099] Step 22: Implementation of sequence segmentation and overlap strategies.

[0100] To achieve parallelization, the reference genome (e.g., human) and the query genome (e.g., mouse) are respectively segmented into appropriately sized fragments (e.g., 10 Mb) with fixed-length overlapping regions (e.g., 1 Mb), such as... Figure 3 As shown on the left. The design of the overlapping region is key to ensuring the integrity of the results. It ensures that any true local alignment that crosses the fragment boundary can be fully captured by at least one parallel task, avoiding the boundary effect problem introduced by data segmentation.

[0101] Steps 2 and 3: Construction of the global IBF index and quantization of repeating sequence regions.

[0102] a) Construct a global IBF index: Logically divide all reference fragments into different "bins" (e.g., ...). Figure 3 The algorithm consists of four bins, each serving as a reference bin. A globally shared Interleaved Bloom Filter (IBF) index is built for all bins. "Global" means that this IBF index is built for all reference fragments. It is a single master index containing information for all reference fragments, rather than an index targeting only a specific region of the genome. "Shared" means that all parallel producer threads access the same IBF index to determine whether a k-mer in the query sequence exists in a particular bin. The IBF index is a highly space-efficient probabilistic data structure that can determine which bins a k-mer might exist in with minimal memory usage and extremely high speed.

[0103] b) Quantitative repeatability: While constructing the global IBF index, the variability score of the reference sequence fragment in each bin is calculated. It is defined as the number of unique k-mers within the reference sequence segment i. Total number of k-mers within the fragment The ratio of . The closer the value is to 1, the more complex the fragment and the lower the repetition rate. The system uses the average variability score of all reference fragments contained in a certain bin as the overall variability score of that bin. The overall variability score of that bin will be used for subsequent intelligent repetitive sequence processing.

[0104] Step 3: Execute hierarchical filtering and comparison based on the "producer-consumer" model:

[0105] Step 3: 1. Producer thread and first-level global filtering.

[0106] The system launches multiple producer threads based on the number of available CPU cores. Each producer thread is responsible for processing a portion of the query fragments. For each fragment, the producer extracts the k-mer and queries the global IBF index, accumulating the shared k-mer hit count between each "query fragment-reference warehouse" pair. If this number exceeds a threshold... If so, the "query fragment - reference repository" pair is considered a candidate task.

[0107] Step 3.2: Dynamic priority scheduling of repeating sequences.

[0108] When determining candidate tasks, if a query segment Q hits more than half of the total number of warehouses, the system determines that query segment Q has high repetition. For example... Figure 3 As shown in the middle section, there are query fragments that simultaneously hit warehouses 1, 3, and 4, with a total of 3 hit warehouses, which is more than half of the total number of warehouses (4). For such query fragments, the system will not simply discard them, but will prioritize pairing the query fragment Q with the reference warehouse (e.g., warehouse 4) with the highest overall variability among all hit warehouses to form a candidate task based on the variability score calculated in steps two and three. This intelligently guides computing resources to the region most likely to generate specific alignments.

[0109] Step 33: Decouple the asynchronous task queue from the system.

[0110] Producers package selected and prioritized candidate tasks into a cart data structure and push them to a global, thread-safe, full-load cart queue. This queue acts as a buffer between producers and consumers, achieving asynchronous decoupling between them, and is key to the system's high throughput and load balancing.

[0111] Steps 3 and 4: Consumer thread and second-level local fine-grained filtering.

[0112] Multiple consumer threads independently retrieve carts from the queue. For each cart's task, the consumer first constructs a local and temporary full-text index I in its thread-private memory, containing only the precise positions of the k-mers required for the current task's query segment. "Local" means that the scope of index I is limited to the batch of query segments contained in the current cart, not all query segments; "temporary" refers to the lifespan of index I, meaning it is constructed in real-time within the consumer thread's private memory. Once the task in this cart is processed, index I is discarded, the memory is released, and "full-text" means that index I records the precise position of each k-mer in the query segment. Subsequently, the SWIFT filter algorithm is executed using this index I to perform a lossless secondary filtering by checking for collinearity of shared k-mers.

[0113] Step 35: Final precise comparison.

[0114] Only after passing the second-level SWIFT filter and identifying a very small number of highly reliable candidate region pairs will the consumer thread invoke the optimized striped Waterman-Eggert dynamic programming algorithm to perform the final base-by-base precise local alignment, such as... Figure 3 As shown on the right. Calculate the comparison score, start and end coordinates, and output the results.

[0115] Step 4: Global integration and output of parallel computing results:

[0116] Step 41: Result Collection and Merging; One main thread is responsible for asynchronously collecting the partial comparison results output by all consumer threads, such as... Figure 2 As shown.

[0117] Step 4.2: Deduplication; Due to overlap in the initial segments, a truly long alignment may be repeatedly discovered by multiple parallel tasks. This step uses an efficient sorting and interval merging algorithm to accurately identify and merge these duplicate alignments based on the genomic coordinates of the alignment, ensuring the uniqueness of the final output.

[0118] 1. Experimental setup:

[0119] Benchmark comparison method: BLAST (v2.12.0), LAST (v1595), and LASTZ (v1.04.22) were selected for comparison.

[0120] Multi-dimensional performance evaluation metrics: core performance metrics (runtime, peak memory usage) and accuracy / sensitivity metrics (error rate).

[0121] To verify the beneficial effects of the method, this invention used two sets of simulated data and one set of real data for evaluation. The experiment was conducted on a high-performance computing server equipped with dual AMD EPYC 9454 48-core processors and 1007GB of available memory. All tools were run at default parameters to ensure the objectivity and reproducibility of the comparison.

[0122] Reference sequence alignment dataset: This invention uses two sets of simulated data, 50MB and 250MB in size, as well as a set of real human (GRCh38) and mouse (GRCm39) genome data.

[0123] Benchmark Comparison Method

[0124] To comprehensively evaluate the performance of this invention, it selects several advanced prior art technologies that are representative in the field of sequence alignment as benchmarks for comparison, including BLAST (v2.12.0), LAST (v1595), and LASTZ (v1.04.22).

[0125] 2. Performance evaluation indicators:

[0126] 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, error rate, and resource consumption.

[0127] Sensitivity: Sensitivity represents the proportion of all real-world variants that the model correctly detects. This metric assesses the model's ability to find real-world variants.

[0128]

[0129] This represents the number of true variants that were correctly detected. This indicates the number of non-variable sites reported in error reports. This indicates the true number of variants that the model failed to detect.

[0130] Miss Rate: The miss rate directly measures the proportion of positive samples missed by the model.

[0131]

[0132] The computational resource consumption indicators include:

[0133] Running time: Records the total time required from the start of the program to the completion of the task.

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

[0135] 3. Experimental Results and Analysis:

[0136] Experimental Results: The present invention evaluated the performance of the method in three typical application scenarios and compared it with the prior art. The experimental results are shown in Tables 1, 2 and 3.

[0137] Table 1. Comparison of local alignment performance between two 50Mb simulated sequences

[0138]

[0139] Under the stringent condition of an 8% error rate, the default BLAST missed nearly one-third of the true alignments, demonstrating its insufficient sensitivity. Other tools that can guarantee sensitivity, such as LAST and LASTZ, have runtimes measured in minutes. In contrast, this invention completes 100% accurate alignments in just 4 seconds, with a memory footprint of only 550Mb, demonstrating outstanding overall performance.

[0140] Table 2. Comparison of local alignment performance between two 250Mb simulated sequences

[0141]

[0142] When the data size increases fivefold, the performance gap is further amplified. This invention completes a 100% accurate comparison in just 63 seconds.

[0143] Table 3. Local alignments between the human (GRCh38) and mouse (GRCm39) genomes.

[0144]

[0145] In the most challenging human-to-mouse whole-genome alignment task, this invention completed the task in 9 minutes and 13 seconds, significantly faster than BLAST and LAST, while other precise methods were completely incapable of achieving the same result. More importantly, there is a difference in the quality of the result set: this invention output 143,000 high-quality matches that strictly conform to the ε-match criterion; while BLAST and LAST produced 44 million and 251 million matches, respectively, the vast majority of which were short, high-error-rate "noisy" alignments.

[0146] The experimental results above demonstrate that this invention not only achieves a breakthrough in speed, but also provides a novel analytical paradigm that can directly extract high-quality, high-confidence biological signals from massive amounts of data.

[0147] 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 parallel precise local sequence alignment method based on hierarchical filtering and adaptive parameter optimization, characterized in that, include: S100. Determine ε-match based on the comparison criteria set by the user and the length of the input sequence, and determine the combination cost based on the false positive rate (FPR) and false negative rate (FNR). By minimizing the combinatorial cost, a set of theoretically optimal k-mer filtering parameters for the current task is determined, including the k-mer length k, shape, and shared k-mer counting threshold. ; S300, two-level hierarchical filtering, and parallel distribution of "producers-consumers": Global coarse-grained filtering: executed in parallel by multiple producer threads; all reference fragments are logically divided into different "warehouses," each warehouse being a reference warehouse, and a globally shared IBF index is built for all warehouses; all parallel producer threads access the same IBF index to determine whether the kmer in the query sequence exists in a certain warehouse, based on a threshold. Filter out some combinations of "query fragment - reference repository"; Task Packaging and Distribution: The producer thread will package the "query fragment - reference warehouse" tasks through the first layer of filtering into a data structure called "cart". The cart will be placed into a global, first-in-first-out task queue. Local fine-grained filtering and precise comparison: This is executed in parallel by multiple consumer threads. Each consumer thread retrieves a "cart" from the task queue and, for all query segments contained within the "cart," immediately builds a local full-text index in memory that records the precise position of the k-mer in the query segment. Using the full-text index, a SWIFT filter is executed. After filtering, the results are fed into the Waterman-Eggert dynamic programming algorithm to find all local comparisons between two sequences that exceed a threshold. S400: Collect the local alignment results output by all consumer threads, identify and merge duplicate alignments based on the aligned genomic coordinates, and finally output the results.

2. The parallel precise local sequence alignment method based on hierarchical filtering and adaptive parameter optimization according to claim 1, characterized in that, The comparison criterion is that the comparison length must be greater than or equal to the minimum comparison length. And the number of edit operations compared must be less than or equal to the maximum number of edit operations. .

3. The parallel precise local sequence alignment method based on hierarchical filtering and adaptive parameter optimization according to claim 2, characterized in that, The number of edit operands must be less than or equal to the maximum number of edit operands. The judgment method is as follows: For a length of ,Include The comparison of each editing operation is performed by maximizing the normalization of editing operation conditions. To make a judgment.

4. The parallel precise local sequence alignment method based on hierarchical filtering and adaptive parameter optimization according to claim 1, characterized in that, The false negative rate (FNR) was obtained using a three-dimensional dynamic programming algorithm, including: A DP table, which is a three-dimensional dynamic programming matrix, is pre-computed and populated to store the number of false negative error configurations. The value of FNR is determined by querying this pre-computed DP table during the process of optimal k-mer filtering parameters. DP table status middle Represents the length of the sequence. Represents the total number of editing operations. Represents the length of k-mer, Represents the shared k-mer count threshold required for filtering; state as follows: The total number of error distributions is determined by the number of combinations. And thus obtain .

5. The parallel precise local sequence alignment method based on hierarchical filtering and adaptive parameter optimization according to claim 1, characterized in that, The false positive rate (FPR) is obtained through the following steps: First, a sequence of length L is modeled as a set of k-mers contained in the sequence; The expected number of occurrences of a specific "canonical k-mer" as follows: in: The effective sequence size correction coefficient is mentioned above; Represents the size of the alphabet. This represents the total number of k-mers in the sequence; The probability of a single k-mer expression producing a "false match". : False positive rate due to hash collisions in IBF data structures and The total false positive probability of a single k-mer is obtained. Thus, the probability is obtained. : in, This represents the total number of k-mers in the pattern; This represents what is observed in the pattern. The cumulative probability of one or fewer k-mer hits; Representative observed The probability of hitting k-mer multiple times is the final estimate of the FPR.

6. A parallel precise local sequence alignment method based on hierarchical filtering and adaptive parameter optimization according to any one of claims 1 to 5, characterized in that, During the task packaging and distribution process, dynamic priority scheduling of the repeating sequence is first performed, and then the producer packages the priority-scheduled candidate tasks into the cart. The process of dynamic priority scheduling for repeating sequences includes: When determining candidate tasks, if the number of hit bins for a query fragment Q exceeds half of the total number of bins, the system determines that the query fragment Q is highly repetitive. For such queries, the system prioritizes pairing the query fragment Q with the reference bin with the highest overall variability among all hit bins to form candidate tasks, thereby directing computational resources to the region most likely to generate specific alignments.

7. The parallel precise local sequence alignment method based on hierarchical filtering and adaptive parameter optimization according to claim 6, characterized in that, The process for calculating the overall variability of the warehouse includes: While constructing the global IBF index, the variability score of the reference sequence fragments in each bin is calculated. It is defined as the number of unique k-mers within the reference sequence segment i. Total number of k-mers within the fragment The ratio; the overall variability score obtained by averaging the variability scores of all reference segments contained in a given warehouse.

8. The parallel precise local sequence alignment method based on hierarchical filtering and adaptive parameter optimization according to claim 7, characterized in that, During the global coarse-grained filtering process, multiple producer threads are started based on the number of available CPU cores. Each producer thread is responsible for processing a portion of the query fragments. For each fragment, the producer extracts the k-mer and queries the global IBF index, accumulating the shared k-mer hit count between each "query fragment-reference warehouse" pair. If this number exceeds the threshold t determined by S100, the "query fragment-reference warehouse" pair is considered a candidate task.

9. The parallel precise local sequence alignment method based on hierarchical filtering and adaptive parameter optimization according to claim 7, characterized in that, The construction method for a full-text index is as follows: Multiple consumer threads independently retrieve carts from the queue. For each task in a cart, the consumer first constructs a local and temporary full-text index I in its thread-private memory, which contains only the k-mer precise position of the query segment required by the current task. Local means that this index I is limited to the batch of query segments contained in the current cart. Index I is constructed in the private memory of the consumer thread in real time. When the task in the cart is finished, index I is discarded and the memory is released.

10. A parallel precise local sequence alignment method based on hierarchical filtering and adaptive parameter optimization according to claim 7, characterized in that, In the process of two-level hierarchical filtering and parallel distribution of "producers-consumers", a symmetrical k-mer pattern with added intervals is adopted.

Images